Marc Deschamps

Michel Castings

Mathieu Renier

Olivier Casula

Karim Jezzine

Keywords: Ultrasonic waves, Non Destructive Testing, Topological Imaging, Waveguide, Structural Health Monitoring

Tubular structures are widely used in a variety of industries. Non Destructive Evaluation (NDE) of these structures plays a crucial role during it's life cycle. In order to test large structures with limited accessibility, guided wave testing was developed as a viable solution. Due to the nature of these waves, they are able to propagate over large distances without losing much of their energy. However, they are also complex in that their velocity is frequency dependent i.e. they are dispersive. Conventionally, guided wave testing require costly finite element simulations. This thesis offers an alternative to such simulations with a quick and robust method to simulate guided wave propagation in tubular structures.

Based on these calculations, the aim of this work is to obtain the 3d topological image of multilayered isotropic tubular structures using ultrasonic guided waves to locate defects.

A mathematical model has been proposed where the wave equation is converted to an ordinary differential equation with respect to radius 'r' using the Fourier and Laplace transforms for the spatial and temporal variables respectively. The partial wave solution, expressed as a combination of Bessel's functions, allows for the creation of a fast robust semianalytical algorithm to compute the Green function in tubular structures. A model to approximate numerical defects is then developed. The defect response is considered as the cumulative response of secondary sources, aiming to negate the incident and diffracted stress field present within it. Next, the numerical model is validated with experimental measurements. Finally, the technique of Topological Imaging is introduced. This method of imaging is based on the idea of performing a correlation between two wave fields for defect localization. The versatility and flexibility of the numerical tool in conjunction with the method of imaging is then successfully demonstrated by localising and imaging a multitude of numerical and experimental defects with dimensions as low as 1=40th of the wavelength.

Résumé

Les structures tubulaires sont largement utilisées dans diverses industries. Le Contrôle Non Destructif (CND) de ces structures joue un rôle crucial au cours de leur cycle de vie. Afin de tester de grandes structures avec une accessibilité limitée, la méthode de CND utilisant des ondes guidées a été développée comme une solution viable. En raison de la nature de ces ondes, elles sont capables de se propager sur de grandes distances sans perdre une grande partie de leur énergie. Cependant, elles sont complexes puisque leur vitesse dépend de la fréquence, c'est-à-dire qu'elles sont dispersives. Classiquement, l'étude de ce type d'ondes nécessite des simulations par éléments finis coûteuses. Cette thèse propose une alternative à de telles simulations avec une méthode rapide et robuste pour simuler la propagation d'ondes guidées dans des structures tubulaires. Partant de ces calculs, pour localiser des défauts, l'objectif de ce travail est d'obtenir des images topologiques 3D de structures tubulaires isotropes multicouches par propagation de ces ondes guidées ultrasonores. Un modèle mathématique est proposé où l'équation d'onde est convertie en une équation différentielle ordinaire par rapport au rayon «r» en utilisant les transformées de Fourier et de Laplace pour les variables spatiales et temporelles respectivement. La solution en ondes partielles, exprimée comme une combinaison des fonctions de Bessel, permet la création d'un algorithme semi-analytique rapide et robuste pour calculer la fonction de Green de structures tubulaires. Un modèle approché en présence de défauts numériques est ensuite développé. La réponse des défauts est considérée comme la réponse cumulative des sources secondaires, visant à annuler le champ de contraintes incident et diffracté présent en son sein. Ensuite, le modèle numérique est validé par des mesures expérimentales. Enfin, la technique de l'imagerie topologique est introduite. Cette méthode d'imagerie est basée sur la corrélation entre les champs ultrasonores sans et avec défaut. La polyvalence et la flexibilité de l'outil numérique en conjonction avec cette méthode d'imagerie sont ensuite démontrées avec succès en localisant et imageant une multitude de défauts numériques et expérimentaux avec des dimensions aussi faibles que 1/40e de la longueur d'onde. The very first versions of sonic NDT was performed by blacksmiths to determine the integrity of bells. Listening to it ring, an experienced blacksmith could hammer it into the right shape or figure out if it was damaged. Developments in the field led to different methods ranging from X-Ray, Magnetic Particle to Radiography testing. However, most of these methods require expensive equipment which are generally non portable or are only able to test small sections of the specimen. Ultrasonic testing emerged in the 1930s as an idea to detect flaws within structures. This method proved to be useful to test large structures due to the high penetration power of sound waves through solids. Guided wave ultrasonic testing technique quickly developed to be used extensively in long range pipeline testing. Guided waves are waves that propagate in thin structures (thickness in the order of wavelength) guided by the boundaries of the structure allowing it to propagate over long distances. These properties of guided waves allow for the detection of inaccessible and distant flaws in structures with just single fixed probe position.

One of the earliest work done with regard to guided wave propagation in hollow cylindrical structures was by D.C Gazis in 1959. His papers describe the analytical and numerical foundations for wave equations in hollow cylinders using Helmoltz potentials [START_REF] Gazis | Three-Dimensional Investigation of the Propagation of Waves in Hollow Circular Cylinders. I. Analytical Foundation[END_REF][START_REF] Gazis | Three-Dimensional Investigation of the Propagation of Waves in Hollow Circular Cylinders. II. Numerical Results[END_REF]. The earliest mention of numerical simulations of guided waves was by Zemanek [START_REF] Zemanek | An Experimental and Theoretical Investigation of Elastic Wave Propagation in a Cylinder[END_REF] who in 1972 performed the computation and verified them with experimental results.

A paper by Rose [START_REF] Rose | Ultrasonic guided waves in structural health monitoring[END_REF] discusses some aspects and applications of guided wave modes. Mohr and Höller [START_REF] Mohr | On Inspection of Thin-Walled Tubes for Transverse and Longitudinal Flaws by Guided Ultrasonic Waves[END_REF] described using longitudinal and torsional modes for defect detection.

To simplify the selection of parameters to be used for ultrasonic guided wave testing, the program "Disperse" (a software to obtain dispersion curves for a vast variety of conditions and materials) was developed by Lowe [53] and has been expanded to cylinders by Pavlakovic [START_REF] Pavlakovic | Review of Progress in Quantitative Nondestructive Evaluation[END_REF].

Ensuing work by Demma [START_REF] Demma | The interaction of guided waves with discontinuities in structures[END_REF][START_REF] Demma | The reflection of guided waves from notches in pipes: A guide for interpreting corrosion measurements[END_REF], Alleyne [START_REF] Alleyne | The Effect of Discontinuities on the Long-range Propagation of Lamb Waves in Pipes[END_REF], Lowe [START_REF] Lowe | Defect detection in pipes using guided waves[END_REF] etc., deal with the interaction and responses of guided waves of various modes with discontinuities. Later works by Cawley and Alleyne [START_REF] Cawley | Practical Long Range Guided Wave Inspection -Managing Complexity[END_REF] also detail the reasons for using axi-symmetric modes in guided wave inspection of pipes. This is due to the ease of excitation of a single pure mode, constant sensitivity in all circumferential positions and to control coherent noise. There have also been discussions on using helical guided modes by Willey et al. [START_REF] Willey | Guided wave tomography of pipes with high-order helical modes[END_REF] and Dehghan [START_REF] Dehghan-Niri | A multi-helical ultrasonic imaging approach for the structural health monitoring of cylindrical structures[END_REF] The fourth chapter introduces the principle of Topological Imaging. This is a technique for identifying and localizing defects using the cross correlation of two fields. The technique is then applied on numerical and experimental defects. Some studies with numerical defects tie into the second chapter to better understand the functioning of the defect model. Finally experimental defects are localized using the principles and techniques developed within this thesis.

Chapter 1

Transient Wave Propagation in a

Stratified Tubular Medium

The complex dynamics of wave propagation has been studied for over a century. One of, if not the earliest, successful attempts to understand and model this was by Rayleigh [START_REF] Rayleigh | On Waves Propagated along the Plane Surface of an Elastic Solid[END_REF] in 1885, where he investigated the behaviour of surface waves propagating along semi-infinite half surfaces. In 1917, Lamb [START_REF] Lamb | On Waves in an Elastic Plate[END_REF] introduced the idea of wave propagation in a single flat layer of finite thickness. He discussed symmetric (S 0 ) and anti-symmetric (A 0 ) modes in plates which are now commonly known as Lamb waves. It was later shown by Love [START_REF]Some Problems of Geodynamics[END_REF] that Shear Horizontal (SH) guided waves could propagate in a medium of finite thickness.

Stoneley [START_REF] Stoneley | Elastic Waves at the Surface of Separation of Two Solids[END_REF] in 1924 described the waves propagating along the interface of two elastic solids. The naming conventions of the equivalent of these modes in hollow cylinders as described by Silk and Bainton [START_REF] Silk | The propagation in metal tubing of ultrasonic wave modes equivalent to Lamb waves[END_REF] is discussed in a later section 3.3. The next step was to understand propagation through multilayered media. This was done in a study by Lindsay [67] in 1939 where the propagation of oblique waves was studied through alternate layers of fluid and solid where the media as a whole behaves as a low pass filter. Thomson [START_REF] Thomson | Transmission of elastic waves through a stratified solid medium[END_REF] in 1950 developed the transfer matrix method for multilayer solid media which was later corrected by Haskell [START_REF] Haskell | The Dispersion of surface waves on multilayered media[END_REF] in 1953. The transfer matrix method consists of relating the amplitudes at a lower interface to that of the uppermost interface and this procedure is carried down through the multiple layers in the medium. Most of the aforementioned methods were developed in the context of geophysics and seismology dealing with large thicknesses and wavelengths. The transfer matrix method had a drawback in that it was unstable in cases where the thickness was much larger than the wavelength of the propagating wave. This was the case where computations were to be performed at high frequencies and came to be known as the "large fd problem". This problem was addressed in the early 1960s by Dunkin [START_REF] Dunkin | Computation of modal solutions in layered, elastic media at high frequencies[END_REF] and others [START_REF] Dunkin | Deformation of a layered, elastic, half-space by uniformly moving line loads[END_REF][START_REF] Schwab | Surface-Wave Dispersion Computations. Knopoff'S Method[END_REF][START_REF] Schwab | Fast Surface Wave and Free Mode Computations[END_REF][START_REF] Kundu | Elastic waves in a multilayered solid due to a dislocation source[END_REF][START_REF] Kundu | Acoustic Material Signature[END_REF][START_REF] Kundu | Calculation of the acoustic material signature of a layered solid[END_REF][START_REF] Kennett | Seismic Wave Propagation In Stratified Media[END_REF], one of them being Knopoff [START_REF] Knopoff | A matrix method for elastic wave problems[END_REF] who in 1964 introduced the Global Matrix method. This method involved assembling the responses of all layers as a single large matrix by applying continuity conditions across the boundaries of the layers. It has the advantage being robust and reliable for variety of solutions and has been implemented ever since by a number of authors [START_REF] Randall | Fast programs for layered half-space problems[END_REF][START_REF] Chin | Matrix methods in synthetic seismograms[END_REF][START_REF] Schmidt | Efficient numerical solution technique for wave propagation in horizontally stratified environments[END_REF][START_REF] Schmidt | Efficient global matrix approach to the computation of synthetic seismograms[END_REF][START_REF]Wave propagation in layered composite laminates under periodic surface loads[END_REF]. A particular disadvantage is that with the addition of more layers to the system, the matrix becomes bulky and computationally expensive. A review of the matrix methods by Lowe [START_REF] Lowe | Matrix techniques for modeling ultrasonic waves in multilayered media[END_REF] gives an extensive overview most of the work done with respect to the Transfer and Global matrix methods since the early stages. The Global Matrix method has been adopted in this work.

The following sections details the mathematical framework developed to compute the transient field in multi-layered transversely isotropic tubular structures.

In the approach followed, the spatial and temporal variables of the guided wave equations in cylinders are expressed in the Fourier and Laplace domains respectively for a given radial position within a layer of the structure. This follows from the work done by Mora et al. [START_REF] Mora | Interaction of a guided wave with a crack in an embedded multilayered anisotropic plate: Global matrix with Laplace transform formalism[END_REF][START_REF] Mora | Transient 3D elastodynamic field in an embedded multilayered anisotropic plate[END_REF] where the equations have been resolved in the transformed domain for guided wave propagation in layered plates. In case of the cylinder, the solution of the equations are expressed as a combination of modified Bessel functions to give six partial waves for each layer of the structure, considering the layers to be unbounded. The layers are then assembled using the Global Matrix method taking into account the continuity of stresses and displacements across the layers. This matrix describes the behaviour of the propagation of guided waves within the assembly and has been details below.

Notations and domains for calculation

In the physical domain, the cylindrical coordinates radius r, azimuth θ and axial position z are defined such that --→ OM = r n r + z n z , where the unit vector n r depends on the azimuth θ as shown in Figure 1.1. Time is denoted as t.

The calculation of the transient field is performed for a given radial position r. This field is calculated in the (n, k, s) domain. The axial wavenumber k corresponds to the variable of the Fourier transform F with respect to the axial position z, given by:

F(u)(k) = ∞ -∞ u(z) exp(i k z) dz ⇐⇒ u(z) = 1 2 π ∞ -∞ F(u)(k) exp(-i k z) dk. (1.1)
The azimuthal or circumferential wavenumber is an integer n, and corresponds to the variable of the Fourier series S with respect to the azimuth θ. A Fourier series expansion is used as the domain is periodic about the circumference and is given by:

S(u)(n) = 1 2 π π -π u(θ) exp(i n θ) dθ ⇐⇒ u(θ) = ∞ n=-∞ S(u)(n) exp(-i n θ). (1.2)
The variable s is the complex variable of the Laplace transform L with respect to time t, given ∀a > 0 by:

L(u)(s) = ∞ 0 u(t) exp(-s t) dt ⇐⇒ u(t) = exp(a t) 2 π ∞ -∞ L(u)(a + i ω) exp(i ω t) dω Bromwich-Mellin Formula
.

(1.3)

Performing the computations in the Laplace domain rather than in the Fourier domain allows to sidestep certain issues caused by aliasing. This aliasing may occur due to the presence of slow or extremely dispersive waves in the frequency bandwidth of computation.

These waves may persist beyond the duration of the computation. Using the complex Laplace variable allows for truncated causal computations which could potentially reduce the time t required for a stable computation. This behaviour has been detailed by Phinney [START_REF] Phinney | Theoretical calculation of the spectrum of first arrivals in layered elastic mediums[END_REF] where he has described it as an imperfectly trapped mode.

It is to be noted that the modeling of the defect would require for the equations to be resolved in a mixed domain of real spatial variables with the Laplace temporal variable. The domains of calculation as well as their respective variables are as described in

Modelling of wave propagation equations

Consider a multilayered medium of N perfect transversely isotropic cylindrical layers stacked together as shown in Figure 1.2. A reference layer is named β with the interface between the layer β and β + 1 also named β located at a radial position given by r β . The structure is assumed to be infinite in the z direction and can have vacuum, isotropic solids or fluids on the inside and outside. The external force can be located anywhere within the region of computation.

Basic equations in the physical space

The mechanical properties, i.e., mass density ρ β and stiffnesses c β ijkm are constant in each layer. With the aforementioned notations the wave equation is written by combining the principle of causality, Newton's second law and Hooke's law.

Newton's second law is expressed as follows:

ρ β ∂ 2 t u - 1 r + ∂ r σ r + 1 r [T + ∂ θ ] σ θ + ∂ z σ z = f β , (1.4)
where u denotes the displacement vector, σ r , σ θ and σ z are the stress vectors in the radial, azimuthal and axial directions, respectively, the field f β denotes the force per unit volume exerted by the part of the source located in layer β, and

T =    0 -1 0 1 0 0 0 0 0    .
Hooke's law expresses the stress σ d in direction d (unit vector) with respect to displacement u as:

σ d = (d β ∇)u + 1 r (d β n θ ) T u , (1.5)
where ∇ is the gradient operator ∂ r , 1 r ∂ θ , ∂ z T and the bilinear product β has been defined by Ducasse in [START_REF] Ducasse | A nonstandard wave decomposition to ensure the convergence of Debye series for modeling wave propagation in an immersed anisotropic elastic plate[END_REF] by a three-by-three matrix (a β b) such that (a β b) im = c β ijkm a j b k , following the Einstein summation convention. This operator has the advantage of being able to concisely represent material properties even for complex cases such as anisotropy.

Consequently, the displacement field u(r, θ, z, t), at any time t and any location (r, θ, z), satisfies the following wave equation expressed in the β layer:

ρ β ∂ 2 t u(r, θ, z, t) - + 1 r n r β + + 1 r n r β n θ T r + T r (n θ β ) + (n θ β n θ ) T r u(r, θ, z, t) = f β (r, θ, z, t),
for t > 0, u(r, θ, z, t) = 0, for t < 0.

(

As we will need to write the continuity of displacement and radial stress at each interface, the radial stress has to be expressed as follows:

σ r (r, θ, z, t) = (n r β ) + (n r β n θ )
T r u(r, θ, z, t) .

(1.7)

The above equations depend on each layer through the values of the elastic constants c β ijkm , i.e. through the operator β .

The stiffness tensor can be represented by a 6-by-6 symmetric matrix, using the Voigt notation:

1 ↔ rr | 2 ↔ θθ | 3 ↔ zz | 4 ↔ θz | 5 ↔ rz | 6 ↔ rθ .
For a transversely isotropic medium, the stiffness tensor is:

             
c 44 0 c 66               , (1.8)
while for an isotropic medium: c 33 = c 11 , c 13 = c 11 -2 c 66 and c 44 = c 66 .

Continuity equations at the interface β are given by the following equations:

  u(r + β , θ, z, t) σ r (r + β , θ, z, t)   -   u(r - β , θ, z, t) σ r (r - β , θ, z, t)   =   0 p β (θ, z, t)   , (1.9)
where r - β and r + β indicate the fact that the field under consideration is calculated in layers β-1 and β, respectively. The interface source term p β (θ, z, t) defines the normal (radial) stress jump at interface β, and corresponds to an applied force per unit area. The condition p β = 0 merely expresses the continuity of displacement and normal stress. It is to be noted that if one of the layers is vacuum (as is the case for the outermost layer of the cylinder in vacuum), the source term (as described in Section 1.4.3) would represent a surfacic source that is used to represent a transducer loading at the surface.

The index β is omitted and will be reintroduced only when necessary to avoid ambiguity. Eq. (1.6) and (1.9) are to be solved using Fourier and Laplace transforms in the invariant dimensions (θ, z) and (t), respectively. In cases with volumic sources, Eq. (1.6) is solved separately in each layer containing the source term. This would define an incident field within each layer, which corresponds to the field that a source radiates into this layer considering it to be unbounded. However, this thesis only deals with sources located at interfaces and hence the volumic source terms are set to 0. In case of the structure placed in vacuum, the outer and inner surfaces are also considered as interfaces thus allowing for the simulation of surface sources. The refracted field is then obtained by the contribution of all the interfaces, and is calculated by taking into account the continuity relationships of Eq. (1.9).

1.3 Equation to be solved in the (r, n, k, s)-domain

General case

After a Fourier transform with respect to the axial position z, a Fourier series expansion with respect to the angular position θ, and a Laplace transform with respect to time t, the wave equation (1.6) becomes the following ordinary differential system in the (r, n, k, s)domain (see Table 1.1):

(n r n r ) Ũ (r) + -i k (n r n z + n z n r ) + 1 r [ (n r n r ) - i ( (n r n θ ) (n I + i T) + (n I + i T) (n θ n r ) ) ] Ũ (r) - ρ s 2 I + k 2 (n z n z ) + k r [ i (n r n z ) + (n I + i T) (n θ n z ) + (n z n θ ) (n I + i T) ] + 1 r 2 (n I + i T) (n θ n θ ) (n I + i T) Ũ(r) = -F(r) , (1.10)
where I denotes the 3-by-3 identity matrix.

The radial stress in the (r, n, k, s)-domain is given by:

Σr (r) = (n r n r ) Ũ (r) -i 1 r (n r n θ ) (n I + i T) + k (n r n z ) Ũ (r) . (1.11)
Equation (1.10) has no analytic solution for the general case. However, an analytical solution is known for transversely isotropic medium and is detailed below. For transversely isotropic media, Equations (1.10) and (1.11) can be simplified by using the properties summarized in Table 1.2 to obtain the following differential equation:

Transversely isotropic medium

Ũ (r) + C 1 (r) Ũ (r) -C 0 (r) Ũ(r) =          -Fr (r) c 11 -Fθ (r) c 66 -Fz (r) c 44         
, where :

C 1 (r) =          1 r -i n (c 11 -c 66 ) r c 11 -i k (c 13 + c 44 ) c 11 -i n (c 11 -c 66 ) r c 66 1 r 0 -i k (c 13 + c 44 ) c 44 0 1 r          and C 0 (r) =          ρ s 2 c 11 + 1 r 2 1 + c 66 c 11 n 2 + c 44 c 11 k 2 -i n (c 11 + c 66 ) r 2 c 11 0 i n (c 11 + c 66 ) r 2 c 66 ρ s 2 c 66 + 1 r 2 1 + c 11 c 66 n 2 + c 44 c 66 k 2 n k (c 13 + c 44 ) r c 66 i k (c 13 + c 44 ) r c 44 n k (c 13 + c 44 ) r c 44 ρ s 2 c 44 + n 2 r 2 + c 33 c 44 k 2          , (1.12)
and the radial stress:

Σr (r) =      c 11 0 0 0 c 66 0 0 0 c 44      Ũ (r) - i      1 r      i (c 11 -2 c 66 ) n (c 11 -2 c 66 ) 0 n c 66 -i c 66 0 0 0 0      + k      0 0 c 13 0 0 0 c 44 0 0           Ũ(r) .
(1.13)

As mentioned earlier, the refracted field results from the refraction of incident waves and emission of sources at interfaces. This thesis discusses the case with surface sources without any volumic sources for the transversely isotropic case.

The exact solution is given by six partial waves that contain the stress and displacement vectors, three of which are ingoing and three outgoing.

Six partial waves: displacement vectors

The displacement vector Ũj corresponds to six partial waves. The three ingoing waves are expressed as combinations of modified Bessel functions of the first kind I i (see e.g., [81, §10.25]). The first two ingoing waves contain axial displacements and are given as:

For j = 1, 2 , Ũj (r) =         I n-1 (η j r) + I n+1 (η j r) 2
-i I n-1 (η j r) -I n+1 (η j r) 2 i b j I n (η j r)         =         I n (η j r) -i n η j r I n (η j r) i b j I n (η j r)         , (1.14)
where η 2 1 and η 2 2 are the square roots of the following polynomial of the second degree in X (where we take only the positive real part of η j for the computation) :

c 11 X -ρ s 2 + c 44 k 2 c 44 X -ρ s 2 + c 33 k 2 + (c 13 + c 44 ) 2 k 2 X , (1.15)
and the coefficient b j satisfies:

b j = c 11 η 2 j -(ρ s 2 + c 44 k 2 ) -(c 13 + c 44 ) k η j = (c 13 + c 44 ) k η j c 44 η 2 j -(ρ s 2 + c 33 k 2 ) , (1.16)
while the third ingoing wave contains no axial displacement:

Ũ3 (r) =         i I n-1 (η 3 r) -I n+1 (η 3 r) 2 I n-1 (η 3 r) + I n+1 (η 3 r) 2 0         =         i n η 3 r I n (η 3 r) I n (η 3 r) 0         , (1.17)
where

η 3 = ρ s 2 + c 44 k 2 c 66 .
(1.18)

The three outgoing waves are expressed as combinations of modified Bessel functions of the second kind K i :

for j = 1, 2 , Ũj+3 (r) =         -K n-1 (η j r) -K n+1 (η j r) 2 -i -K n-1 (η j r) + K n+1 (η j r) 2 i b j K n (η j r)         =         K n (η j r) -i n η j r K n (η j r) i b j K n (η j r)         , (1.19)
and .20) Note that η 1,2,3 are not exactly radial wavenumbers, even if they have the same unit, because η 2 j = -κ 2 , where κ is a radial wavenumber [see dispersion equations (1.15) and (1.18)].

Ũ6 (r) =         i -K n-1 (η 3 r) + K n+1 (η 3 r) 2 -K n-1 (η 3 r) -K n+1 (η 3 r) 2 0         =         i n η 3 r K n (η 3 r) K n (η 3 r) 0         . ( 1

Six partial waves: radial stress vectors

The radial stress vectors Σr j correspond to six partial waves and are expressed as follows:

for j = 1, 2 , Σr j (r) =         (c 11 η j + c 13 b j k) I n (η j r) + c 66 η j I n-2 (η j r) -2 I n (η j r) + I n+2 (η j r) 2 -i c 66 η j I n-2 (η j r) -I n+2 (η j r) 2 i c 44 (b j η j -k) I n-1 (η j r) + I n+1 (η j r) 2         , (1.21
)

Σr 3 (r) =         i c 66 η 3 I n-2 (η 3 r) -I n+2 (η 3 r) 2
c 66 η 3 I n-2 (η 3 r) + I n+2 (η 3 r) 2 c 44 k I n-1 (η 3 r) -I n+1 (η 3 r) 2         , (1.22)
represent the three ingoing waves and as is the case of the displacement vectors are expressed as a combination of the modified Bessel function of the first kind I i .

Furthermore, the three outgoing waves are represented as a combination of the modified Bessel function of the second kind K i given as:

for j = 1, 2 , Σr j+3 (r) =         (c 11 η j + c 13 b j k) K n (η j r) + c 66 η j K n-2 (η j r) -2 K n (η j r) + K n+2 (η j r) 2 -i c 66 η j K n-2 (η j r) -K n+2 (η j r) 2 i c 44 (b j η j -k) -K n-1 (η j r) -K n+1 (η j r) 2         , (1.23
)

Σr 6 (r) =         i c 66 η 3 K n-2 (η 3 r) -K n+2 (η 3 r) 2
c 66 η 3 K n-2 (η 3 r) + K n+2 (η 3 r) 2 c 44 k -K n-1 (η 3 r) + K n+1 (η 3 r) 2         . (1.24)

Computation of dispersion curves

The formulation described above allows for the computation of the dispersion curves. This could serve as a useful tool for mode selection as well as comparing experimental velocity measurements with theoretical values.

The dispersion curves can be computed from Equations (1.21 -1.24) which describe the radial stress vectors. These stress vectors are assembled after replacing s by i ω and are written as:

Σr (r) = Σr1 (r) Σr2 (r) Σr3 (r) Σr4 (r) Σr5 (r) Σr6 (r)

Software implementation of the model

The mathematical model developed in the earlier section is used to create a simulation tool programmed with Python. First the methodology of discretisation is discussed. This is then followed by the explanation of the implementation of the multiprocessing model.

Then, the numerical source, used for the entirety of the thesis, is defined. Finally, the length of the calculation is obtained based on the dimensions of the computation area.

Discrete calculation (creation of grid)

The continuous physical and temporal domain is discretised for the computation which is performed for a given radial position r. This makes the computation semi-analytic

[analytic in r and numerical in (z, θ, t)]. Numerically, all computations are made using the Fast Fourier Transform as described by Cooley and Tukey [START_REF] Cooley | An Algorithm for the Machine Calculation of Complex Fourier Series[END_REF].

For the axial direction, the real domain variables are given as the maximum axial length z max , discretisation step size dz and number of discretisation points n z = 2z max /dz.

The factor 2 is taken so as to account for the positive and negative axial wavenumber values which is inherent to the Fourier transform, thus leading to a computation for a real domain extending from -z max to z max . In the transformed domain, the axial wavenumber discretisation step size is given as dk z = π/z max and the maximum axial wavenumber is given as k zmax =n z dk z /2. An array of computation points is created as

[0, dk, 2 dk, . . . , k max , -k max +dk, . . . , -dk].
The circumferential wavenumber values are given by discrete integer values due to the periodicity of this physical domain and are obtained in a similar fashion to the axial wavenumber variables. The only difference is that the maximum circumferential wavenumber value is used to obtain the discretisation step size. The array of computation

points are given by [0, 1, 2, . . . , N max , -N max +1, -N max +2, . . . , -1].
The temporal variable t is transformed to the Laplace variable s corresponding to frequency space. The sampling frequency f s satisfies the Nyquist-Shannon criterion (f s > 2 f max ). The discretisation step is dt = 1/f s and the discretised array of s values is [γ, γ + i dω, . . . , γ+i (π f s -dω), γ+i π f s ], where dω = 2 π f s /n t = 2 π/d, the numbers d and n t denoting the duration and the (even) number of time values, respectively. The number of frequency values is n t /2+1 (due to the symmetry of the Fourier transform of a real-valued function) and the positive number γ satisfies exp(-γ d) 1 (10 -5 typically, see [START_REF] Mora | Transient 3D elastodynamic field in an embedded multilayered anisotropic plate[END_REF]).

Multiprocessing

The equations are formulated in a way that allows for multiple independent computations.

The problem is also computationally heavy due to complex calculations performed for each value of N, k and s. For these reasons, a Multiprocessing functionality was decided to be implemented. The computer used for the simulations has an Intel ® Xeon ® Processor E5-1650 v4 which consists of 6 cores each capable of running 2 threads. This allows for the computation to be performed in 12 parallel threads. The computer also has 96 Gb of RAM allowing for extremely large computations.

This implementation was done using the multiprocessing package available in Python that drastically reduced computation times by a factor of 10.

Definition of source

The source in all cases mentioned in this thesis (unless otherwise mentioned) is defined by a Gaussian function (similar to the one defined in Section 2.2) multiplied by a tone burst signal. The Gaussian function in 2 dimensions (axial and circumferential) describing the spatial form of the source is given by:

G s (z -z s , a -a s ) = exp[-2/9 (z -z s ) 2 /(δz 2 )] exp[-2/9 (a -a s ) 2 /(δa 2 )],
where, (z s , a s ) are the axial and circumferential positions of the source, δz and δa are the discretisation steps sizes as shown in Figure 2.1. The coefficient 2/9 is chosen arbitrarily to ensure that the shape extends over atleast 8 points in the grid to satisfy the Nyquist-Shannon criterion. The tone burst signal of the source is given by:

T s (t) = ( √ 2π/n c ) sin(2 π f c t) exp[-2/n c (π t) 2 ],
where n c is the number of cycles and f c is the central frequency of the tone burst.

The source is hence given as,

S(z, a, t) = G s (z, a) T s (t)

Determining the length of the calculation

This section deals with identifying the length of the computation which would ensure that aliasing is avoided due to periodicity of the spatial domain, arising due to intrinsic periodicity of the Fourier transform.

The cylindrical domain is unwrapped as a rectangle. The upper and lower surfaces represent adjacent circumferential positions of the cylinder representing -π and π. The left and right edges represent the -z max and z max axial positions.

Consider an infinitely long tubular structure with a discrete source as defined in Section 1.4.3 and no damping. In the real case, guided waves generated by this source would radiate over the extreme axial edges of the domain and continue propagate outward to infinity. The waves that propagate over the -π edge would appear over the π edge.

In case of an axially finite medium, the waves propagating to the edges would reflect and continue to echo within the region.

In the case of the computation performed for a finite length of the pipe however, In the first case as seen in Figure 1.3, the source is located close to a circumferential edge at (z s = 0 m, a s = 0.8 π). Here, the periodicity of the computation plays to our advantage as the real domain is circumferentially periodic. Hence, waves propagating over the lower circumferential edge (a = -π) returns to the domain through the opposite edge (a = π). To better understand this phenomenon, a simulation is performed with a similar setup where the source is located at the center of the domain (z = 0, θ = 0). This behaviour determines the minimum axial length of the region of computation.

In the most common cases of Non Destructive Testing of tubular medium, the structure is illuminated by a principle axi-symmetric mode. Due to the nature of guided waves, there may also exist other modes generated which provide more information regarding the illuminated medium. These modes usually have different velocities.

The determination of the axial length for two common cases are described. The first case describes a situation where it is sufficient to illuminate a part of or the entire structure to visualize and understand the mechanism of propagation of the modes within it. In the second case, the interaction of the incident field with a defect and the subsequent measurement of the diffracted field by receivers may be envisaged.

For both cases, the source is considered to be at the center of the domain (z = 0)

and the incident field is assumed to contain two modes where the velocity of propagation of the faster mode is v f and that of the slower mode is v s . The region of interest has a length z max . The minimum axial length of the domain to avoid aliasing is z comp . For the second case, the defect is located at z def and the emitter transducers are also considered to be the receivers. These cases have been visualised in Figure 1.6.

Computation of axial length to illuminate a domain

To ensure that the domain is illuminated by the slower mode without aliasing of the faster mode, the minimum duration t min (time taken by the slower mode to span the region of interest) is given by:

t min = z int /v s .
Hence, the minimum length of the domain is:

z comp = t min v f .

Computation of axial length to illuminate and measure a defect

For the case with a defect, it is also necessary to consider the velocity of the slowest useful mode that may be diffracted by the defect v s_def . The time taken by the slowest mode of the incident field to reach the defect is,

t min_inc = z def /v s ,
and the time taken by slowest mode of the diffracted field to propagate back to the transducers would be:

t min_dif = z def /v s_def .
Finally, the minimum required axial length of the computation domain is given as,

z comp = (t min_inc + t min_dif ) v f .
However, with minute defects, it may be hard to estimate the modes that would be diffracted. To simplify this, a conservative estimate could be taken where

v s_def = 3 × v s ,
giving the axial length of the domain as:

z comp = 4/3 t min_inc v f ,
ensuring minimal to no aliasing in the medium.

Summary of model creation

This chapter introduces a new semi-analytical mathematical model used to simulate ultrasonic guided wave propagation in tubular structures. The first section introduces the model which is based on solving the problem in the Fourier domain for the spatial variables and the Laplace domain for the temporal variable. This spectral method offers a quicker alternative to FEM simulations mainly due to the simplicity of the domain as well as the nature of the computation to be performed (wave propagation). It allows for the problem to be solved as a function of independent spatial and temporal variables, which in turn allows for the computation to be performed in a distributed fashion. It also sidesteps the issue of finely discretizing the real domain, leading to shorter computation times. The section concludes by giving the partial wave solution of the wave propagation equations expressed as a sum of Bessel's functions of the first and second kind.

The discretisation domain is first introduced. The axial wavenumber domain is

given by Fourier transform variable, the circumferential wave number domain by a Fourier series expansion due to its periodicity and the temporal variable by the Laplace transform variable. The problem is solved analytically with respect to the radial position.

The next section introduced the idea of using a multiprocessing functionality to take advantage of the independence of the variables of the problem. The definition the source is then discussed. This definition is used during the entirety of this thesis. The source is defined as a product of Gaussian shapes in the axial and circumferential directions and an n cycle toneburst signal.

Finally, the factors determining the dimensions of the computation are discussed.

First, a limitation of the method related to the intrinsic periodicity of the Fourier and Laplace domains is considered. This is then followed by methods of calculating the dimensions of the domain for two simplified cases.

The model described in this chapter was used to create a fast numerical simulation tool on Python. This tool is used for all further computations. The validation of this tool is discussed in a further chapter of this thesis along with the modelling, simulation, identification and localization of numerical and experimental defects.

Chapter 2

Diffraction of a Field by a Defect

This chapter discusses the modelling of a numerical defect. The cylinder is modelled as two layers bonded together and the defect is considered as a 2d delamination, represented as a localised discontinuity of stresses and displacements at the interface of the two layers.

This defect is at radial position r I and is represented by surface D as shown in Figure 2.1. It can have any shape and has been shown here as a circle. Figure 2.1a shows the defect at the interface of the two layers. For ease of visualisation, the unwrapped surface containing the defect within the discrete grid of points has been shown in Figure 2.1b.

The axial dimension has a discretisation step of δz and that of circumferential direction is δa where δa = r I δθ. This representation of the defect in the grid is used as a reference for further computations mentioned in this thesis.

The field diffracted by the defect is approximated by a field generated by a secondary source located at the defect position. This source is modelled as a displacement jump at the interface r = r I which in turn is modelled as a sum of Gaussian functions and their derivatives. This secondary source depends on the radial stresses generated by the incident field on the defect surface D and is determined by the minimization of the total normal radial stresses on D.

Equations to solve

Consider a delamination defect at the interface of two layers at radius r I and centered at axial and circumferential positions z and θ respectively (refer Figure 2.2).

An incident field at the interfacial layer is given as σ

[inc] r (z, θ, r, s). To ensure that the axial and circumferential variables share the same units (meters instead of radians)

the angular variable θ is written in terms of a circumferential variable a where a = r I θ and r I represents the radius of the defect interface. It follows that the discretisation step is δa = r I δθ.

It must be noted that the method described below can also be used to simulate diffraction by a defect in plates. To equate the derivation for an arbitrary co-ordinate system, it would suffice to replace the variables (z, a, r) with (x, y, z) and the variables of the transformed domain (k, n, κ) with (k x , k y , k z ). However, to remain within the pedagogy of the current work, the co-ordinates of the cylindrical system are used.

The secondary source that corresponds to the defect response is described as a displacement jump at the interface and is written as:

∆u(z, s) = u(z, r + I , s) -u(z, r - I , s), (2.1)
where r + I and r - I refer to the position just above and below the interface and z denotes the spatial variable which could represent (z) for the 2d, axisymmetric cases and (z, a)

for the 3d, non-axisymmetric cases. The variable z has been separated from r as the computations are performed numerically with respect to z and analytically with respect to r. In the absence of a defect, the value of ∆u would be cancelled, thereby maintaining continuity between the two layers. This method of modelling a delamination defect is also valid when the two layers have different material properties.

The radial stress field diffracted by the defect is given by:

σ [dif] r (z, r, s) = G(z -ξ, r, s) ∆u(ξ, s) dξ, (2.2)
that represents a convolution in space, and the displacement field diffracted by the defect is given by:

u [dif] (z, r, s) = H(z -ξ, r, s) ∆u(ξ, s) dξ, (2.3)
where G(z -ξ, r, s) and H(z -ξ, r, s) are the Green tensors of radial stresses and displacement respectively. These Green tensors are the responses of equations 1.11 and 1.10 to a unit vertical displacement jump of a Dirac shape at the interface r = r I , located at an arbitrary position ξ. Mathematically, this is written as:

∆H(z, s) = δ(z) I 3 ,
where I 3 denotes the 3d identity matrix. The vectors Ũ(r), Σr (r) and F(r) are replaced by the matrices G(r), H(r) and δ(z) I respectively.

In its simplest form, the equation to be solved is :

∀z ∈ D , σ r (z, s) = σ [inc] r (z, s) + σ [dif] r (z, s) = 0, (2.4)
which says that the sum of the incident and diffracted field at the defect D is zero. This implies the absence of continuity of displacement at that location. Here, σ

[inc] r (z, s) is computed using the mathematical model introduced in Section 1.2 and σ

[dif]
r (z, s) is given by Equation 2.2.

The key element in this equation is ∆u(ξ, s) which corresponds to the sum of mul-tiple secondary displacement jump sources given as:

∆u(ξ, s) = i c i (s) b i (ξ) (2.5)
where, c i (s) are the amplitude coefficients of each of the secondary sources with shapes b i (ξ). This has further been detailed in Section 2.2.

Numerically however, the function to be minimized by the method of least squares is given as:

σ r (z, s) • σ r (z, s) φ(z) dz, (2.6)
where φ is the function that describes the shape of the defect. This shape function φ is typically given by:

φ(z) = 1 ; z ∈ D, φ(z) = 0 ; z / ∈ D.

Inner product associated with the defect shape

The minimization of Equation (2.6) requires us to have the values of the radial stress vectors at any point z. However, the numerical model only has the values at the grid of points of the discrete spatial domain. This can be resolved by the Shannon's interpolation formula giving the field at any position z as:

v(z) = n sinc(z -z n ) v n , (2.7)
where:

v n = v(z n ) , sinc(z) = sin π z δz δz π z (2d) , sinc[(z, a)] = sin π z δz sin π a δa δz δa π 2 z a (3d) ,
i.e. the field at a location anywhere on the grid (seen in Figure 2.1b) is given in terms of the sum of the field at all grid locations z n . The variable a = r I θ.

Next, an inner product of two discretised fields can be written as:

v, w = v(z) • w(z) φ(z) dz = n,m sinc(z -z n ) sinc(z -z m ) φ(z) dz an,m (v n • w m ) = V • (A W ) , (2.8)
where the real-valued symmetric positive square matrix

A = sinc(z -z n ) sinc(z -z m ) φ(z) dz n,m
, depends only on the shape of the defect and the grid of calculation. This signifies that the value of the field at a given location depends on the field at all discrete points on the grid albeit with magnitudes depending on its distance from the point. It is also important to notice that the interpolation need not be performed for the fields v and w.

The size of the defect is usually small compared to the size of the grid. As a result, the dimension of the subspace, where the inner product •, • has a significant effect on the response of the defect, is small compared to the dimensions of the grid N (number of grid points). Consequently, the number q of significantly positive eigenvalues of the square matrix A (i.e. its rank) is small compared to the order N of A.

The diagonalization of A gives q orthogonal unit eigenvectors b (with respect to the standard Euclidean inner product), with real-valued components, such that A b = λ b (λ > 0). If we consider the vectors q = √ λ b stack together to form the q-by-N matrix Q, we obtain:

A = Q T Q , i.e. v, w = (Q V) • (Q W ) . (2.9)
Finally, the inner product •, • is completely characterized by the Q matrix which is substantially less expensive to store than the A matrix (q×N against N 2 without using sparse matrix storage).

Thus, the use of the Q matrix ensures that the inner product is positive, while numerical errors in the computation, due to the avoidance of smaller inner product values of the A matrix, can generate tiny negative eigenvalues.

Creation of orthonormal basis: general case

As understood in Equation (2.4), in the ideal case, the function or functions that approximate the displacement jump to completely cancel the incident stress field within the contour of the defect is known. In reality however, it is almost impossible to exactly realize these functions. Hence, the functions are approximated by the Gaussian function and its derivatives. These functions that approximate the displacement jump which in turn represent the defect, are henceforth referred to as the "functions that approximate the defect" to maintain conciseness. It is also important to understand that these functions may not be linearly independent. The method of decomposition using wavelets has been discussed by multiple authors in different contexts [START_REF] Ososkov | Gaussian wavelet features and their applications for analysis of discretized signals[END_REF][START_REF] Rosyidi | Signal reconstruction of surface waves on SASW measurement using Gaussian Derivative wavelet transform[END_REF][START_REF] Wang | Modeling and Simulation of Nonstationary Processes Utilizing Wavelet and Hilbert Transforms[END_REF].

Due to the approximate nature of the defect modelling, these functions will not completely cancel the incident field i.e. there is no exact solution for Equation (2.4). The problem is thus converted to one of minimizing Equation (2.6). This requires for the functions that approximate the defect to be linearly independent. These linearly independent functions can be achieved using the Gram-Schmidt orthonormalisation process to form an orthonormal basis and is described below.

Consider a free family of m functions (θ α ) 1 α m (hereby referred to as free functions). These free functions may be linearly dependent on each other. The Gram-Schmidt process allows us to create a family of orthonormal functions

(ψ α ) 1 α m from (θ α ) 1 α m
iteratively and is built as follows:

         ψ 1 = d 1 θ 1 ; d 1 = θ 1 , θ 1 -1/2 and ∀α > 1 , ψ α = d α θ α - α-1 β=1 θ α , ψ β ψ β ; d α = θ α , θ α - α-1 β=1 θ α , ψ β 2 -1/2 .
(2.10)

The first vector of the orthonormal basis is built starting with a reference free function θ 1 and normalising it by taking its product with d 1 . The addition of a subsequent basis function ψ α is done by subtracting the projection α-1 β=1 θ α , ψ β ψ β of the newly added free function θ α on the pre-existing basis from the free function itself [refer Equation (2.10)]. The normalising factor d α is obtained by taking the inverse of the norm of the rejection of the new free function θ α on the pre-existing orthonormal basis. The rejection of a vector v on vector w is given by: Special care must be taken at this step as it is possible that the norm of the rejection is very small. This would make the new function θ α quasi-linearly dependent (Figure 2.4b) on the pre-existing orthonormal basis (ψ 1 , ψ 2 , ..., ψ α-1 ) and would lead to the newly added basis function, ψ α , to have an amplitude that is extremely large. In the real domain, this manifests as a divergent solution where the amplitude of the diffracted field is much larger than the incident field. In an ideal case, the newly added basis function would be linearly independent of the pre-existing basis as seen in Figure 2.4a.

rej v ( w) = v -proj v ( w) (2.11)
(a) Free function that is linearly independent of the pre-existing orthonormal basis due to a large rejection.

(b) Free function that is quasi-linearly dependent on the pre-existing orthonormal basis due to a small rejection. The extreme case where the rejection is 0 would mean that the new basis function is linearly dependent on the pre-existing orthonormal basis. This would signify the addition of a basis function with an infinite amplitude leading to an infinite amplitude of the field diffracted by the defect.

Any function v can be projected onto the linear span of the set of orthonormal functions (ψ α ) 1 α m to obtain:

P(v) = m α=1 v, ψ α ψ α , (2.12)
where the inner product v, ψ α represents the component c α (hereby referred to as coefficient) of v on the ψ α basis vector, which then minimizes

v -P(v), v -P(v) = |v(x) -P(v)(x)| 2 φ(x) dx, (2.13)
where P(v) represents the sum of the projections of the function v onto the vectors of the orthonormal basis and φ(x) represents the shape function described earlier in this section.

Gaussian secondary source and its derivatives

The numerical defect that this work deals with is equivalent to a delamination between two layers of identical or different material properties. If the defect D is not too large compared to the dimensions of the grid (seen as the circle in Figure 2.1b) the secondary source that generates the diffracted field can be approximated by a Gaussian function, positioned at the center of the defect, and its derivatives. In 1d, these higher order defects are represented as shown in Figure 2.5. It is to be noted that increasing the order of the derivative, increases the spatial resolution of the shape that describes the secondary source. This property can be used to simulate defects of smaller dimensions.

An important limitation however, is that increasing the order, moves the peaks of the shape function outside the contour of the defect. The sum of the negative values of the projections approximates the field diffracted by the defect.

A visual example has been shown in Figure 2.7 where the defect can be exactly approximated by a base with two basis vectors. This implies that the ideal subspace is 2d. Hence, the vector that describes the field diffracted by the defect, lies on the plane containing the two basis vectors (also the Defect Subspace). The incident field vector, shown as σ inc in the figure, is projected onto the defect subspace. The sum of the projections onto the different orthonormal basis vectors gives the vector that describes negative diffracted field -σ dif (shown as the light blue vector). The purple vector thus corresponds to the field diffracted by the defect.

In the ideal case, a perfect basis can be created such that the sum of the incident and diffracted field within the defect is zero i.e. the total field is orthogonal to the orthonormal basis. Here, the vector σ total (shown in dark green) describing the total field is linearly independent of the basis i.e. perpendicular to the 2d plane containing the defect response.

In the real case, however, it is not possible to obtain the exact vectors that create the ideal basis. This leads to an incomplete basis (shown as computed case) where the sum of the incident and diffracted fields within the defect is non zero. The total field vector σ total has an unresolved component in the defect subspace which means that the total incident field is not completely cancelled within the defect. This diffracted field acts as a secondary source with the input as a displacement jump at the position of the defect given by:

∆u(z, s) = n i=0 c i (s) ∂ i z g(z -z d ) (2.14)
where for the 2d case:

g(z) = exp[-2 z 2 /(9 δz 2 )];
∂ i z is the i th -order derivative with respect to z; c i (s) is a vector which depends on the incident field.

and for the 3d case:

g[(z, a)] = exp[-2 z 2 /(9 δz 2 )] exp[-2 a 2 /(9 δa 2 )]
∂ i z is the i th -order derivatives with respect to z and a:

∂ 1 z = ∂ z ∂ a , ∂ 2 z = ∂ 2 z ∂ z ∂ a ∂ 2 a , . . . ∂ i z ∂ i z is a vector with (i+1) components; c i (s) is a 3-by-(i+1
) matrix which depends on the incident field.

The coefficient 2/9 is arbitrarily chosen to create a Gaussian shape that spans 8 grid points to satisfy the Nyquist-Shannon criterion.

Equations (2.2) and (2.14) yield the diffracted field (vertical stress) on the interface:

σ [dif] r (z, s) = n i=0 G(z -ξ, r I , s) ∂ i z g(ξ -z d ) dξ c i (s) = n i=0 ∂ i z [G g (z -z d , s)]c i (s) = n i=0 ∂ i z G g z (z -z d , s) c z i (s) + ∂ i z G g a (z -z d , s) c a i (s) + ∂ i z G g r (z -z d , s) c r i (s) , (2.15)
where

G g z (z -z d , s), G g a (z -z d , s) and G g r (z -z d , s
) are the responses of the multilayered structure to a unit displacement step at the interface, of shape g(z -z d ), in the axial, circumferential and radial directions, respectively. These three responses are directly computed in the (k, r, s)-domain before performing a discrete Fourier transform to return to the real (z, r, s)-domain.

Consequently, we want to find the three components of c i (s), namely c z i (s), c a i (s)

and c r i (s), that minimize

σ [inc] r (•, s) + σ [dif] r (•, s), σ [inc] r (•, s) + σ [dif] r (•, s) i.e., σ [dif]
r (•, s) is the negative of the sum of the projections of the normal stress σ

[inc] r (•, s)
of the incident field, on the basis

∂ i z G g z (• -z d , s), ∂ i z G g a (• -z d , s), ∂ i z G g r (• -z d , s) 0 i n .

Definition of secondary source: Summary

The process of defining a defect given by a secondary source can be summarised as follows.

First, the equation for the total stress to be minimised σ r (z, s) has been described simply by equation 2.4. The secondary source defined as ∆u(ξ, s) is given as a sum of multiple secondary sources b i (ξ) with their corresponding coefficients c i (s). The aim is to find the coefficients for this set of secondary source shapes which minimise the total stress field within the defect D. The method of least square minimisation is then employed to minimise the total stress within the defect and is described by equation 2.6.

Next, the Shannon interpolation formula is introduced (Equation 2.7). This formula allows to perform a computation of a continuous space given a set of discrete grid points hence converting a continuous computation to a discrete one. The development leads to the creation of the Q matrix that only depends on the shape of defect and the points of computation i.e. the Q matrix contains the defect properties.

The process of creating an orthonormal basis is then introduced and is shown by equation 2.10. This process allows for the creation of a family of orthonormal functions (ψ α ) 1 α m from an arbitrary set of free functions (θ α ) 1 α m . In this thesis, the set of free functions are selected as the Gaussian function and its derivatives. Finally the coefficients c i (s) are represented by σ inc , ψ i (given by equation 2.12) implying that the coefficients are found by projecting the incident radial stress field on each of the functions of the orthonormal basis.

Defect studies based on defect size

In this section, studies are performed to understand the response of the defect based on its size as well as the order of the derivative of the Gaussian function that approximates it. It is to be noted that the Gaussian function (considered to be order 0) exists in the axial and circumferential directions. Throughout the rest of this work, any mention of an order includes all the orders that exists below it. For example, order 1 would consist of the basis (0, 0), (0, 1), (1, 0) for each of the 3 Green tensors, i.e. basis with 9 elements.

A basis consists of a set of vectors such that every element in the basis can be described The spectrum of the input signal has a central frequency of 250 kHz. However, due to the excitability curve of this setup, the response of the medium is shifted towards the lower frequencies. This concept has been detailed well in a number of articles [START_REF] Wilcox | Dispersion and excitability of guided acoustic waves in isotropic beams with arbitrary cross section[END_REF][START_REF] Tran | Excitability of ultrasonic lamb waves in a cortical bone plate: A simulation study[END_REF][START_REF] Samaitis | Influence of the spatial dimensions of ultrasonic transducers on the frequency spectrum of guided waves[END_REF][START_REF] Treyssède | Numerical and analytical calculation of modal excitability for elastic wave generation in lossy waveguides[END_REF].

The study is performed for radii of 2 mm, 4 mm and 6 mm and orders ranging between 0 (Base Gaussian) and 5. Order 0 consists of an orthonormal basis of size 3 (one for each of the 3 directions (r, θ, z)) and order 6 consists of 63 element basis. In each of these cases, the total axial stress field (incident + diffracted) at the center of the defect The total stress field at the center of the defect has been plotted as a function of frequency for each of the different derivatives and is shown in Figure 2.9. This is done to study the effectiveness of the different order derivatives in cancelling out the incident stress field within the defect. Also overlayed on this plot is the incident stress field obtained at the same point in a defectless medium as a reference.

Figure 2.9: The spectrum of the total field measured at the center of the 6 mm radius defect for derivatives of different orders. The dotted line denoting 10 % of the maximum amplitude below which the defect approximation may be considered effective.

The series of images in Figure 2.10 gives a snapshot of the total field for four different orders of derivatives for the same time step of 45 µs. From Image 2.9, it would seem that the 4th order derivative may be sufficient to model the defect by minimising the total field. However, from Figure 2.10c, it is clear that the incident field is not completely cancelled near the edges of the defect. The total field within the contour for the 5th order approximation is sufficiently close to zero and could be considered as a suitable condition for a defect of this dimension.

Case 2: radius of defect 4 mm

As in the first case, the images in Figure 2.12: The spectrum of the total field measured at the center of the 4 mm radius defect for derivatives of different orders. The dotted line denoting 10 % of the maximum amplitude below which the defect approximation may be considered effective.

From Figure 2.12 it can be observed that like the earlier case, using the 5th order derivative is effective in minimising the total field within the contour. However, due to the smaller radius of the defect, a smaller order of derivative i.e. 4th order would also be sufficient. This can be verified in the snapshots in Figure 2.11. The 3rd order derivative is also almost effective but fails to cancel the stress near the contours. Unlike the two earlier cases discussed in this study, the response of the defect from order 4 and 5 diverges. This is explained by the fact that the peaks of the higher order functions are further away from the contour of the defect. To cancel the incident stress field in the defect, they would need to have a higher amplitude, leading to a divergent solution. The images in Figure 2. [START_REF] Davies | Synthetic focusing for high resolution guided wave pipe inspection: Further results and robustness studies[END_REF] shows the snapshot of the total field for 4 different orders at the same time step. The 1st order (Figure 2.16b) is effective in cancelling the incident field within the defect. In case of the 3rd order derivative seen in Figure 2.16c, the field diffracted by the defect has an amplitude comparable to the incident field. This signifies that the solution is on the verge of diverging. Additionally, increasing the order to 4, diverges the solution for reasons mentioned in section 2.1.2. This is seen with the amplitude of the diffracted field being much larger than the incident field.

Observations and conclusions

The small sample of defect studies performed, give a large variety of results, observations and conclusions. The most direct observation is that the number of derivatives used to approximate the defect depends on the size of the defect with respect to the grid. A larger defect requires higher order derivatives. This can be explained by the shapes of the functions that are used for the approximation. A higher order derivative not only covers a larger area but also contains a greater number of peaks. This allows the model to effectively cancel out incident stress fields with wavelengths of comparable length to that of the defect size.

There is, however, a mathematical upper limit to the maximum order that can be used to simulate the defect based on its size. In the case of the 2 mm defect, the solution diverges at the 4th order whereas in the case of the 6 mm defect, the solution stays stable until the 5th order. The diverging solution is due to the method by which the orthonormal basis is built. As discussed in section 2.1.2, the normalisation of a newly added basis requires it to be multiplied by the inverse of the norm of the rejection. This rejection value depends on the shape of the defect as well as the discretisation step length.

A smaller defect size leads to the rejection having a smaller value at lower orders and vice versa.

There is also a numerical limit to the maximum order of the derivative that can be used. This limit depends on the discretisation step size as well as the standard deviation of the Gaussian basis function and arises from the Nyquist-Shannon criterion. If the number of grid points are not sufficient to accurately describe the curved shape of the function, aliasing may occur. In the case of the examples above, this limit is reached around the 6th or 7th order. This can be resolved by placing multiple Gaussian functions within the defect, where each Gaussian function and its derivatives are added to the orthonormal basis.

From the above studies, certain conclusions can be drawn regarding the conditions suitable to accurately approximate the field diffracted by sufficiently large circular defects.

The peaks of the Gaussian functions and its derivatives used should lie either within or at the immediate vicinity of the defect contour and the absolute amplitude of the diffracted field should remain lesser than or equal to that of the incident field.

The basis can also be built with a lesser number of derivatives depending on the shape of the defect and how well it conforms to the shape of the initial Gaussian basis function. For example, consider an elongated defect (refer Figure 2.17). In such a case, the 0th order derivative would suffice for the axial direction. It would be necessary to take multiple derivatives in the circumferential direction so as to have the peaks of the secondary sources that cover the surface within the contour. However, a thorough study needs to be performed to understand the basis requirements for these sort of defect shapes which only partially conform to the shape of the Gaussian basis function function.

The lower limit of the model to approximate small defects depends on the discretisation step and the Gaussian basis functions used to approximate it. Extremely small defects would be approximated by secondary sources which are considerably further away from the defect itself. This could be resolved by decreasing the discretisation step size which comes at the price of a costlier computation. It is also possible that an alternate basis function is used instead of the Gaussian. A new narrower optimised function has been developed that would simulate smaller defects more accurately than the Gaussian.

However, this line of work requires further studies to better understand and implement it.

The contour of the defect also plays a part on the functions used to approximate it. In the examples discussed above, the contour of the defect is regular and conforms to the shape of the Gaussian function. However, there may be cases where the defect is irregularly shaped, as seen in Figure 2. [START_REF] Davies | The application of synthetic focusing for imaging cracklike defects in pipelines using guided waves[END_REF]. In such cases it could also be possible to approximate the defect with two or multiple Gaussian functions (seen in Figure 2.18b).

Each of the Gaussian functions and their derivatives would be added as a separate basis function to the overall orthonormal basis.

Experimental setup

Details of experiment

The experimental setup consists of a 6 m long commercially obtained 2017A Aluminum pipe (procured from Almet) with inner and outer diameters being 56 mm and 60 mm respectively. The exact physical characterisation of the material is discussed in a later section.

The transducers used to excite the guided modes come from Guided Ultrasonics Ltd. The transducers from the emitter belt are then connected to the multichannel signal acquisition system (by LeCoeur). The impedance of the transducers and the acquisition system are not matched leading to poor performance in the pulse echo mode with a very low signal to noise ratio. As a workaround to this problem, a second belt of transducers is used in tandem with the first to work in the pitch catch mode. The receivers are connected to the acquisition system through an impedance adapter to improve the signal to noise ratio. An advantage of this setup is that the signals can be measured at various distances from the emitter ring which can then be compared with simulations. Figure 3.3 shows the schematic of this setup. The acquisition system is run through a Matlab script.

Normalisation of transducers

The normalisation of the test setup is an important step to ensure the generation of pure guided modes as well as to have the exact amplitudes at the receivers. As is the case of maintaining an equi-angular spacing between the transducers, an inconsistency in the amplitudes in the emitter belt would lead to the generation of unwanted modes.

There is an unavoidable variation of amplitudes applied by each transducer onto the structure of the structure. This may be attributed to slight differences in manufacturing, age of transducers or the clamping mechanism. Application of the shear coupling gel may also lead to a disparity between generated amplitudes. It is thus essential to normalise all transducers in the emitter ring to ensure identical amplitudes are generated. It is equally important to normalise the transducers in the receiving ring so as to have measurements of comparable amplitudes. These measurements are to be used at a later stage for the localisation of the defects with Topological Imaging.

The normalisation of all transducers is executed by a calibration step performed before each round of experimental readings. This calibration allows for consistent and repeatable measurements and is valid when there are an equal number n of emitters and receivers and when spatial circumferential periodicity is maintained.

After setting up the experiment with an emitter and a receiver ring, each of the transducers in the emitter ring is excited in turn and their individual response recorded at each of the receiver locations. The maximum amplitudes received in each of the cases is then stored in the following square matrix:

M =             m 00 m 01 . . . . . . m 0(n-1) m 10 . . . . . . m ij . . . . . . m (n-1)0 m (n-1)(n-1)            
where, i and j refer to the emitter and receiver positions respectively.

Each of the emitters and receivers have gains of e i and r j associated with them.

These coefficients are applied to the respective emitter and receiver to ensure that they are normalised. There is also a transmission coefficient t j-i between each emitter receiver pair which describes the ideal coefficients that exists between them. A visual representation of these relations is shown in Figure 3.4. In the ideal case it is clear that the transmission coefficient between the emitter e 0 and receiver r 0 should be equal to that between e 1 and r 1 and so on. In general, the transmission coefficient between emitter e i and receiver r j is given by t k where k is j -i. This is represented as :

m ij = e i t j-i[n] r j .
Taking the natural logarithm of this equation we get:

log(m ij ) = log(e i ) + log(t j-i[n] ) + log(r j ),
which is represented simply as:

µ ij = i + τ j-i[n] + ρ j , (3.1)
where τ j-i can be represented by a circular matrix as shown below:

Emitter i → Receiver j ↓                τ 0 τ 1 • • • τ j-1 • • • τ n-1 τ n-1 τ 0 • • • τ j-2 • • • τ n-2 . . . . . . . . . . . . . . . . . . τ n-i+1 τ n-i+2 • • • τ j-i • • • τ n-i . . . . . . . . . . . . . . . . . . τ 1 τ 2 • • • τ j • • • τ 0               
.

This equation has an infinite number of solutions. However, we make an assumption that the product of the amplitudes of all transducers in either the emitter or receiver belts is unity. This leads to a unique solution and in logarithmic terms, corresponds respectively to:

µ ij = n ij τ j-i[n] = n k=1 nτ k = n 2 μ,
where μ is the average value of µ ij Due to the circularity of the matrix, the vector τ k is repeated n times. Hence, the transmission coefficients can be written as,

τ k = 1 n j µ j+k[n],j (3.3)
i.e. mean of µ ij over the (i, j) pairs such that i -j ≡ k mod (n).

1

n j µ ij -μ ; ∀i . (3.5)
The equations (3.3), (3.4) and (3.5) are inverted from the log scale to the real scale by taking their exponential to get the actual coefficients of transmissions and gains to be applied to the emitters and receivers.

This method ensures that the emitter and receiver belts are normalised to have amplitudes close to unity with respect to each other. Before the experiments are performed, a calibration step is carried out where the response of each emitter is recorded at each receiver position. Using the aforementioned method, all transducers in the setup are normalised before proceeding with the actual experiments.

Attenuative influence of multiple transducers on incident field

The effect of transducer contact on the surface of the cylinder has an attenuative effect on the generated signal and was been studied here. This effect was studied at the position of the emitter ring. The material attenuation of the signal is seen regardless of the number of added transducers. It is also clear that the addition of multiple transducers reduces the amplitude of the circumferentially propagating wave. Attenuation due to material and transducers is not taken into account in the simulations and the significance of this is better understood in the chapter of Topological Imaging (Chapter 4). There is, however, a possiblity of including damping in the properties through the diamond product as defined in Section 1.2 and has been described in a paper by Mora et al. [START_REF] Mora | Transient 3D elastodynamic field in an embedded multilayered anisotropic plate[END_REF].

3.2 Dimensioning of numerical model

Material characterisation

Before running the simulations, it is essential to characterise the material to obtain the properties for the numerical simulations. The Hydro-static weighing technique, being one of the most accurate and convenient, was employed to determine the density of the test specimen. There are, however, various methods to obtain the bulk velocities of a material. Some of these have been detailed in [START_REF] Selfridge | Approximate Material Properties in Isotropic Materials[END_REF][START_REF] Hsu | Simultaneous ultrasonic velocity and sample thickness measurement and application in composites[END_REF][START_REF] Franco | Measurement of Elastic Properties of Materials By the Ultrasonic Through-Transmission Technique Medición De Las Propiedades Elásticas De Materiales Por El Método De Transmisión Ultrasónica[END_REF][START_REF] Kiefer | Simultaneous Ultrasonic Measurement of Thickness and Speed of Sound in Elastic Plates Using Coded Excitation Signals[END_REF][START_REF] Grünsteidl | Determination of thickness and bulk sound velocities of isotropic plates using zero-group-velocity Lamb waves[END_REF]. Most of these studies deal with flat plates allowing for through transmission measurements. In this study a method involving the simplest methods are employed to obtain the material properties.

A small portion of the test specimen is cut out from the pipe. According to the material sheet, the pipe has an inner diameter of 56 mm and an outer diameter of 60 mm giving it a thickness of 2 mm. On careful measurement though, the thickness is found to vary between 1. Taking the Fourier transform of this signal, we get the frequency spectrum as in Figure 3.10. This spectrum consists of a number of equidistant peaks (frequency combs)

which contain the information of the time period of oscillation. To obtain the period from these regularly spaced frequency peaks (∆f ), another Fourier transform is performed on the absolute value of this spectrum (similar to the method in [START_REF] Lafargue | Localized lasing modes of triangular organic microlasers[END_REF][START_REF] Maksymov | Synthesis of discrete phase-coherent optical spectra from nonlinear ultrasound[END_REF]). The result is as shown in Figure 3.11, where the first peak gives an accurate estimate of the time period between the multiple reflections. Similar measurements were taken multiple times and averaged. Taking thickness to be 1.95 mm, we get a longitudinal bulk velocity of 6254 m/s.

A similar study was performed to obtain the transverse bulk velocity to be 3052 m/s. There are a couple of sources of error for these bulk velocity values. Placing the flat surface of the transducer exactly normal to the surface of the cylinder is a challenge. In practice, a conceivable variation of about 5°from the normal would give a difference of about ± 0.3% error. To minimize the effect, five separate measurements with 10 averages for each case were taken. Another source of error is the measurement of the wall thickness.

As mentioned earlier, the thickness varies between 1.92 mm to 2.02 mm. This gives a possible error of about ± 2.5% in the velocity value. The test was performed with a sampling frequency of 40 MHz. The peak of the wave packets are measured within an accuracy of 0.025 µs. A total of 8 reflections are used for the estimation of time t with the final estimate of time being 0.62 µs. This in turn causes an error of about ± 0.67%. The total estimated error is in the order of ± 3.5%. It is necessary to take the final velocity values with a pinch of salt.

These calculated values of C L and C T seem to be on the lower end of the spectrum of the bulk velocities in literature, however an article by Turnbull et al. [START_REF] Ginzel | Determining approximate acoustic properties of materials[END_REF] show lower experimentally calculated bulk velocity values and explains the multitude of reasons this could be true.

These material properties are fixed for all future studies and listed out in Table 3

Modelling of source: transducer sizing

A numerical study was performed to understand the variation of the field as a function of the dimensions of the transducer. The transducer is modeled as a Gaussian distribution along the circumferential and axial directions (described in Section 1.4.3). The actual length and width of the transducer contact area are 13 mm and 3 mm, respectively. To test the effects of the dimensions of the Gaussian describing the transducer, 3 sizes were chosen in each of the two directions of the transducer. These sizes were 80 %, 100 %, 120 % of the actual axial length and circumferential width. The dimensions chosen are listed in Table 3 The lower limit of the width was fixed by approaching the limit of the calculation of the Bessel function. At extremely high circumferential wave numbers, this function leads to singularities by reaching the limits of the computer beyond which the number is considered as infinity. This limits the maximum circumferential wave number which in turn limits the minimum dimension of the circumferential width of the simulated transducer source. Simulating smaller transducer sizes would require a smaller discretisation step, thus increasing the number of points of calculation for an already costly computation. It is thus necessary to balance the accuracy and cost of the computation.

The incident displacement field for the L(0, 2) mode was simulated by taking all combinations of the lengths and widths. The Z-directional displacement component of these fields at the surface were measured at a distance of 50 cm from the emitter belt and at an angular position of θ = 0. Comparison of these fields show their susceptibility to various dimensions of the transducer. The simulations were run with a 6-cycle Gaussian window input signal with a central frequency of 100 kHz to match the experiment. At this frequency, the primary L(0, 2) mode used for this study is only slightly dispersive. This is a useful outcome as the simulation can be performed with a larger transducer width, reducing the cost of the computation. However, to remain as close as possible to the experimental case, it was decided that a width of 3.6 mm (120 %) to be used for all further simulations.

Effect of circumferential width on amplitude and wavefront

Generation and comparison of modes : experiment vs simulation

In this section, the simulation model is validated experimentally by comparing the behaviour of three modes, namely the L(0, 2), F (1, 3) and L(4, 2), to experimental results.

The modes are chosen so as to verify purely axi-symmetric as well as flexural behaviours of lower and higher order circumferential modes.

According to the naming convention of guided wave modes in cylinders as described by Silk and Bainton in 1979 [START_REF] Silk | The propagation in metal tubing of ultrasonic wave modes equivalent to Lamb waves[END_REF] and Nishino in 2001 [START_REF] Nishino | Modal Analysis of Hollow Cylindrical Guided Waves and Applications[END_REF], the first index refers to the number of complete wavelengths around the circumference of the cylinder and the second index determines the mode number. As a side note, trying to generate a pure L(0, 2) mode also leads to the generation of the F (8, 1) and L(0, 1) modes as seen in the earlier section. The F (8, 1) is created as there are 8 equally spaced emitters around the circumference and the fundamental mode of the circumferential wavenumber N = 8 has a cut-off frequency lower than the testing frequency of 100 kHz. Figure 3.15 shows the fundamental modes of each of the circumferential wavenumbers from 8 -15. These dispersion curves have been computed using the method described in Section 1.3.2. The modes with circumferential order 0 -13 have a cut-off frequency less than 100 kHz and hence using upto 13 transducers all in phase in an attempt to generate a pure axi-symmetric mode such as L(0, 2) would lead to the unavoidable generation of one of these modes corresponding to the number of transducers in the emitter belt. On the other hand, using 14 transducers would generate an extremely dispersive mode with a group velocity of almost zero, making the signal incoherent. Ideally, for this given system, using 15 or more equiangularly placed transducers activated in phase would generate a pure L(0, 2) mode. However, the physical width of the transducers hinders the positioning of more than 14 transducers around the circumference.

The transducers contact each other, leading to insufficient contact with the surface of the cylinder. This in-turn leads to the generation of other modes which would unnecessarily complicate the wave field. The L(0, 1) is generated due to the axi-symmetric nature of excitation as well as the mode having a small but significant z-directional displacement component at the surface. The L(0, 1) and F (8, 1) modes belong to the same family and are similar to the anti-symmetric(A 0 ) mode in plates. They also have identical variations of stresses and displacements over the thickness. The major differences are that the F (8, 1) has a slower group velocity and the amplitude varies with respect to the angular position, consisting of 8 complete wavelengths around the circumference. Figure 3.16 shows the spectrum of the calculated field in the frequency/axial wavenumber domain. The highlighted regions show the extent to which the three modes are activated in the simulation and this is overlayed on the analytical dispersion curves. This explains the two packets that follow the L(0, 2) mode.

Figure 3.16: Frequency-axial-wavenumber spectrum of the computed field overlaid on relevant dispersion curves showing a good agreement between the calculated and theoretical solutions.

The numerical tool also closely simulates the response of an undamaged cylindrical structure. This is verified by comparing the experimental and theoretical signals of various modes measured by 8 equiangularly spaced transducers positioned at 50 cm and 100 cm from the transducer belt. Figure 3.17 shows that the axi-symmetric mode L(0, 2) is exactly replicated by the experiment. Figure 3.18 shows that the tool also simulates the response of a flexural mode F (1, 3) pretty accurately and Figure 3.19 shows that the tool can not only simulate extremely dispersive modes (in this case the F (4, 2) mode) but also retains the exact group velocity of said dispersive mode. In case of the simulation of the F (4, 2) mode measured at 50cm, the response of the circumferentially propagating mode also exists and is seen by the repeating signals after the primary packet of the mode. This response is not seen in the experiment as there is a strong damping factor due to the material and interaction of the fields with the emitters at the source location as discussed in section 3.1.3, both of which are not taken into account in the simulation.

Capabilities and limitations of the numerical model

Before moving on to the defect detection using the method of Topological Imaging, a brief summary is discussed regarding the capabilities and limitations of the tool and experimental setup.

The model allows for a fast and accurate simulation of guided wave propagation in tubular structures. This arises from the fact that the computation is performed in the Fourier and Laplace domains, thus avoiding the complexity of solving the problem with higher order derivatives that exists in the real domain. Using the Laplace variable for the time domain instead of the Fourier variable also allows for truncated causal computations.

In addition, the independence of the variables allows for the computation to be performed using multiple processors which further reduces the total computation time. The tool is also capable of simulating multi-layered structures which could have different material properties. This may, for example, be used for simulation of guided wave propagation in insulated pipes. It can also simulate wave propagation from single or multiple discrete localized sources hence simulating individual transducers.

The method of computation also allows for the simulation of defects. This thesis discusses diffraction caused by delamination defects of different sizes at any position in the medium. Irregularly shaped defects are also approximated well by using multiple Gaussian functions and their derivatives. Further work is ongoing to simulate volumetric defects such as indentations and cracks.

The simulation also holds up well in relation to analytical models and experimental measurements. The frequency wave-number spectrum images show the spectrum of the simulated fields in line with the analytical dispersion curves. Comparing the simulation and experimental measurements of various modes shows an accurate reproduction of the phase, velocity and waveform. This is also true for higher order and dispersive modes.

Notwithstanding the accuracy and validity of the numerical model, there are a couple of numerical and experimental limitations. One of the numerical limitations is with regard to the circumferential wavenumber due to the modified Bessel's function. The simulated transducer is approximated by a Gaussian distribution that requires a minimum of 8 points to be well defined. This controls the circumferential discretisation step which in-turn decides the maximum circumferential wave number value to accurately run the simulation satisfying the Nyquist-Criterion. The Bessel's function is unstable at large values of N . However, as seen in the study regarding the transducer sizing, this circumferential width of the transducer has a comparatively small effect on the calculated field and the transducer can be approximated by a Gaussian with a larger width lowering the maximum circumferential wavenumber value. This limitation also affects the minimum dimension of the numerical defect. Higher order derivatives of the Gaussian functions have peaks that lie outside of smaller defect contours leading to an instability in the diffracted field.

Another limitation is with regard to the dimension of the computation region. Due to the periodic nature of the Fourier domain, care has to be taken to select the physical dimension and duration of the computation to ensure aliasing does not occur.

There were also some experimental limitations. One such limitation is the upper limit to the number of transducers which can physically be placed around the pipe without contacting each other. This in combination with a limit of the transducer central frequency used meant that it was impossible to only generate pure axi-symmetric modes (N = 0) without also generating the highly dispersive higher order modes (N = 8).

One other issue was the inherrent dissimilarity in the amplitudes generated by each source due to factors such as transducer strength, channel differences, couplant. This led to the generation of unwanted trailing signals and was resolved by the normalisation of the transducers.

On the whole, the versatility of the tool in simulating a vast extent of practically useful cases has been demonstrated.

Chapter 4 Topological Imaging

The numerous methods realised to simulate wave propagation go hand in hand with a multitude of post-processing techniques developed to image defects. There exists a number of techniques for a vast variety of applications within the fields of NDT, Seismology, Medicine etc. Most of these techniques are usually cross compatible, hence allowing them to be adapted to other fields.

One of the earliest and robust methods developed for imaging defects is the Synthetic Aperture Focusing Technique (SAFT). This method involves the capture of data from multiple positions and the application of synthetic focusing algorithms to arrive at the defect location. A review by Busse et al. [START_REF] Busse | Review and discussion of the development of synthetic aperture focusing technique for ultrasonic testing (SAFT-UT)[END_REF] details the method and it has been effectively used since then by a number of authors [START_REF] Burch | Ultrasonic Synthetic Aperture Focusing Using Planar Pulse-Echo Transducers[END_REF][START_REF] Doctor | SAFT -the evolution of a signal processing technology for ultrasonic testing[END_REF][START_REF]Ultrasonic characterisation of defects in steel using Multi-SAFT imaging and neural networks[END_REF][START_REF] Karaman | Synthetic Aperture Imaging for Small Scale Systems[END_REF][START_REF] Blouin | Improved resolution and signal-to-noise ratio in laser-ultrasonics by SAFT processing[END_REF][START_REF] Sicard | A SAFT algorithm for lamb wave imaging of isotropic plate-like structures[END_REF][START_REF] Wilcox | Omni-directional guided wave transducer arrays for the rapid inspection of large areas of plate structures[END_REF][START_REF] Pitkänen | SAFT -Is it a Tool for Improved Sizing in Ultrasonic Testing?[END_REF][START_REF] Davies | Synthetic focusing for high resolution guided wave pipe inspection: Further results and robustness studies[END_REF][START_REF] Gaul | Localization of Defects in Pipes Using Guided Waves and Synthetic Aperture Focussing Technique (SAFT)[END_REF]. A study by Kitze et al. [START_REF]SAFT and TOFD -Compared in Ultrasonic Defect Detection[END_REF] compares SAFT with Time of Flight Diffraction (TOFD) technique.

The TOFD method uses the reflection from the back wall and the direct path signal as a reference to determine the positions and dimensions of the defects. The SAFT technique has also been used for defect detection in tubular structures [START_REF] Davies | The application of synthetic focusing for imaging cracklike defects in pipelines using guided waves[END_REF]. A review by Davis et al. [START_REF] Davies | Review of synthetically focused guided wave imaging techniques with application to defect sizing[END_REF] summarises and compares some Synthetic Focusing Techniques such as SAFT, Common Source Method (CSM) (A technique very similar to SAFT which is especially effective in dealing with tubular structures) [START_REF] Fletcher | Detection of axial cracks using focused guided waves[END_REF][START_REF] Fletcher | Detection of axial cracks in pipes using focused guided waves[END_REF] and Total Focusing Method (TFM).

The TFM method works with the Full Matrix Capture (FMC) acquisition technique where all available time traces between transmitters and receivers are stored and processed for imaging. These techniques were made possible by the development of Phased Array systems, one of the earliest mentions of which, were by Mahaut et al. [START_REF] Mahaut | Development of phased array techniques to improve characterization of defect located in a component of complex geometry[END_REF]. An article by Holmes [START_REF] Holmes | Post-processing of the full matrix of ultrasonic transmit-receive array data for non-destructive evaluation[END_REF] summarises some of these FMC techniques concisely with some interest-ing results showing the effectiveness of the TFM. This technique has been successfully implemented by various authors [START_REF] Holmes | Advanced post-processing for scanned ultrasonic arrays: Application to defect detection and classification in nondestructive evaluation[END_REF][START_REF] Njiki | Total focusing method for non destructive evaluation: Toward realtime imaging systems[END_REF][START_REF] Stratoudaki | Adapting the full matrix capture and the total focusing method to laser ultrasonics for remote non destructive testing[END_REF][START_REF] Dumas | Adaptive ultrasonic imaging with the total focusing method for inspection of complex components immersed in water[END_REF][START_REF] Deutsch | High Resolution Phased Array Imaging using the Total Focusing Method[END_REF][START_REF] Lin | Total focusing method imaging of multidirectional CFRP laminate with model-based time delay correction[END_REF].

Some of the other imaging techniques include probabilistic methods [START_REF] Bagheri | Ultrasonic imaging algorithm for the health monitoring of pipes[END_REF][START_REF] Dehghan-Niri | A multi-helical ultrasonic imaging approach for the structural health monitoring of cylindrical structures[END_REF][START_REF] Wang | A sparse sensor network topologized for cylindrical wave-based identification of damage in pipeline structures[END_REF], combination of the TFM and Sign Coherence Factor [START_REF] Liu | Damage localization in aluminum plate with compact rectangular phased piezoelectric transducer array[END_REF] and Helical Ultrasound Tomography [START_REF] Leonard | Guided wave helical ultrasonic tomography of pipes[END_REF].

Topological Imaging (TI) has its roots in the field of shape optimisation. The idea of shape optimisation was first discussed by Eschenauer, Schumacher et al. [START_REF] Eschenauer | Bubble method for topology and shape optimization of structures[END_REF].

The method was then generalised by Sokolowski [START_REF] Sokolowski | On the Topological Derivative in Shape Optimization[END_REF] and extended to heat transfer and linear elasticity by others [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF][START_REF] Novotny | Topological sensitivity analysis[END_REF][START_REF] Samet | The Topological Asymptotic for the Helmholtz Equation[END_REF]. The idea of this method is to study the point wise sensitivity of a cost function in the presence of holes at those points. Some of the notable works with regard to imaging defects using this method are by Amstutz [START_REF] Amstutz | Topological sensitivity analysis in the context of ultrasonic non-destructive testing[END_REF], Dominguez [START_REF] Dominguez | Time domain topological gradient and time reversal analogy: An inverse method for ultrasonic target detection[END_REF][START_REF] Dominguez | Non-destructive imaging using the time domain topological energy method[END_REF], Guzina [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], Bellis [START_REF] Bellis | Crack identification by 3D time-domain elastic or acoustic topological sensitivity[END_REF] and Bonnet [8] going by the name of Topological Sensitivity, Topological Gradient and Topological Derivative interchangably.

The study by Dominguez et al. [START_REF] Dominguez | Time domain topological gradient and time reversal analogy: An inverse method for ultrasonic target detection[END_REF] and Bonnet [START_REF] Bonnet | Topological sensitivity for 3D elastodynamic and acoustic inverse scattering in the time domain[END_REF] discusses the Time Domain Topological Gradient (TDTG) method where the computation of the direct and adjoint problems for target detection based on the idea of minimising a cost function has been successfully

implemented. These methods are iterative and require multiple computations.

A crucial concept for the functioning of Topological Imaging is the time reversal technique. This technique for focusing on reflective targets was first mentioned by Fink et al. [START_REF] Fink | Self focusing in inhomogeneous media with time reversal acoustic mirrors[END_REF][START_REF] Fink | Time Reversal of Ultrasonic Fields-Part[END_REF] named the Time Reversal Mirror (TRM). In essence, this technique shows that that diverging diffracted signals from a scatterer recorded at multiple positions transmitted back into the medium after time reversal tend to converge at the position of the scatterer. An extreme case of the reversibility of the wave-fields was proven by Derode et al. [START_REF] Derode | Robust acoustic time reversal with high-order multiple scattering[END_REF] in an experiment with 2000 scatterers where the signal converges back to the source with the original wave form. In the same year Prada et al. [START_REF] Prada | Decomposition of the time reversal operator: Detection and selective focusing on two scatterers[END_REF] introduced the Decomposition of the Time Reversal Operator (French abbreviation DORT) to separate out signatures and selectively focus on two scatterers with the same TRM technique.

Leutenegger et al. [START_REF] Leutenegger | Detection of defects in cylindrical structures using a time reverse method and a finite-difference approach[END_REF] later successfully used this method to localise defects by obtaining the spatio-temporal maxima of the time reversed field. The principle of time reversal to solve inverse problems in multiple scattering media has also been studied by authors such as Blomgren [START_REF] Blomgren | Super-resolution in time-reversal acoustics[END_REF] and Devaney [START_REF] Devaney | Super-resolution processing of multi-static data using time reversal and MUSIC[END_REF] using the Multiple Signal Classification (MUSIC) imaging algorithm.

The Time Domain Topological Energy (TDTE) method was introduced by Dominguez, Sahuguet and Gibiat [START_REF] Dominguez | Non-destructive imaging using the time domain topological energy method[END_REF][START_REF] Sahuguet | Biological tissues imaging with time domain topological energy[END_REF][START_REF] Gibiat | Wave guide imaging through time domain topological energy[END_REF]. This was soon followed by Rodriguez et al. [START_REF] Rodriguez | Fast topological imaging[END_REF] who introduced frequency domain topological imaging [Fast Topological Imaging Method (FTIM)]. The latter method works with a small bounded frequency domain instead of the complete time domain by phase-conjugating the residual signal in the frequency domain. This allows for faster computation of images than the TDTE method.

However, FTIM cannot be used in the context of this thesis due to the fact that the temporal variable of the computation is in the Laplace domain which contains the exponential decay window. These topological energy methods in simple media were the next iteration of the TDTG method in that it was possible to localise discontinuities with just a single computation using a direct and adjoint field instead of having to perform multiple iterations. The method was then used to image defects in bounded media [START_REF] Tokmashev | Experimental validation of the topological sensitivity approach to elastic-wave imaging[END_REF][START_REF] Lubeigt | Topological imaging in bounded elastic media[END_REF],

isotropic waveguides [START_REF] Rodriguez | Topological imaging of defects by Lamb wave inspection[END_REF], anisotropic waveguides [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Rodriguez | Topological Imaging of Defects in Anisotropic Plates[END_REF][START_REF] He | Damage identification for composite structures using a cross-correlation reverse-time migration technique[END_REF] as well as highly dispersive media [START_REF] Rodriguez | Guided wave topological imaging of isotropic plates[END_REF]. It also proved to to reasonably effective in dealing with defects which did not have a direct line of sight to the transducer [START_REF] Rodriguez | Locating damages in a complex structure using one or two fixed ultrasonic transducers[END_REF] and in cylindrical structures [START_REF] Nguyen | Defect mapping in pipes by ultrasonic wavefield cross-correlation: A synthetic verification[END_REF].

This chapter deals with the concept of Topological Imaging and presents some interesting experimental and numerical results. The first section introduces and explains the principle of this imaging technique. The next section illustrates the effectiveness of the method by imaging various numerical defects as defined in Chapter 2 taking into account factors such as size, order of the defect, frequency of the excitation signal and number of transducers. The penultimate section deals with the imaging of a defect using measurements of the diffracted field from a Finite Element Simulation on CIVA. The final section demonstrates the potential of the method of Topological Imaging combined with the numerical model developed within the scope of this thesis, to successfully image multiple experimental defects to a high degree of precision.

Principle of the method

Consider a discontinuity at an arbitrary unknown position (r, θ, z) within the structure. A source consisting of one or more transducers illuminates the structure with certain guided modes. These waves propagate over the structure and interact with the discontinuity, causing it to be diffracted in all directions. The reflected diffracted fields measured by the receiver contains information regarding the position and dimensions of the defect. This information used in conjunction with the known incident field can help determine the position of the defect.

Topological imaging is a method based on the cross-correlation of two wave-fields to accurately identify and locate these medium discontinuities. This method requires the computation of two wave-fields, namely the direct (or forward) and the adjoint fields. A cross correlation of the two fields is performed to localise the defects.

The direct field corresponds to the wave-field propagated in a healthy medium. This field (as represented in Figure 2.2) is computed by considering the numerical source (combination of transducer dimensions and signal as defined in Section 1.4.3) to be identical to that of the experimental source used to illuminate the tested medium. This simulation should accurately mirror the real world experiment (as is the case in Section 3.3) albeit in a medium without defects. This is also referred to as the incident field.

The adjoint field is computed using the time reversed residual of the experimentally measured diffracted field (represented in Figure 2.3) as the source for a simulation. The residual is the difference between the measured diffraction response and the solution of the direct problem at the receiver locations. In experimental terminology, this corresponds to baseline subtraction. The residual obtained at each of the transducers is time reversed and used as a source to create the adjoint field.

Finally, a cross correlation between the time reversed direct field and forward propagating adjoint field as described in Equation (4.1) is performed. This physically corresponds to the adjoint field converging at the positions of discontinuities at the same time as the reversed direct field passing over it. The series of images below (refer Figure 4.1)

shows method at work. The idea of the adjoint field is the same as the aforementioned time reversal methods which would require looking for the spatio-temporal maxima. However, correlating it with the direct field directly reveals the positions of the discontinuity. Thus, Topological Imaging is a non-iterative, numerical, post-processing technique that takes advantage of time reversal without having to look for the spatio-temporal maxima to localise defects in complex media.

The method of TI used in this thesis can be distilled down to one equation:

I(z, r) = T 0 u(r, z, t) • v(r, z, T -t) dt , (4.1)
where u(r, z, t) corresponds to the adjoint or time reversed residual field and v(r, z, T -t)

to the direct or forward field propagated backward in time (T -t). The integral over time T corresponds to considering the cumulative effect over all time steps. Taking the absolute value of the calculated image corresponds to taking the spatial envelope [START_REF] Rodriguez | Fast topological imaging[END_REF] which removes the oscillation behaviour leading to a smooth continuous image. Figure 4.2 shows the image obtained by taking the correlation between the adjoint field and time reversed direct field. Figure 4.2a shows the actual computed image that consists of striations that result from the correlation of the peaks and troughs of the two fields at the defect position. Figure 4.2b shows the same computed image but with the addition of having the spatial envelope. This envelope leads to a smoother image, devoid of striations thus making it easier to identify the patch that signifies the defect. The direct field is calculated in a healthy medium devoid of any defects or discontinuities and can have any arbitrary source. However, in testing of tubular structures, axisymmetric and minimally dispersive modes are usually chosen to ensure uniform illumination over the circumference, to reduce the probability of artifacts and to ensure the visibility of all defects in the final image.

The source for the adjoint field is the time reversed residual measured at the receiver positions from the numerical simulation or the experiment. In practice, a belt of transducers is placed circumferentially around the cylinder and act as both transmitter and receiver.

Numerical experiment

This section aims to demonstrate the technique of Topological Imaging to image numerical defects as modelled in Chapter 2. Studies are then performed for various orders of derivatives, frequencies and number of transducers. This is followed by results of some unconventional defect forms.

For all studies presented, the physical and material properties of the cylinder are shown in Table 4.1.

For the first three studies, defects are located at (r , θ, z) = (29 mm, 75°, 450 mm).

The incident source is generated by a circumferentially continuous belt located at z = 0. The transducers used for the computation of the adjoint field have been modelled with diameters of 55 mm. This dimension was chosen to allow for a lighter computation due to the lower circumferential wavenumber requirement. Imaging the numerical defects also expands the understanding of the mechanism of the defect model.

The first section aims to localise the defects of three different sizes introduced in Chapter 2. Imaging these defects reveals the method by which the defect model functions by highlighting the positions of the peaks of the dominant order of the Gaussian that approximates it. In the second section, a defect has been imaged at different frequencies to understand its effect on image sharpness. Next, a study is done to identify the minimum number of transducers to clearly image defects. Following this, two unconventional defects have been modelled and imaged to further explain the bases as described in Chapter 2.

Finally, a defect defined in the simulation tool CIVA [START_REF] Extende | CIVA -Ultrasonic testing[END_REF] has been imaged using only the data recorded from the simulation.

Imaging defects of multiple sizes

In this section, the three defects of radii 2 mm, 4 mm and 6 mm as described in the earlier As mentioned in Section 2.3.3, the defect of radius 2 mm seems to be well simulated using the 1st order derivative and diverges while using the 4th order. In the case of the defect of radius 4 mm, the earlier studies show that using orders 3, 4 and 5 maybe optimum to simulate the defect. Figure 4.4a again shows the 1st order derivative being the dominant secondary source (as in 2 mm case). Moving on to the 3rd order derivative (seen in Figure 4.4b), it would seem like a combination of the base Gaussian function and a second order derivative in the axial directions together minimise the stress field within the contour and the 1st order derivative is less dominant.

The defect studies also show that the 4th order derivative can effectively approximate the defect. Using the 2nd order derivatives for the 6 mm defect seems to approximate the defect similar to the 1st order and has an additional secondary source that acts in unison with the 1st order derivative to further enhance the minimisation of the incident stress. The 4th order derivative behaves similar to the earlier 4 mm case and in this case, the source of this derivative is pronounced.

Studying these cases of topological images of numerical defects approximated by the gaussian function and its derivatives gives us a brief understanding of the mechanism by which the defect model works. Considering the cases in Figure 4.4a and 4.3a, we clearly see the dominant first order derivative of the gaussian behaving as a secondary source. The crest and trough of this function act as secondary sources to minimize the normal incident stress field within the contour of the defect. In these cases the secondary sources of the defect lie outside the defect contour. In other cases, the secondary source is obtained as a combination of multiple dominant orders. However, this method of approximating defects is newly developed and there is still scope for improving and generalising the model based on picking the right base function, the correct order of derivatives, acceptable distances of the peaks in relation to the contour etc.

One other interesting aspect to note is that the method of topological imaging is able to accurately bring out the positions of the secondary sources as well as provide an insight regarding certain aspects of the defect.

Effect of number of transducers on image

This section deals with understanding the effect of the number of transducers on the accuracy and sharpness of the Topological Image. The 4 mm defect is chosen and approximated by the 3rd order derivative. Unlike the earlier case, the incident field is generated by the same transducers that would act as receivers effectively simulating a pulse echo setup. This aides in understanding the advantages and drawbacks for the number of transducers chosen for the study. In each of these cases, the transducers are positioned at equal circumferential distances from each other and the field diffracted by the defect is measured at and back propagated from each of these positions to create the adjoint field.

The images are then constructed, compared and analysed. In an unbounded medium, three non-dispersive wave packets should be sufficient to localise the defect. However, due to the periodic nature of the structure about the circumference, the wave packets of the adjoint field using three transducers have two ways to propagate towards the location of the defect leading to a few prominent patches around the position of the defect, one representing the defect and the others being artifacts. On the whole, this is not an ideal setup for imaging.

The image constructed using four transducers is already accurate in localising the whereabouts of the defect both axially and circumferentially (Figure 4.6b). The formation of patches indicates the presence of a discontinuity but it fails in accurately mapping the approximate form of the defect.

Using six transducers gives an acceptable localisation of the defect (Figure 4.6c) albeit longer axially and thinner circumferentially. There still are some artifacts around the actual position of the defect and there maybe certain cases where this could lead to missing out a defect (For example, a defect present at one of these artifact positions)

Eight transducers seems to be an ideal minimum required to satisfactorily localise a defect for this test configuration (Figure 4.6d). The intensity of the artifacts are negligible in comparison to the patch indicating the defect. The patch is also slightly off centered but this accuracy would be sufficient for most practical purposes at this frequency .

Finally, the image created using twelve transducers is devoid of any visible artifacts around the defect and would more than suffice for most low frequency guided wave nondestructive testing (Figure 4.6b).

There are however some evident artifacts in the region close to the transducers. This is because of the N=12 mode (i.e. 12 wavelengths around the circumference) that exists due to the number of transducers. As explained in Chapter 3.1, it is generally unavoidable to generate a mode with a circumferential order that corresponds to the number of transducers if the cut-off frequency of that mode is near the bandwidth of the excitation spectrum. This behaviour is also evident in the images with 4 or 6 transducers and is persistent in the image due to lack of attenuation in the simulation.

On the whole, the sharpness and accuracy of the image improves with increase in number of transducers.

Effect of frequency on image

This section studies the effect of frequency on the sharpness of the image. To better visualise the effect, the spatial discretisation step has been halved compared to the earlier examples. This leads to the base gaussian function that is proportionally narrower. A defect of radius 4 mm has been approximated by the gaussian and its first derivative. In The series of images (Figure 4.7) clearly show that the sharpness of the image improves as the frequency is increased. This is because the resolution of an image depends on the wavelength of the incident and diffracted field and a smaller wavelength would allow for more localised imaging. It is also seen that for the incident field at 100 kHz, the dominant secondary source is the the 1st order derivative. This is because the secondary source of the defect model depends on the incident field and at 100 kHz, the wavelength (52 mm) is much larger than the diameter of the defect (8 mm).

Alternate defect forms

This section demonstrates the imaging of certain unconventional defect. The examples are specifically chosen to emphasise certain ideas to keep in mind while deciding the base functions and derivatives to be used to create an ideal orthonormal basis for a given defect.

The first example is the elongated defect as mentioned in Chapter 2. This defect is centered at θ=0°and z=0.45 m, having a circumferential length of 100 mm and axial width The simulation is performed with the same specifications as in the case above (5cycle, 200 kHz, 24 transducers). The topological image of this setup is shown in Figure 4.9b. This image not only verifies the functioning of the method of TI but also validates the method of modelling defects.

CIVA

In this final section dealing with numerical defects, a defect has been modelled in the simulation software CIVA [START_REF] Extende | CIVA -Ultrasonic testing[END_REF]. CIVA is an extensively used analysis and simulation software in the field of NDT. It has a comprehensive list of functionalities involving bulk and guided waves, allowing the user to simulate processes such as pulse echo or TOFD, using a variety of probes (Phased Arrays, EMATS etc.) to simulate a large range of flaws using techniques frequently used in the NDT world. The numerical tool created during the course of this thesis would be implemented as a plug-in for the simulation of guided waves in cylinders in CIVA.

A simulation was performed on CIVA with properties as shown in Table 4. This shows cross platform validation of the tool that is not only able to accurately replicate the incident field generated by the simulation tool CIVA, but also use this information and the residual signals to localise the defect.

Physical experiment

This section deals with the TI of experimental discontinuities. Throughout the course of the studies that follow, the words defect and discontinuity have been used interchangeable. This is done due to the way in which the experimental defects were created. During the course of certain experiments, it was found that a thin foil of Metglas ® (generally used in magnetostriction) affixed onto the surface of the cylinder using a shear couplant and kept in place with scotch tape was extremely reflective. This behaviour worked to our advantage as the "defect" could be cut to the required size and affixed at different positions without physically damaging the structure. This also allows for the extraction of pure experimental residuals containing only the signature of the defect and limited only by the sensitivity of the equipment.

The experimental setup used is as described in Section 3.1. It consists of a transmitter and a receiver belt. The position of the receiver belt is at z = 0 cm. The transmitter belt is placed behind (in relation to the side containing the defects) located at z = -3 cm.

A 5-cycle toneburst signal with a central frequency of 100 kHz has been used. Unless otherwise specified, the L(0, 2) longitudinal mode has been adopted as the source for defect localisation. At 100 kHz, the wavelength of the L(0, 2) guided wave is around 52.8 mm.

Effect of baseline subtraction on imaging sensitivity

This study presents the Topological Images obtained using theoretical and experimental residuals and demonstrates the versatility of the method. Three square defects of sides In theory however, as the principle mode (the minimally dispersive L(0, 2) in our case) propagates away, there is no deflection of the signal from the mean amplitude except for those caused due to any discontinuities (defects or back wall). To avoid any complications due to the reasons mentioned above, the baseline of the theoretical signal is taken to be zero. This would mean that the theoretical residual is taken to be the measured experimental signal.

Topological Imaging based on defect size using theoretical residuals

The following topological images are obtained using theoretical residuals. Figure 4.16a

shows the TI of the defect of size 15 mm (around 1 /3rd the wavelength). The defect can be fairly accurately distinguished from the background noise. There is also some coherent noise present at the beginning of the image. This is due to the interaction between the circumferentially propagating direct field of the simulation and the backward propagated adjoint field and is persistent due to the absence of attenuation in the simulation.

Moving on to the image of the defect of length 10 mm (around 1 /5th the wavelength) Figure 4.16b shows that the patch that indicates the defect is still distinguishable from the background coherent noise. However, the intensity of the patch is slightly diminished compared to the earlier case of the 15 mm defect. This could be attributed to the fact that the smaller defect has a smaller signature within the diffracted signal compared to the baseline.

Observing the topological image of the 5 mm defect ( 1 /10th the wavelength), it is clear that we are approaching the lower limit of the defect size that can be imaged using only the theoretical residual. The defect patch is still visible, however its intensity is comparable to the noise around it. This is due to the fact that the signature of the defect has an amplitude very similar to the baseline, almost masking it. In the examples shown above, the theoretical residual helps locate the defect with varying degrees of intensities. Regardless of the amplitude of the signal diffracted by the defect, during the creation of the topological image by the method of correlation of the time reversed direct field with the adjoint field, the time reversed residuals emitted at each transducer, converge at the defect location at the exact moment the direct field passes over it. It is true that the backward propagating direct field also interacts with the adjoint field generated by all other wave packets and disturbances. However, the amplitude of the correlation of the fields at the defect location is far greater than that around it. This is reminiscent of the spatio-temporal maxima technique used in earlier Time Reversal methods with the added step of not having to search through the data for this maxima.

Topological imaging of special defects with experimental residual

In the earlier study, the theoretical residual successfully illuminates the defect in all cases with little to no hassle. This is simply due to the fact that in subtracting the baseline from the signal, all that remains is the signature of the defect. In this section, a couple of examples have been discussed that demonstrates the scope of using experimental residuals for TI. This technique would fall under the purview of Structural Health Monitoring (SHM).

This study first explores the imaging of extremely small perturbations that approach the sensitivity of the equipment. Next, defects that are completely masked are made visible by utilising the experimental residual.

Imaging of extremely small perturbations

To test the limits of the system, extremely small perturbations were applied in the form of the tip of a blunt pencil and of a compass. These objects were placed on the surface and sufficient force was applied so as to avoid slipping. The tip of the pencil was measured to have a diameter of around 1.37 mm (close to 1 /40th the wavelength and that of the compass had a diameter of around 0.35 mm (close to 1 /150th the wavelength).

Figure 4.17 compares the Topological Images of the pencil point defect using the theoretical and experimental residuals. From Figure 4.17a, it is clear that it is no longer possible to image perturbations of such small dimensions purely by using the theoretical residual (as defined for this chapter). The signature of the defect is completely hidden within the coherent noise of the baseline signal, rendering it effectively invisible in the image. On the other hand, using the experimental residual still allows the defect to be imaged as can be seen in Figure 4.17b. It should be noted that this sensitivity in imaging may only be realized when there is very little to practically no change between the measurements of the specimen before and after the introduction of the discontinuity. This includes the positioning of the transducers, the strength of couplant, temperature (which can change not just the velocity of the guided wave but also the viscosity of couplant), equipment used for signal acquisition, etc. The sensitivity for such an image is so delicate that subsequent topological images attempting to use the baseline signal data from an earlier date failed to reproduced an image of the quality as presented in Figure 4.17b. Nevertheless, the result by itself is interesting and goes to show the capabilities of using the experimental residual for topological imaging.

As a final test of sensitivity of the equipment, Figure 4.18 below shows the topological image of the compass tip as the perturbation using the experimental residual. As was the case earlier with the 5 mm square defect (Figure 4.16c) the patch signifying the defect may be mistaken for noise due to its intensity. This maybe considered as the lowest limit of defect detection with the equipment, frequency bandwidth, number of transducers and acquisition system that was used. On the other hand, this study opens up the possibility of verifying the limits of the method with other defects such as corrosion, slits, pits, through holes, etc.

Imaging of hidden defects

In the above section, the potential of the experimental residual has been explored.

In this section, an additional use has been examined. In pipeline testing, situations can be anticipated with obstacles such as support structures or welds. In such cases, the reflected signal may contain the information of a defect (if present), having an amplitude far smaller than the reflection from the support or weld. It is also possible that the defect is masked by the obstacle. Usually, a part of the signal is transmitted and a part reflected.

The transmitted portion of the wave may interact with a defect beyond the obstacle and get reflected back to the transducer. Due to losses in energy due to material attenuation and an impedance mismatch, the reflected secondary signal (from the defect) may have an amplitude comparable to coherent noise. It would suffice to store the signals of the undamaged structure beforehand to use as the baseline.

One such example has been shown below. A strip of Metglas ® is wrapped around the pipe to simulate such an obstacle. A square patch of Metglas ® with side 15 mm has been affixed 100 mm behind the circumferential strip to simulate the hidden defect. Clearly, the intensity of the continuous strip obscures the defect patch behind it. This is To image the hidden defect, the experimental residual can be computed considering the baseline signal that already contains the signature of the obstacle (circumferential strip in our case). This baseline could be obtained either from a previous experimental measurement containing the obstacle or from an accurate numerical computation that simulates the effect of the obstacle. Figure 4.20 shows the image of the defect obtained using the experimental residual. The defect is clearly visible and the signature of the belt has been completely eliminated.

On the whole, using the experimental residual is a powerful tool in revealing extremely small perturbations as well as defects which may otherwise be hidden by necessary obstacles having much larger signatures. However, it is to be noted that the method by which the defects were created experimentally allows for an accurate determination of the experimental residual. Obtaining a perfect experimental residual in real world NDT applications would be limited. Although, taking the product of the two (or more) images would work in some cases it may lead to certain unwanted behaviour which is made evident in the second case.

Case 2: Imaging of defects at azimuths 0°and -90°I n this case, the second defect is positioned at an angle of -90°. This has been done

to show the shortcoming of using the product method to improve the image. As in the earlier case, the image (as seen in Figure 4.24) obtained by the L(0, 2) mode exposes the two defects and maps out their axial and circumferential positions. On imaging the same setup with the F (1, 3) mode however, the second defect is rendered invisible (seen in Figure 4.25). This is due to the shape of the wavefront of the F (1, 3) mode. This mode (and all of the N = 1 modes) have two nodes when viewed around the circumference. In our case, the nodes are positioned at the 90°and -90°a

zimutal positions leading to a lack of illumination at these zones over the entire length of the cylinder. This inturn leads to an almost zero reflection from this defect and hence its signature is not measured at the transducers. Taking the pointwise product of the intensities in such cases would lead to the creation of a Topological Image that is incomplete as can be seen in Figure 4.26. Here, there is an improved clarity of the first defect but the second defect is not imaged. An alternate method to improve the accuracy of the image would be to take the pointwise mean of the intensities. Figure 4.27 shows the image obtained by takint the pointwise mean of the images obtained using using the two modes (N = 0 and N = 1).

Evidently, the clarity of the first defect has been slightly improved and the information of the second defect has not been completely lost. However, this may not be the ideal method to improve the clarity and of images as taking a mean of multiple images may lead to such defects being lost within the coherent noise of the images. Further studies maybe carried out to obtain a model of combining the images obtained from multiple modes.

Summary of Topological Imaging

This chapter discusses the Topological Imaging of various kinds of defects. First, the principle of TI was introduced and explained. Next, defects created by the modelling technique as discussed in Section 2.3 were imaged. This also opened up a window into understanding the method by which the model functions. The chapter goes on to discuss a couple of studies based on the effect of number of transducer and frequency on the sharpness of the image concluding that in general a higher frequency and a greater number of transducers would lead to a sharper and more localised image. Modelling and imaging of certain alternate numerical defect forms were then explored showing the effectiveness of the method in obtaining the general shape of the defect as well as to identify multiple defects.

The next part discusses the imaging of a defect modelled in CIVA, using only the output of the simulation software (signals at the receiver). The imaging of this defect showed that the simulation tool and method of TI is cross compatible with external FE simulation software.

The chapter then examines the imaging of experimental defects. The studies explore the effect of baseline subtraction on the quality of imaging defects. The first study uses theoretical residual and shows the sharpness of images based on the dimension of the defect. Here defects as low as 1 /10th the wavelength were successfully imaged. Applying experimental baseline pushes the lower limit of imaging of extremely small defects (in orders down to 1 /40th of the wavelength). It also allows for the imaging of masked defects.

Finally a couple of studies show the imaging of multiple defects and improving image sharpness using multiple modes. In particular, two modes were used for imaging and two methods of image combination were examined. A reasonable conclusion is that in most cases, it is probably a safer option to use axi-symmetric modes for imaging with multiple modes as this allows to cover the entire circumferential extent of the tubular structure.

In all cases discussed, the defects have been successfully localised in the axial and circumferential extents. Situations could be envisaged where it may be important to know the radial position of the defect. For example, if a corrosion is at the inner or outer surface of the structure or a defect is located within the material. This would call for higher order modes that preferentially illuminate or shadow certain radial positions of the structure.

These studies maybe a part of future work using the models and techniques discussed in this thesis.

In general, this chapter shows that Topological Imaging is an effective method of accurately and robustly localising numerical and experimental defects in complex media.

Conclusions

This thesis deals with the development and implementation of a novel numerical tool for the simulation of guided wave propagation in tubular structures as well as a new technique to approximate numerical defects. The method of Topological Imaging is then introduced and effectively applied for the localisation of numerical and experimental defects. Chapter 2 introduces the method of approximating numerical defects in layered media. These defects are described as a discontinuity between the different layers, simulating a delamination. The defect is considered as a secondary source that creates a diffracted field such that the incident field is completely negated within its contour. It is approximated as the sum of the responses of a Gaussian function and its derivatives. There are also possibilities of modelling 3d defects such as notches and inclusions and work has already begun towards the same. One other possibility is to broaden the model to include structures with liquid within the pipe (Ex: Oil pipeline) or on the outside (Ex:

Pipeline submerged under the sea) taking into account leaky guided waves. There may also be potential to further the model to include pipe with bends or even structures with varying cross sections such as wind turbine support structures or blades.

  Mots clés : Ondes Ultrasonores, Contrôle Non Destructif, Imagerie Topologique, Guide d'ondes, Contrôle de Santé Intégré. Introduction 1 Transient Wave Propagation in a Stratified Tubular Medium 1.1 Notations and domains for calculation . . . . . . . . . . . . . . . . . . . . 1.2 Modelling of wave propagation equations . . . . . . . . . . . . . . . . . . . 1.2.1 Basic equations in the physical space . . . . . . . . . . . . . . . . . 1.3 Equation to be solved in the (r, n, k, s)-domain . . . . . . . . . . . . . . . 1.3.1 General case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Transversely isotropic medium . . . . . . . . . . . . . . . . . . . . . 1.4 Software implementation of the model . . . . . . . . . . . . . . . . . . . . 1.4.1 Discrete calculation (creation of grid) . . . . . . . . . . . . . . . . . 1.4.2 Multiprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3 Definition of source . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.4 Determining the length of the calculation . . . . . . . . . . . . . . . 1.5 Summary of model creation . . . . . . . . . . . . . . . . . . . . . . . . . .

4 . 2 . 5 2 . 4

 42524 CIVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Physical experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Effect of baseline subtraction on imaging sensitivity . . . . . . . . . 4.3.2 Multi-defect imaging with multiple modes . . . . . . . . . . . . . . 4.4 Summary of Topological Imaging . . . . . . . . . . . . . . . . . . . . . . . Conclusions ix List of Figures 1.1 Cylindrical coordinates system. . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Multilayered infinite "n" layered cylinder. . . . . . . . . . . . . . . . . . . . 1.3 Snapshots over time showing the displacement field radiated by a localized source to demonstrate circumferential periodicity. The waves radiate from a source located at z = 0 m, θ = 0.8 π and is visible from the θ = -π edge starting from t = 30 µ s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Snapshots over time showing the displacement field radiated by a localized source to demonstrate axial periodicity. The waves radiate from a source located at z = -0.38 m, θ = 0 and is visible from the z = 0.4 m edge starting from t = 25 µs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Graphic showing the propagation of two principle wave packets over time (in the Y axis) over the axial length (centered at 0) for θ = 0 and extending from -z max to z max . The region of valid computation is in blue and corrupted region of computation signals due to periodicity is in red. . . . . 1.6 Schematic showing the domain of interest, domain to be computed, two modes with velocities v f and v s and transducer position at z = 0. . . . . . 2.1 The 3d and its respective 2d representation of the defect in the grid and its position at the interface of the two layered cylinder. . . . . . . . . . . . 2.2 Defect at axial distance z, circumferential position θ and radial position r I being radiated by an axisymmetric incident field generated by the transducers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Field diffracted by the defect, computed as a secondary source to cancel the incident field (represented as a radial stress vector) within its contour. x Representation of dependency conditions of new elements to an orthonormal basis of dimension α = 3. . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 1D representation of the Gaussian function and its derivatives. . . . . . . . 2.6 2d representation of the Gaussian function and its derivatives placed on the grid. Dashed circular line represents the contour of the defect. . . . . . 2.7 Visual representation of the ideal and computed cases of the conditions that define the defect model. Image shows the incident field vector projected onto the complete (left) and incomplete basis (right) and its components that define the diffracted and total fields. . . . . . . . . . . . . . . . . . . . 2.8 Gaussian function and derivatives overlaid on the 6 mm defect contour showing the positions of the peak of the functions in relation to the defect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 The spectrum of the total field measured at the center of the 6 mm radius defect for derivatives of different orders. The dotted line denoting 10 % of the maximum amplitude below which the defect approximation may be considered effective. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Snapshot at time 45 µs showing the axial component of the sum of the incident field (created by a circular source located at the center) and the diffracted field by the 6 mm radius defect for various order of derivatives. . 2.11 Gaussian function and derivatives overlaid on the 4 mm defect contour showing the positions of the peak of the functions in relation to the defect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 The spectrum of the total field measured at the center of the 4 mm radius defect for derivatives of different orders. The dotted line denoting 10 % of the maximum amplitude below which the defect approximation may be considered effective. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Snapshot at time 45 µs showing the axial component of the sum of the incident field (created by a circular source located at the center) and the diffracted field by the 4 mm radius defect for various order of derivatives. . xi 2.14 Gaussian function and derivatives overlaid on the 2mm defect contour showing the positions of the peak of the functions in relation to the defect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 The spectrum of the total field measured at the center of the 2 mm radius defect for derivatives of different orders. The dotted line denoting 10 % of the maximum amplitude below which the defect approximation may be considered effective. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Snapshot at time 45 µs showing the axial component of the sum of the incident field (created by a circular source located at the center) and the diffracted field by the 2 mm radius defect for various order of derivatives. . 2.17 Images showing all the derivatives required to approximate an elongated defect to cover the area within the contour. . . . . . . . . . . . . . . . . . . 2.18 Images showing the placement of the Gaussian basis functions to partially and more effectively approximate an irregularly shaped. . . . . . . . . . . . 3.1 Direction of transducer polarisation. . . . . . . . . . . . . . . . . . . . . . . 3.2 Transducers placed in the ring to maintain spacing. . . . . . . . . . . . . . 3.3 Schematic representation of the experimental test setup. . . . . . . . . . . 3.4 Coefficients of emitters (e i ), transmission (t k ) and receivers (r j ) in relation to their relative positions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Un-calibrated and calibrated signals at 8 receivers at a distance of 4.7 m from the transmitter belt. . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 6

 6 Setup to study circumferential attenuation due to the addition of transducers. The wave propagates from the emitter (orange) to the receiver (blue).The inactive transducers (grey) are added sequentially to measure their additional impedance at the surface. . . . . . . . . . . . . . . . . . . . . . 3.7 Signals measured at the receiver showing increased attenuation with the addition of transducers circumferentially. . . . . . . . . . . . . . . . . . . . 3.8 Measurement of velocity with multiple through thickness reflections. . . . . 3.9 Signal of multiple reflections of the longitudinal bulk wave over the thickness. xii 3.10 Frequency spectrum of reflected signal as obtained in Figure 3.9. . . . . . . 3.11 Spectrum of frequency spectrum (as obtained in Figure 3.10) obtained by taking a Fourier transform of the absolute value of the spectrum giving the time between successive reflection. . . . . . . . . . . . . . . . . . . . . . . . 3.12 Normalised displacement field signal for Gaussian sources of different widths measured at 50 cm. Figure also shows the packets that represent various modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Normalised displacement field signal for gaussian sources of different lengths measured at 50 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Amplitude distribution of two modes around the circumference of the pipe. These may represent stress or displacement amplitudes. . . . . . . . . . . . 3.15 Frequency-axial-wavenumber dispersion curve showing the fundamental modes from circumferential wavenumber 8 to 15 in relation to 100 kHz central frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Frequency-axial-wavenumber spectrum of the computed field overlaid on relevant dispersion curves showing a good agreement between the calculated and theoretical solutions. . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Comparison of the computed (orange) and experimentally measured (blue) group velocities of the axisymmetric (N = 0) mode [L(0,2)] measured at 50 and 100 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Comparison of the computed (orange) and experimentally measured (blue) group velocities of a non axisymmetric (N = 1) mode [F(1,3)] measured at 50 and 100 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Comparison of the computed (orange) and experimentally measured (blue) group velocities of a non axisymmetric (N = 4) mode [F(4,2)] measured at 50 and 100 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Spectrum of the axial displacement component of the computed field obtained by exciting the L(0,2) mode. . . . . . . . . . . . . . . . . . . . . . . 3.21 Spectrum of the axial displacement component of the computed field obtained by exciting the F(1,3) mode. . . . . . . . . . . . . . . . . . . . . . . xiii 3.22 Spectrum of the axial displacement component of the computed field obtained by exciting the F(4,2) mode. . . . . . . . . . . . . . . . . . . . . . . 3.23 Comparison of theoretical group velocities (solid lines) of the three studied modes [L(0,2), F(1,3) and F(4,2)] with various trials of experimentally obtained group velocities showing good agreement. . . . . . . . . . . . . . 4.1 Snapshots over time showing the forward propagated adjoint field (top), the backward propagated direct field (middle) and the topological image of the defect forming by the correlation of the two aforementioned fields (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Comparison of the topological images obtained by the cross correlation of the adjoint and time reversed direct fields to identify the defect (orange contour). The image on the left is the actual image and that on the right is computed by taking the spatial envelope of the first image. . . . . . . . . 4.3 Topological images of the 2 mm defect modelled using different orders of the base gaussian function. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Topological Images of 4mm defect modelled using 1st, 3rd, 4th and 5th orders. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Topological Images of 6mm defect modelled using 2nd and 4th orders. . . . 4.6 Topological images formed using a number of transducers for adjoint field computations showing that increasing the number of transducers improves accuracy and sharpness of the image. . . . . . . . . . . . . . . . . . . . . . 4.7 Series of images showing the topological image of the same defect imaged at four different frequencies ; 100, 150, 200 and 250 kHz. . . . . . . . . . . . 4.8 Topological image of a defect 5.5 mm wide axially and 100 mm long circumferentially. The orange contour represents the delamination defect. . . 4.9 Figures showing the single contour to describe axially and circumferentially separated defects (left) and the topological image (right) identification of the two defects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 4.10 Schematic showing the axial and circumferential position as well as the angular extent of the through-wall defect as defined in the simulation software CIVA. The rectangles in red show the circumferential positions of the transducers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Snapshot of the incident axial displacement field at the surface of the pipe at 140 µs. The incident field is generated by 8 equiangular transducers exciting a 15-cycle tone burst signal at 100 kHz. . . . . . . . . . . . . . . . 4.12 The signals of the diffracted field received at the 8 transducer positions measured in CIVA. These signals, temporally inverted are used as the input to create the adjoint field in the numerical tool. . . . . . . . . . . . . . . . 4.13 Topological image of the defect modeled in CIVA and localised by the numerical model using only the reflected signals obtained from the CIVA simulation shown in Figure 4.12. . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Experimental residuals of signals diffracted from 5 mm (blue), 10 mm (orange) and 15 mm (green) defects measured at a single receiver obtained by experimental baseline subtraction. . . . . . . . . . . . . . . . . . . . . . . . 4.15 Theoretical residuals of signals diffracted from 5 mm (blue), 10 mm (orange) and 15 mm (green) defects measured at a single receiver computed with the numerical simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Topological images of the 15 mm, 10 mm and 5 mm Metglas ® defects showing the decrease in amplitude of the patch representing the defect, signifying the decrease in sensitivity based on defect dimensions. . . . . . . . . . . . . 4.17 Comparison of the topological images of the pencil tip used to mimic a small defect obtained using the theoretical and experimental residuals. . . 4.18 Topological image of a compass tip used as an extremely small diffractor computed using the experimental residual. . . . . . . . . . . . . . . . . . . 4.19 Topological image of a defect behind a continuous Metglas ® belt without experimental baseline subtraction showing the inability to image the defect. 4.20 Topological image of a defect behind a continuous Metglas ® belt using experimental baseline subtraction clearly showing defect. . . . . . . . . . . xv Introduction Soon after the culmination of the industrial revolution in the 1860s, oil and natural gas were used as primary sources of fuel to run different industries and machinery. Over time, this led to the development of vast networks of pipelines running over 3 million kilometers worldwide. Due to its ubiquitous presence and the nature of items they transport, it is of utmost importance to ensure their integrity and health. This required a solution which came in the form of Non Destructive Testing (NDT) and Structural Health Monitoring (SHM).

  to improve depth estimation accuracy and to use probablistic based imaging algorithms. This PhD deals with the development, verification and implementation of a new mathematical model to simulate Guided wave propagation in Tubular structures and to localize numerical and experimental defects with reduced computation times. The first chapter develops the mathematical model which is based on solving the problem in the Laplace and Fourier domains. The solutions are given in the form of modified Bessel's functions of the first and second order. The implementation of the model is then discussed. The second chapter details a method of modelling defects as secondary sources. These defects are approximated by Gaussian functions and their derivatives and a few defect studies are performed to understand the behaviour of the model. The third chapter deals with the experimental verification of the numerical model. It introduces the experimental setup and a method to normalise the transducers. It then discusses the parameters of the numerical model that most accurately simulate the experimental setup. The model is then validated with experimental results using different modes. Finally the capabilities and limitations of the model are discussed.

Figure 1 . 1 :

 11 Figure 1.1: Cylindrical coordinates system.

Figure 1 . 2 :

 12 Figure 1.2: Multilayered infinite "n" layered cylinder.

c 11 c 11 -2 c 66 c 13

 1113

Table 1 . 2 :

 12 Diamond product of the cylindrical basis vectors, for transversely isotropic media, where: (b a) = (a b) T .

  zero six-dimensional vector a is searched such that traction at the outer and inner radii (r max and r min respectively) of the pipe wall are zero, i.e., Σr (r min ) = 0 and Σr (r max ) = 0(1.26)Finally the standard dispersion equation is given by: min ) Σr2 (r min ) Σr3 (r min ) Σr4 (r min ) Σr5 (r min ) Σr6 (r min ) Σr1 (r max ) Σr2 (r max ) Σr3 (r max ) Σr4 (r max ) Σr5 (r max ) Σr6 (r max )

  periodicity of the real domain is observed due to the periodic nature of the wavenumber domain of the Fourier transform. This periodicity is demonstrated by two simple simulation examples with extreme placement of the source. For these simulations, a domain with axial length extending from -0.4 m to 0.4 m is considered. The radius of the pipe is 0.3 m with a thickness of 2 mm. The axial and circumferential widths of the transducer source are 20 mm and 10 mm respectively and the signal is a 5 cycle tone burst with a central frequency of 100 kHz.

Figure 1 . 3 :

 13 Figure 1.3: Snapshots over time showing the displacement field radiated by a localized source to demonstrate circumferential periodicity. The waves radiate from a source located at z = 0 m, θ = 0.8 π and is visible from the θ = -π edge starting from t = 30 µ s.

Figure 1 . 4 :

 14 Figure 1.4: Snapshots over time showing the displacement field radiated by a localized source to demonstrate axial periodicity. The waves radiate from a source located at z = -0.38 m, θ = 0 and is visible from the z = 0.4 m edge starting from t = 25 µs.

Figure 1 . 5 represents

 15 the field propagating at spatial co-ordinates (r = r rmax , θ = 0, z) with respect to time (represented in the Y axis), i.e. field propagation represented in two dimensions, where the X axis represents a line from -0.4 m to 0.4 m for a fixed θ position (θ = 0) and the Y axis represents time. The graphic shows the evolution of two mode packets of different velocities propagating outward from z = 0 (source position). The region in blue is the valid region of computation. The region in red contains a corrupted signal of the faster mode due to the periodicity of the domain inherent to the Fourier domain computation. This does not correspond to the actual physics of propagation. The second slower wave packet has not had time to illuminate the entire physical domain before the faster wave packet traverses the edge of the domain and corrupts the wave field. This phenomenon is extremely important to consider especially to obtain the Topological Image of the defect.

Figure 1 . 5 :

 15 Figure 1.5: Graphic showing the propagation of two principle wave packets over time (in the Y axis) over the axial length (centered at 0) for θ = 0 and extending from -z max to z max . The region of valid computation is in blue and corrupted region of computation signals due to periodicity is in red.

Figure 1 . 6 :

 16 Figure 1.6: Schematic showing the domain of interest, domain to be computed, two modes with velocities v f and v s and transducer position at z = 0.

  The incident field generated by the transducers (each modelled as shown in Section 1.4.3) represented in Figure2.2 propagates towards the defect positioned at an arbitrary axial position z, azimuthal position θ and radial position at the interface r I .

  (a) 3d representation of the layered cylinder and the defect at the interfacial surface placed on the grid of discretisation points (b) Unwrapped view of the cylindrical domain showing the defect contour and discretisation points.

Figure 2 . 1 :

 21 Figure 2.1: The 3d and its respective 2d representation of the defect in the grid and its position at the interface of the two layered cylinder.

Figure 2 . 2 :

 22 Figure 2.2: Defect at axial distance z, circumferential position θ and radial position r I being radiated by an axisymmetric incident field generated by the transducers.

Figure 2 . 3 :

 23 Figure 2.3: Field diffracted by the defect, computed as a secondary source to cancel the incident field (represented as a radial stress vector) within its contour.

Figure 2 .

 2 Figure 2.4 shows two cases of the intermediate step of the addition of a third basis vector ψ 3 to a pre-existing orthonormal basis consisting of the two basis functions ψ 1 and ψ 2 . This newly added basis function is created by first taking the new free function θ 3(in green) and subtracting from it its projection (in red) on the subspace containing the pre-existing orthonormal basis. This creates a third basis function which is orthogonal to the pre-existing basis. The next step of normalizing this newly added basis function is done by taking the inverse of the norm of the rejection (in blue).

Figure 2 . 4 :

 24 Figure 2.4: Representation of dependency conditions of new elements to an orthonormal basis of dimension α = 3.

Figure 2 .

 2 [START_REF] Blouin | Improved resolution and signal-to-noise ratio in laser-ultrasonics by SAFT processing[END_REF] shows the different functions which are used together to try and simulate the defect. The peaks of the first and second orders lie outside the contour of the defect.Ideally, a Dirac delta function should be used to describe point size defects. However, representing a Dirac accurately would require an infinitesimally small step size of the grid which is impractical in terms of the size of the computation. A Gaussian function and its derivatives, that approximates a Dirac delta function and its derivatives, is chosen due to the discrete nature of the problem and to satisfy the Nyquist-Shannon criterion (axial wavenumber less than π/δz and circumferential wavenumber less than π/δθ, where δz and δθ denote the axial and circumferential discretization step sizes).

Figure 2 . 5 :

 25 Figure 2.5: 1D representation of the Gaussian function and its derivatives.

Figure 2 .

 2 Figure 2.6 shows the form of the 2d Gaussian function its derivatives, used to create the orthonormal basis, overlaid on the grid and defect contour. The Gaussian functionand its family of derivatives correspond to the family of free functions θ as mentioned in Section 2.1.2. The shape function φ of the defect can also be defined as any arbitrary form inscribed within the shape of delamination and the minimisation is carried out within the area described by the contour of the function.

  (a) Gaussian Function (b) First order derivative in z (c) First order derivative in a (d) Second order derivative in z and a

Figure 2 . 6 :

 26 Figure 2.6: 2d representation of the Gaussian function and its derivatives placed on the grid. Dashed circular line represents the contour of the defect.

Figure 2 . 7 :

 27 Figure 2.7: Visual representation of the ideal and computed cases of the conditions that define the defect model. Image shows the incident field vector projected onto the complete (left) and incomplete basis (right) and its components that define the diffracted and total fields.

  uniquely as a linear combination of the vectors. The two indices represent the order of the partial derivatives of the Gaussian function in the axial and circumferential directions respectively. The objective of the model is to define secondary sources present in the vicinity of the defect which work together in a way to minimise the incident stress field within the contour of the defect.Consider an isotropic aluminum cylinder of inner and outer radius 28 mm and 30 mm respectively. The length of the cylinder (188.5 mm) has been chosen to be equal to the circumference. A circular source of radius 15 mm centered at (0, 0) placed on the surface, exerts a axial force of 5 cycles at 250 kHz. The "delamination" defect is placed at the midplane of this setup at radius 29 mm and centered at axial position 47 mm and angular position 90°. This placement was chosen so as to take into account the effect of the axial and circumferential wave vectors of the computation. The step size in the axial direction is 2.5 mm and in the circumferential direction is 2.35 mm (circumferential wavenumber of 40) to try and maintain a constant scale in the two directions. The coefficient of the Gaussian shape function is chosen as 2/9 so as to satisfy the Nyquist Shannon criterion.

( 2 . 3 . 1

 231 red point) has been plotted as a function of frequency. The study aims to understand the effectiveness of the different order derivatives of the base Gaussian function in cancelling out the incident stress field within the contour for defects of varying dimensions. Snapshots of the total field have been shown to visualise the response of the defect to the incident field. The study ends by describing the conditions to model a defect with this technique. Case 1: radius of defect 6 mm This section deals with a defect of radius 6 mm. The contour of the defect encloses 21 points in the grid as seen in the series of images in Figure 2.8. Due to the dimension of the defect in comparison to the Gaussian function, higher order derivatives are required to cover the area within the contour. The images also show the Gaussian function and its derivatives in relation to the defect contour.

Figure 2 . 8 :

 28 Figure 2.8: Gaussian function and derivatives overlaid on the 6 mm defect contour showing the positions of the peak of the functions in relation to the defect.

Figure 2 .

 2 [START_REF] Burch | Ultrasonic Synthetic Aperture Focusing Using Planar Pulse-Echo Transducers[END_REF] show the defect of radius 4 mm placed within the grid overlaid on the Gaussian function and the different derivatives. The contour of the defect encloses 9 points in the grid.

  (a) 1st order derivative approximation (b) 2nd order derivative approximation (c) 4th order derivative approximation (d) 5th order derivative approximation

Figure 2 . 10 :

 210 Figure 2.10: Snapshot at time 45 µs showing the axial component of the sum of the incident field (created by a circular source located at the center) and the diffracted field by the 6 mm radius defect for various order of derivatives.

Figure 2 . 11 :

 211 Figure 2.11: Gaussian function and derivatives overlaid on the 4 mm defect contour showing the positions of the peak of the functions in relation to the defect.

  (a) 1st order derivative approximation (b) 3rd order derivative approximation (c) 4th order derivative approximation (d) 5th order derivative approximation

Figure 2 . 13 :

 213 Figure 2.13: Snapshot at time 45 µs showing the axial component of the sum of the incident field (created by a circular source located at the center) and the diffracted field by the 4 mm radius defect for various order of derivatives.

Figure 2 .

 2 Figure 2.14 shows the defect of radius 2 mm within the grid overlaid on the Gaussian function and its derivatives. The contour of the defect encloses only 1 point in the grid.

Figure 2 . 14 :

 214 Figure 2.14: Gaussian function and derivatives overlaid on the 2mm defect contour showing the positions of the peak of the functions in relation to the defect.

Figure 2 . 15 :

 215 Figure 2.15: The spectrum of the total field measured at the center of the 2 mm radius defect for derivatives of different orders. The dotted line denoting 10 % of the maximum amplitude below which the defect approximation may be considered effective.

Figure 2 .

 2 Figure 2.15 shows the total field measured only for the first 3 orders of derivatives. The response of the 4th and 5th derivatives lie beyond the specified range.

Figure 2 . 16 :

 216 Figure 2.16: Snapshot at time 45 µs showing the axial component of the sum of the incident field (created by a circular source located at the center) and the diffracted field by the 2 mm radius defect for various order of derivatives.

Figure 2 . 17 :

 217 Figure 2.17: Images showing all the derivatives required to approximate an elongated defect to cover the area within the contour.

  (a) Irregularly shaped defect partially approximated by a Gaussian basis function. (b) Irregularly shaped defect approximated better with two Gaussians basis functions.

Figure 2 . 18 :

 218 Figure 2.18: Images showing the placement of the Gaussian basis functions to partially and more effectively approximate an irregularly shaped.

  and are shear polarised as shown in Figure3.1. The contact surface of the transducer has dimensions 3 mm×13 mm. Eight of these transducers are placed equi-angularly around the cylinder so as to have the direction of polarisation of the transducers along the axial direction. Failing to maintain this axial symmetry leads to the generation of unwanted modes. A ring was designed and 3D printed to maintain the angular separation and to ensure repeatability. This ring can be slid into position and the transducers slotted into the pockets. The setup is mounted as shown in Figure3.2.

Figure 3 . 1 :Figure 3 . 2 :

 3132 Figure 3.1: Direction of transducer polarisation.Figure 3.2: Transducers placed in the ring to maintain spacing.

Figure 3 . 3 :

 33 Figure 3.3: Schematic representation of the experimental test setup.

Figure 3 . 4 :

 34 Figure 3.4: Coefficients of emitters (e i ), transmission (t k ) and receivers (r j ) in relation to their relative positions.

2 )

 2 Summing the equation (3.1) over all i and j values and applying the simplification as shown in Equations (3.2), ij

Figure ( 3

 3 .5) below shows the effect of this calibration. The two images show the combined signals at all 8 receivers in response to an axi-symmetric source. The first of the two images shows the situation where the calibration has not been applied. The variation in amplitudes at each receiver is visible in the inset figure. There are also trailing signals which for our purpose is considered as noise. The second image shows the signals at the receivers after the calibration step.Here, the amplitudes at each receiver have the same amplitude with a very negligible variation. The trailing noise is also eliminated. On the whole, this technique of calibration is effective in normalising the transmitters and receivers, hence creating a pure axi-symmetric mode.

Figure 3 . 5 :

 35 Figure 3.5: Un-calibrated and calibrated signals at 8 receivers at a distance of 4.7 m from the transmitter belt.

Figure 3 .

 3 [START_REF] Blouin | Improved resolution and signal-to-noise ratio in laser-ultrasonics by SAFT processing[END_REF] shows the emitter (orange) and receiver (blue) pair placed at diametrically opposite slots in the transducer ring used in the pitch catch setup. The series of images also shows the addition of inactive transducers (grey) in progression.

Figure 3 . 6 :

 36 Figure 3.6: Setup to study circumferential attenuation due to the addition of transducers. The wave propagates from the emitter (orange) to the receiver (blue). The inactive transducers (grey) are added sequentially to measure their additional impedance at the surface.

Figure 3 . 7 :

 37 Figure 3.7: Signals measured at the receiver showing increased attenuation with the addition of transducers circumferentially.

Figure 3 . 8 .

 38 Figure 3.8. The multiple reflections received by this transducer is shown in Figure 3.9.

Figure 3 . 8 :

 38 Figure 3.8: Measurement of velocity with multiple through thickness reflections.

Figure 3 . 9 :

 39 Figure 3.9: Signal of multiple reflections of the longitudinal bulk wave over the thickness.

Figure 3 . 10 :

 310 Figure 3.10: Frequency spectrum of reflected signal as obtained in Figure 3.9.

Figure 3 . 11 :

 311 Figure 3.11: Spectrum of frequency spectrum (as obtained in Figure 3.10) obtained by taking a Fourier transform of the absolute value of the spectrum giving the time between successive reflection.

Figure 3 .

 3 Figure3.12 compares the normalised axial displacement field created by transducers with a fixed axial length of 15.6 mm and varying widths. An important point to note is that the amplitudes of the wavefronts are directly related to the dimension of the transducer i.e. doubling any one of the dimensions leads to a doubling of the amplitude. It is clear that the circumferential width of the transducer has a negligible effect on the wave fronts. The difference between the normalised wave-fronts can only be found on closer observation of the second and third packets. The differences in magnitudes of these packets are in the order of 0.4 %.

Figure 3 . 12 :

 312 Figure 3.12: Normalised displacement field signal for Gaussian sources of different widths measured at 50 cm. Figure also shows the packets that represent various modes.

  Figures 3.14a and 3.14b show the amplitude distribution of stress or displacement loads of the N = 0 and N = 4 modes around the circumference of the cylinder.

  (a) Amplitude distribution for N = 0. (b) Amplitude distribution for N = 4.

Figure 3 . 14 :

 314 Figure 3.14: Amplitude distribution of two modes around the circumference of the pipe. These may represent stress or displacement amplitudes.

  Figures 3.12 and 3.13 also inadvertently end up validating the presence and arrival of the F (8, 1) mode.

Figure 3 . 15 :

 315 Figure 3.15: Frequency-axial-wavenumber dispersion curve showing the fundamental modes from circumferential wavenumber 8 to 15 in relation to 100 kHz central frequency.

Figure 3 . 17 :

 317 Figure 3.17: Comparison of the computed (orange) and experimentally measured (blue) group velocities of the axisymmetric (N = 0) mode [L(0,2)] measured at 50 and 100 cm.

Figure 3 . 18 :

 318 Figure 3.18: Comparison of the computed (orange) and experimentally measured (blue) group velocities of a non axisymmetric (N = 1) mode [F(1,3)] measured at 50 and 100 cm.

Figure 3 . 19 :

 319 Figure 3.19: Comparison of the computed (orange) and experimentally measured (blue) group velocities of a non axisymmetric (N = 4) mode [F(4,2)] measured at 50 and 100 cm.

Figure 3 . 20 :

 320 Figure 3.20: Spectrum of the axial displacement component of the computed field obtained by exciting the L(0,2) mode.

Figure 3 . 21 :

 321 Figure 3.21: Spectrum of the axial displacement component of the computed field obtained by exciting the F(1,3) mode.

Figure 3 . 22 :

 322 Figure 3.22: Spectrum of the axial displacement component of the computed field obtained by exciting the F(4,2) mode.

Figure 3 . 23 :

 323 Figure 3.23: Comparison of theoretical group velocities (solid lines) of the three studied modes [L(0,2), F(1,3) and F(4,2)] with various trials of experimentally obtained group velocities showing good agreement.

Figure 4 . 1 :

 41 Figure 4.1: Snapshots over time showing the forward propagated adjoint field (top), the backward propagated direct field (middle) and the topological image of the defect forming by the correlation of the two aforementioned fields (bottom).

( a )

 a Topological image computed without using the spatial envelope. (b) Topological image computed with the spatial envelope.

Figure 4 . 2 :

 42 Figure 4.2: Comparison of the topological images obtained by the cross correlation of the adjoint and time reversed direct fields to identify the defect (orange contour). The image on the left is the actual image and that on the right is computed by taking the spatial envelope of the first image.

Chapter 2

 2 have been imaged. The adjoint field is created by measuring the field diffracted by the defect at 24 equiangular positions around the circumference of the cylinder at the same axial position (z = 0) as the continuous source belt. These signals are then temporally inverted and used as the source for the adjoint field. Certain examples have been specifically chosen so as to provide an insight to the functioning of the defect model.

( a )

 a Topological image of the defect approximated by the 1st order derivative.(b) Topological image of the defect approximated by the 4th order derivative.

Figure 4 . 3 :

 43 Figure 4.3: Topological images of the 2 mm defect modelled using different orders of the base gaussian function.

Figure 4 .

 4 Figure 4.3a may seem surprising at first due to the two patches on either side of the defect. On further analysis, these patches are understood to be the one of the principle secondary sources that approximates the defect i.e. the secondary sources as modelled in Chapter 2. This secondary source acts to minimise the incident field within the contour of the defect. In this case, it represents the combination of the first order derivative in the axial direction and the zeroth order in the circumferential direction. The peaks of this secondary source lie well beyond the contour of the defect.Using the 4th order derivative diverges the solution and this is observed in its topological image as well (Figure4.3b) where the image of the defect is not coherent.

  4th Order 4 mm (d) 5th Order 4 mm

Figure 4 . 4 :

 44 Figure 4.4: Topological Images of 4mm defect modelled using 1st, 3rd, 4th and 5th orders.

Figure 4 .

 4 4c suggests that the dominant secondary source consists of a 4th order derivative in the circumferential direction with the 0th derivative along z. This makes the secondary source slimmer in the θ direction. Taking the 5th order derivative (Figure4.4d), the image clearly shows a secondary source which may be a combination of two different derivatives. One of these sources is not in the vicinity of the defect contour. This leads to a grey area in the definition of the new defect model and is discussed further in the conclusions of this section.

Figure 4 . 5 :

 45 Figure 4.5: Topological Images of 6mm defect modelled using 2nd and 4th orders.

  (a) TI formed with 2 transducers showing a lot of artifacts. (b) TI formed with 4 transducers showing some artifacts around the defect. (c) TI formed with 6 transducers showing clear defect with minimal artifacts. (d) TI formed with 8 transducers showing clear defect image. (e) TI formed with 12 transducers showing clear defect with additional artifacts near transducers.

Figure 4 . 6 :

 46 Figure 4.6: Topological images formed using a number of transducers for adjoint field computations showing that increasing the number of transducers improves accuracy and sharpness of the image.

Figure 4 .

 4 Figure 4.6 shows the topological images of the defect using a varying number of transducers in pulse echo mode. Using two transducers (Figure 4.6a) gives an idea of the position of the defect. In an unbounded medium, three non-dispersive wave packets should

  each of the cases, 12 transducers have been used to create the adjoint field and the defect has been imaged at 100, 150, 200 and 250 kHz.

  (a) TI formed at 100 kHz (b) TI formed at 150 kHz (c) TI formed at 200 kHz (d) TI formed at 250 kHz

Figure 4 . 7 :

 47 Figure 4.7: Series of images showing the topological image of the same defect imaged at four different frequencies ; 100, 150, 200 and 250 kHz.

of 5 .

 5 6 mm. The incident field is a 5 cycle 200 kHz axisymmetric guided wave generated by a continuous circumferential belt at z=0. The diffracted field is measured at 24 equiangular positions around the circumference which is then used to create the adjoint field. The topological image of this defect has been shown in Image 4.8.

Figure 4 . 8 :

 48 Figure 4.8: Topological image of a defect 5.5 mm wide axially and 100 mm long circumferentially. The orange contour represents the delamination defect.

( a )

 a Contour of the defect created to represent two defects separated axially and circumferentially. (b) Topological image of the two separated defects lying within the same defect contour showing that the method of modelling the defect can be used to simulate multiple defects.

Figure 4 . 9 :

 49 Figure 4.9: Figures showing the single contour to describe axially and circumferentially separated defects (left) and the topological image (right) identification of the two defects.

  2. The configuration consists of a steel pipe with inner and outer diameters of 56 mm and 60 mm, respectively. A transducer belt with 8 individual elements equally spaced is located at z=0. A through thickness crack defect spanning 50°is located at a distance of z=750 mm.The setup along with the defect specifications is given in Figure4.10.

Figure 4 . 11 :

 411 Figure 4.11: Snapshot of the incident axial displacement field at the surface of the pipe at 140 µs. The incident field is generated by 8 equiangular transducers exciting a 15-cycle tone burst signal at 100 kHz.

Figure 4 . 12 :

 412 Figure 4.12: The signals of the diffracted field received at the 8 transducer positions measured in CIVA. These signals, temporally inverted are used as the input to create the adjoint field in the numerical tool.

Figure 4 . 13 :

 413 Figure 4.13: Topological image of the defect modeled in CIVA and localised by the numerical model using only the reflected signals obtained from the CIVA simulation shown in Figure 4.12.

  First, a study is performed to image defects of various sizes with theoretical residuals. Next, the potential of the experimental residual is observed by imaging concealed and extremely small defects. The residual is obtained by taking the difference between signals from a medium with a discontinuity and a healthy medium. The signal without defect contains the information of the healthy medium and is generally numerically computed. The signal with defect additionally contains information of any discontinuity. The concepts of theoretical and experimental residuals have been further detailed in Section 4.3.1. Multiple modes are then used to improve the quality of the TI. Finally, certain limitations of the experimental setup have been discussed to conclude the chapter.

5 mm, 10

 10 mm and 15 mm have been studied. Each of these defects are placed at an axial distance of 81.5 cm from the receiver belt and at an azimuth corresponding to 0°.Experimental residuals (shown for a single transducer in Figure4.14) are obtained by subtracting the response of the damaged medium from that of a healthy medium (before affixing the defects) measured at the same receiver positions. It is clear that as the size of the defect increases, so does the amplitude of the reflected signal. It is also clear that the defect signature is clearly distinguished from the background noise by taking the experimental residual.

Figure 4 . 14 :

 414 Figure 4.14: Experimental residuals of signals diffracted from 5 mm (blue), 10 mm (orange) and 15 mm (green) defects measured at a single receiver obtained by experimental baseline subtraction.

Figure 4 . 15 :

 415 Figure 4.15: Theoretical residuals of signals diffracted from 5 mm (blue), 10 mm (orange) and 15 mm (green) defects measured at a single receiver computed with the numerical simulation.

  (a) Topological image of the 15 mm defect (b) Topological image of the 10 mm defect (c) Topological image of the 5 mm defect

Figure 4 . 16 :

 416 Figure 4.16: Topological images of the 15 mm, 10 mm and 5 mm Metglas ® defects showing the decrease in amplitude of the patch representing the defect, signifying the decrease in sensitivity based on defect dimensions.

( a )

 a Topological image of the pencil tip without experimental baseline subtraction (Theoretical residual). (b) Topological image of the pencil tip using experimental baseline subtraction (Experimental residual).

Figure 4 . 17 :

 417 Figure 4.17: Comparison of the topological images of the pencil tip used to mimic a small defect obtained using the theoretical and experimental residuals.

Figure 4 . 18 :

 418 Figure 4.18: Topological image of a compass tip used as an extremely small diffractor computed using the experimental residual.

Figure 4 .

 4 Figure 4.19 shows the Topological Image of the aforementioned setup using the baseline of a healthy medium devoid of either the defect or the circumferential strip.

Figure 4 . 19 :

 419 Figure 4.19: Topological image of a defect behind a continuous Metglas ® belt without experimental baseline subtraction showing the inability to image the defect.

Figure 4 . 20 :

 420 Figure 4.20: Topological image of a defect behind a continuous Metglas ® belt using experimental baseline subtraction clearly showing defect.

Figure 4 . 21 :

 421 Figure 4.21: Topological image of two defects at an angular position of 0°, separated axially by 100 mm radiated using the L(0, 2) mode.

Figure 4 . 22 :

 422 Figure 4.22: Topological image of two defects at an angular position of 0°, separated axially by 100 mm radiated using the F (1, 3) mode.

Figure 4 . 23 :

 423 Figure 4.23: Combined image of two defects at an angular position of 0°, separated axially by 100 mm by taking the pointwise product of the imaged obtained by radiating the specimen with the L(0, 2) and F (1, 3) modes.

Figure 4 . 24 :

 424 Figure 4.24: Topological image of two defects at angular positions 0°and -90°, separated axially by 100 mm radiated using the L(0, 2) mode.

Figure 4 . 25 :

 425 Figure 4.25: Topological image of two defects at angular positions 0°and -90°, separated axially by 100 mm radiated using the F (1, 3) mode.

Figure 4 . 26 :

 426 Figure 4.26: Combined image of two defects at angular positions 0°and -90°, separated axially by 100 mm by taking the pointwise product of the imaged obtained by radiating the specimen with the L(0, 2) and F (1, 3) modes.

Figure 4 . 27 :

 427 Figure 4.27: Combined image of two defects at angular positions 0°and -90°, separated axially by 100 mm by taking the pointwise mean of the imaged obtained by radiating the specimen with the L(0, 2) and F (1, 3) modes.

Chapter 1

 1 establishes a robust mathematical model for the fast simulation of guided waves in cylindrical structures. This model is based on solving the guided wave propagation equations in layered media in the Fourier and Laplace domains. The layers are combined using the Global Matrix method. The solutions derived are given in the form of partial waves containing stress and displacement vectors, written as combinations of Bessel's functions of the first and second kind. The software implementation of this model is then described. The new semi-analytical mathematical model developed allows for a fast computation of the Green function in tubular structures which could be a viable alternative for conventional finite element simulations.

  These functions are used to create an orthonormal basis onto which the incident field is then projected. The negative sum of components of the field on each of the basis vectors gives the field diffracted by the defect. Defects studies are performed to understand the mechanism by which the model functions to help pick the right parameters. This novel method aids in the modelling and simulation of small defects having dimensions lesser than the wavelength of computation.Chapter 3 discusses the experimental setup used to verify the numerical model and perform tests with defects. The setup, equipment and configuration are described. A technique for normalising transducers is then discussed to ensure the generation of pure modes. Experimental parameters such as material properties and transducer dimensions to be used in the simulation are then obtained. The numerical model is then compared and verified with experimental measurements. Finally some limitations and capabilities of the model are discussed.Chapter 4 introduces the method of Topological Imaging and its implementation in the context of this thesis. After introducing the principle of the method, numerical defects as defined in Chapter 2 are imaged. This section also sheads some light onto the functioning of the defect model. Next, a defect defined on an external simulation tool is imaged using only the signals collected at the receivers. Physical experiments are then done showing the effectiveness of using the theoretical and experimental resides in creating the image. Defects having dimensions as low as 1 /40th of the wavelength of the guided wave have been imaged. This versatile method of imaging numerical and experimental discontinuities of various dimensions may also be used with other numerical techniques regardless of the complexity of the media or propagating wave for fast and accurate defect localisation.The mathematical models, methods and techniques described in this thesis offer a quick and robust way of simulating guided wave propagation in layered isotropic tubular structures as well as the modelling of numerical delamination defects and localising said numerical and experimental defects accurately. The techniques also open up possibilities for interesting future work.The model can be expanded to simulate guided wave propagation in anisotropic structures. This would allow it to be used for wave propagation simulation in composite structures such as airplane fuselages. There is also scope for improving and expanding the defect model. It was developed during the thesis and has not yet been fully explored.

  Creation of orthonormal basis: general case . . . . . . . . . . . . . 2.2 Gaussian secondary source and its derivatives . . . . . . . . . . . . . . . . 2.2.1 Definition of secondary source: Summary . . . . . . . . . . . . . . 2.3 Defect studies based on defect size . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Case 1: radius of defect 6 mm . . . . . . . . . . . . . . . . . . . . . 2.3.2 Case 2: radius of defect 4 mm . . . . . . . . . . . . . . . . . . . . . 2.3.3 Case 3: radius of defect 2 mm . . . . . . . . . . . . . . . . . . . . . 2.4 Observations and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . viii 3 Experimental Identification of Parameters for the Numerical Model 3.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Details of experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Normalisation of transducers . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Attenuative influence of multiple transducers on incident field . . . 3.2 Dimensioning of numerical model . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Material characterisation . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Modelling of source: transducer sizing . . . . . . . . . . . . . . . . . 3.3 Generation and comparison of modes : experiment vs simulation . . . . . . 3.4 Capabilities and limitations of the numerical model . . . . . . . . . . . . . Principle of the method . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	4 Topological Imaging
	4.1

2 Diffraction of a Field by a Defect 2.1 Equations to solve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Inner product associated with the defect shape . . . . . . . . . . . . 2.1.2 4.2 Numerical experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Imaging defects of multiple sizes . . . . . . . . . . . . . . . . . . . . 4.2.2 Effect of number of transducers on image . . . . . . . . . . . . . . . 4.2.3 Effect of frequency on image . . . . . . . . . . . . . . . . . . . . . . 4.2.4 Alternate defect forms . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1

 1

	.1.
		Time domain t notation Laplace domain s notation
	Physical space r, θ, z	(r, θ, z, t)	u	(r, θ, z, s)	U
	Azimuthal and axial wavenumbers n, k	(r, n, k, t)	ũ	(r, n, k, s)	Ũ
	Radial position r
	Azimuthal and axial wavenumbers n, k	(κ, n, k, t)	û	(κ, n, k, s)	Û
	Radial wavenumber κ

Table 1 .

 1

1: Notations.

Table 3 .

 3 .1.

	Mass density ρ	Longitudinal velocity C L Transverse velocity C T
	2.703 mg/mm 3	6.254 mm/µs	3.052 mm/µs

1: Measured mass density and sound velocities of aluminum.

Table 3 .

 3 .2 below. 2: Dimensions tested for sizing the Gaussian distribution representing the transducer.

	Percentage Length(mm) Width(mm)
	80 %	10.4	2.4
	100 %	13.0	3.0
	120 %	15.6	3.6

Table 4 .

 4 1: Cylinder dimensions and material properties.This setup generates purely axi-symmetric modes and is done to avoid the effects of highly dispersive (refer Image 3.15), slowly propagating or unwanted modes (which usually cause artifacts), ensuring that the study only deals with the factors being discussed.

	Length	1000 mm
	Inner radius	28 mm
	Outer radius	30 mm
	Density ρ	2.703 mg/mm 3
	Longitudinal bulk velocity C L 6.254 mm/µs
	Transverse bulk velocity C T	3.052 mm/µs

Acknowledgements

Acknowledgements vi

This thesis was carried out at the Institute of Mechanics and Engineering [l'Institut de mécanique et d'ingénierie] (I2M) in the University of Bordeaux in collaboration with and funded by CEA (Commissariat à l'énergie atomique et

Next, a reference amplitude is taken as the average of the entire M matrix and is represented as μ. This would give:

Summing Equation (3.1) for a fixed receiver position j, i µ ij = nμ + nρ j ; ∀j gives the coefficients for the receivers as:

Similarly, fixing i, the coefficients for the transmitters are given as:

Effect of axial length on amplitude and wavefront

The variation of the axial length has a substantial impact on the wave-fronts of the second and third packets. The fields are computed keeping the width as a constant of 3.6 mm as decided in the earlier section and varying the lengths.

On comparing the normalised signals as shown in Figure 3.13, a smaller length of the transducer has the effect of an increased amplitude of the second and third wave packets. This is because a smaller transducer has a larger angle of directivity in the circumferential direction. Intuitively, the mode L(0, 1) being axially symmetric should scale similar to the L(0, 2) mode. However, this changing behaviour can be attributed to the fact that the tested lengths of the transducer are comparable to the wavelength of the L(0, 1) mode and as the length approaches the wavelength, there is a decrease in amplitude. Even though these modes have a small part to play in future studies, this particular behaviour is useful in sizing the length of the transducer. The incident field is generated by exciting all eight transducers in phase with a 15-cycle 100 kHz toneburst signal. A snapshot of this field is shown in Figure 4.11. The advantage of using a large number of cycles is that it limits the frequency bandwidth of the spectrum, reducing the chances of other modes being generated. The drawback is that the signal is temporally large which inturn leads to the field stretching to around 400mm spatially.

Multi-defect imaging with multiple modes

This section aims to demonstrate some techniques to improve clarity of images using multiple modes, while also presenting certain limitations. The idea of using multiple modes arises from the fact that each of the modes have unique stress and displacement distributions not just across the thickness but also about the circumference. For example, the L(0, 1) (similar to the anti-symmetric mode in plates) has a radial displacement component that is constant over the thickness whereas the L(0, 2) mode (similar to the symmetric mode in plates) has a radial component at the upper and lower surfaces which are out of phase by π radians. This means that each of the modes illuminate different sections of the test specimen with different intensities, opening avenues to improving the accuracy of the images.

In this study, the L(0, 2) and the F (1, 3) modes have been utilised to improve the clarity of the images of the defects at two critical positions. These modes belong to the same family with the only exception being that the F (1, 3) mode has an amplitude that varies as a sinusoid around the circumference. The position of the first defect is fixed at θ=0°and z=81.5 cm. The second defect is positioned axially at z=91.5 cm. This defect is located at an azimuth of 0°and -90°for the first and second cases respectively. The pipe is illuminated by each of the two modes before affixing the defects to obtain the experimental baseline signal. It should be noted that the direct field for the Topological Image would be different for each of the two modes.