Transient Propagation and Scattering of Quasi-Rayleigh Waves in Plates: Quantitative comparison between Pulsed TV-Holography Measurements and FC(Gram) elastodynamic simulations

Transient Propagation and Scattering of Quasi-Rayleigh Waves in Plates: Quantitative comparison between Pulsed TV-Holography Measurements and FC(Gram) elastodynamic simulations

Faisal Amlani Department of Aerospace and Mechanical Engineering, University of Southern California, Los Angeles, CA 90007 (USA) Oscar P. Bruno Applied and Computational Mathematics, California Institute of Technology 217-50, Pasadena, California 91125 (USA) José Carlos López-Vázquez jclopez@uvigo.es Cristina Trillo Ángel F. Doval José L. Fernández Pablo Rodríguez-Gómez Departamento de Física Aplicada, Universidade de Vigo, Escola de Enxeñería Industrial, Rua Maxwell, s/n, Campus Universitario, E36310 Vigo, Spain
Abstract

We study the scattering of transient, high-frequency, narrow-band quasi-Rayleigh elastic waves by through-thickness holes in aluminum plates, in the framework of ultrasonic nondestructive testing (NDT) based on full-field optical detection. Sequences of the instantaneous two-dimensional (2-D) out-of-plane displacement scattering maps are measured with a self-developed Pulsed TV Holography (PTVH) system. The corresponding simulated sequences are obtained by means of a FC(Gram) elastodynamic solver introduced recently, which implements a full three-dimensional (3D) vector formulation of the direct linear-elasticity scattering problem. A detailed quantitative comparison between these experimental and numerical sequences, which is presented here for the first time, shows very good agreement both in the amplitude and the phase of the acoustic field in the forward, lateral and backscattering areas. It is thus suggested that the combination of the PTVH system and the FC(Gram) elastodynamic solver provides an effective ultrasonic inspection tool for plate-like structures, with a significant potential for ultrasonic NDT applications.

keywords:
Ultrasonic nondestructive testing, Elastic waves in plates, Pulsed TV holography, Spectral simulation, Fourier continuation, High order accuracy

Abbreviations and acronyms111Charged Coupled Device (CCD), Courant-Friedrichs-Levy (CFL), Fast Fourier Transform (FFT), Finite Differences (FD), Finite Element Method (FEM), Fourier Continuation (FC), Initial Boundary Value Problem (IBVP), Lead Zirconate Titanate (PZT), Ordinary Differential Equation (ODE), Partial differential Equation (PDE), Pulsed TV Holography (PTVH), Region Of Interest (ROI), Singular Value Decomposition (SVD), Ultrasonic Nondestructive Testing (NDT).
journal: Mechanical Systems and Signal Processing

1 Introduction

We study the scattering of transient, high-frequency, narrow-band quasi-Rayleigh elastic waves by through-thickness holes in aluminum plates, in the framework of ultrasonic nondestructive testing (NDT) based on full-field optical detection. Our work combines an innovative pulsed TV holography (PTVH) experimental system fdezijo07 () and a state-of-the-art three-dimensional (3D) FC(Gram) numerical elastodynamic solver amlani2016fc (). For the first time, a detailed quantitative comparison between experimental and numerical sequences of the instantaneous out-of-plane displacements at the fully two-dimensional (2D) top surface of the plates is presented, showing very good agreement both for the amplitude and the phase of the acoustic field in the forward, lateral and backscattering areas. It is thus suggested that the combination of the PTVH system and the FC(Gram) elastodynamic solver provides an effective ultrasonic inspection tool for plate-like structures.

As it is well-known, ultrasonic techniques are a mature and powerful tool for NDT and evaluation in industry HNDT_Bricks () with well-known benefits (deep-penetration capability with high degree of interaction with flaws or inhomogeneities, wide temporal bandwidth, high information content,…) and a wide range of different applications, particularly, in the case of ultrasonic guided waves NDT Rose_hot (); ews302 (). Its basic principle can be described in three basic steps: first, an adequate ultrasonic field is insonified in the part to be inspected; second, interaction with defects in the material produce scattering of the insonified field; third, some physical quantity associated with the acoustic field is measured providing information about the associated scattering phenomena.

This general scheme is implemented in classical pointwise techniques employing pulse-echo or pitch-catch configurations HNDT_Bricks () and using excitation and detection at one or two points with contact transducers (typically PZT), which renders outputs signals with high resolution in time and frequency, but normally with low spatial information content. From this well-established classical schemes, ultrasonic NDT technology has progressed in the last decades in two main complementary directions: on the one hand, developing non-contact inspection approaches noncontact () (f.i. air-coupled devices aircoupled (), electromagnetic acoustic transducers EMATS (); inplane () or pointwise optical generation and probing of ultrasound Scruby (); int35 ()) and, on the other hand, obtaining more spatial information content about the ultrasonic field in contact techniques (f.i. pointwise techniques with phased array transducers to control the directivity of excited or detected waves phase () or combined with some kind of scanning in detection Knuuttila_IJO () or excitation tom35 () or techniques with transducers spatially distributed or scanned over the part to be inspected combined with tomographic approaches ews57b ()).

In particular, optical probing of ultrasound with full-field techniques (f.i. interferometry Nakano_IJO (), holographic interferometry Blackshire_IJO (); ews59 (), speckle photography kum01 (), shearography ews68 () or TV-holography trilloao03 ()), are intrinsically non-contact and specially well-suited for obtaining spatial information by registering images with high lateral spatial resolution and irradiance distributions correlated with the acoustic field at a given time interval, and allow the inspection or remote or inaccessible areas maintaining the previously mentioned classical benefits of NDT with ultrasound. Within this framework, those of us at the University of Vigo have demonstrated that ultrasonic non-destructive inspection can be performed in plates using a self-developed pulsed TV-holography (PTVH) system fdezijo07 (); trilloao03 (); trillomst03 (); cernadasmssp06 () that records the two-dimensional (2D) spatial distribution of the acoustic field corresponding to the instantaneous out-of-plane displacements over the surface of the plate. The measured ultrasonic 2D field maps propagation and scattering patterns from which characterization of defects (position, dimensions, orientation, etc.) could be extracted by solving the corresponding inverse scattering problem. As it is characteristic of full-field optical techniques, our PTVH output has high spatial resolution with a large number of spatial samples (about ), even though the number of temporal samples is low (typically lower than 128). Also, in our case, the field of view for the detected 2D scattering pattern includes tens of ultrasonic wavelengths both in the near-field and far-field zones, which means that the corresponding scattering phenomena are in the high frequency range scs213 ().

All ultrasound-based NDT approaches require an additional fourth assessment step, namely, processing and analysis of the output signal—leading to subsequent detection, location and characterization of any existing flaws ews287 (). This assessment process implies, implicitly or explicitly, some kind of modelling of wave propagation and scattering review_Bond (), which is usually developed in the framework of the rigorous three-dimensional (3D) vector linear elasticity theory, in conjunction with a great variety of analytical or numerical schemes (e.g. normal mode expansion method ews91 (), finite difference method ews291 (), finite element method (FEM) ews292 (), boundary element method Cho_Rose_ews71 (), discrete point source method ews79 (), local interaction simulation approach ews290 (), various types of combinations or hybrid methods Liu_ews72 (); ews254 (), etc.). In order to avoid the high computational cost typically associated to standard 3D simulation approaches in the high frequency regime, approximate theories (i.e. plate theories ews197 (); Wang_Chang_2005 () or other even more highly simplified appoximations mast-gordon ()), which are valid only for a limited number of situations, but that occur often in practice, have been utilized frequently ews199 (); OEng2010 (); OEng2013 ().

The computational solvers employed in this contribution are based on a novel “spectral” approach: the recently introduced Fourier-Continuation (FC) method presented in amlani2016fc (); bruno2010high (); albin2011spectral (); albin2012fourier (). Like regular spectral solvers boyd2001chebyshev (), the FC methods represent solutions by means of Fourier series, and they therefore make it possible to dramatically reduce the problematic dispersion errors—which accumulate in the finite-difference or finite-element approximation of the solution over domains spanning many wavelengths. But, unlike the classical Fourier-based method, the FC approach is applicable to general domains. And, unlike the Chebyshev-based methods, the FC approach utilizes regular Cartesian grids, and it therefore does not suffer from the severe Courant-Friedrichs-Levy (CFL) constraints that stem from the point-clustering inherent in Chebyshev-based explicit solvers.

As detailed in the aforementioned references, the ability of the FC method to provide high-order, low-dispersion, spectral accuracy, for general domains, stems from two main strategies, namely, 1) decomposition of the computational domain into a number of overlapping patches (following henshaw2008parallel ()); and 2) Fourier expansion along each coordinate axis, in each one of the patches, by creation of periodic functions (the continuation functions) defined in adequately chosen domains beyond the physical patch. In view of its low dispersion, the elastic FC solver has provided improvements by factors of hundreds, or even thousands and beyond, over the computing times required by other methods, for a given accuracy, for the problem of propagation of high-frequency elastic waves amlani2016fc ().

As mentioned above, this contribution compares experimental values measured by PTVH with results obtained from the aforementioned FC-based elastic-wave solver, for the problem of scattering of transient quasi-Rayleigh waves by through-thickness holes. The paper is organized as follows. Section 2 establishes the nomenclature and the theoretical formulation of the direct scattering problem within the framework of the linear elasticity. Section 3 then describes details concerning the experimental methods (for excitation of quasi-Rayleigh transient wavetrains, measurement of the displacement field by PTVH and the two-step spatio-temporal Fourier transform method for deriving the final complex displacement field) as well as the numerical implementation of the FC elastic-wave solver approach. In both cases only the essential points are included for completeness: the PTVH technique was described in detail in previous contributions (see references 1,  21,  22 and  46), while a detailed description of the FC(Gram) elastic solver was presented in amlani2016fc (). The experimental and simulated maps are presented in section 4, including comparisons over various specific areas in the 2-D image of the acoustic field, and the agreement between theory and experiment is discussed taking into account the pixel-wise matching of the spatial field distribution, the profiles of amplitude and phase, and the corresponding values of the global (root-mean-square) relative errors. To the best of our knowledge this is the first time that such a quantitative comparison is reported for transient waves both in near- and far-field zones, apart from our corresponding recent studies in which this work is based amlani2016fc (); Tesis_PRG (). A full 3D sequence generated by the FC(Gram) solver, displaying the temporal and spatial variation of all three components of the displacement vector both in the interior and the lateral boundaries of the plate, is also included in this section—providing useful insights on rarely observed mechanisms concerning elastic scattering. Section 5, finally, present our closing remarks.

2 Theoretical framework

We assume that for modelling wave propagation and scattering in our case it is adequate to describe plates as linear isotropic solids in the framework of linear elasticity theory within small displacement regime and linear stress-strain relationship rose00 (); Achenbach_R (); Graff (). In this context, first in subsection 2.1 we introduce the notation and establish the basic equations that define the initial boundary value problem (IBVP) to be solved for the scattering of waves in plates with stress-free boundaries in transient regime. After that, in subsection 2.2 we consider the particular properties that can be used in the case of the scattering of quasi-Rayleigh wavetrains.

2.1 Formulation of the direct scattering problem

Using the Cartesian reference frame of figure 1, with origin at point and orthonormal base vectors , and , the position of a material point of the plate in the reference configuration is characterized by its Cartesian coordinates which are the components of its position vector (for the shake brevity we simply write ). When a wave perturbation propagates in the plate, the associated movement of with respect to its reference position at a given time instant is specified by the displacement field

(1)

with components that can always be written as

(2)

where and are the acoustic amplitude and the acoustic phase, respectively ( is referred as the total acoustic phase). Despite the fact expression (2) is formally valid in general Born_Wolf (), it is particularly useful in two important particular cases: on the one hand, for perturbations with harmonic time-dependence at a given temporal frequency with displacement field and harmonic components , because in this case the associated acoustic amplitudes and acoustic phases are truly time-independent fields; on the other hand, for transient narrow-band perturbations with a central temporal frequency , as those employed in our experiments, because the associated acoustic amplitudes and acoustic phases are nearly constant within any time interval of the order of the central period of the wavetrain in such a way that they can be understood in close analogy with the harmonic case, except perhaps in areas that are in the leading and trailing edges of the wavetrain (see figure 2).

(a) (b)
Figure 1: Scheme of a plate of thickness with a cylindrical hole of diameter and residual depth : (a) section of the plate showing the detail of the instantaneous surface displacement associated to a harmonic guided wave with in-plane component and out-of-plane component that propagates along the axis, (b) top side of the plate, showing the position of point , the intersection of the virtual line-source parallel to axis associated to the quasi-cylindrical wavetrain of the incident field.

Considering that the wave perturbation can be described using the small strain regime, we have an associated strain tensor with components

(3)
(a) (b)
Figure 2: (a) Typical temporal record of the out-of-plane component of the incident field at a given point over the top surface of the plate for a narrow-band quasi-Rayleigh wavetrain with central excitation frequency MHz and bandwidth (), (b) amplitude of the temporal spectrum of the profile in (a). In both cases we show the normalized profiles as a function of in s and in MHz, respectively.

Assuming that the plate is a linear homogeneous isotropic solid with volume mass density and a linear elastic response characterized by the Lamé constants and , the linear stress-strain relationship can be described writing the components of the stress tensor as

(4)

where the Kronecker delta. In this context, the dynamic behavior in the interior volume of the plate is described by the Lamé-Navier equations

(5)

which can be derived employing equations (3) and (4) in the balance of momentum equation Pao_rand () and assuming that the interior of the plate is a region free of body forces. The components of the displacement that are solutions of equation (5) can be discomposed in longitudinal (L) and transversal (T) wave components that propagate, respectively, with phase velocities

(6)

Using (6) we get and in such a way that we can express the constants of the Lamé-Navier equations exclusively in terms of the phase velocities.

Any solution of the Lamé-Navier equation in plates with stress-free boundaries (as it is our case) has to verify over its top side , its bottom side and the hole border the stress-free boundary condition

(7)

being and the components of the normal unit vector at point pointing outwards the plate interior . Equations (5-7), with additional initial conditions for displacement and velocity in the plate interior volume , are the essential elements that define the initial boundary value problem (IBVP) to be solved for the scattering of waves in plates with stress-free boundaries in the transient case.

2.2 Scattering of quasi-Rayleigh waves by through-thickness holes

In the following we circumscribe the analysis to the particular case of propagation and scattering of narrow-band quasi-Rayleigh wavetrains by through-thickness holes, with residual depth (figure 1) considering how to obtain the incident field inside the plate from the value of the out-of-plane component of the incident field at the top surface of the plate. More specifically, we use the complex representation of the displacement field (see Appendix A) to present the strategy to reconstruct the complex Cartesian components of the displacement of the incident field at any point inside the plate from the value of the complex out-of-plane component of the displacement at the corresponding point a the top surface of the plate, which is a quantity that can be measured by our PTVH system (see figure 2.a and subsection 3.1).

Figure 3: Frequency spectrum of Lamb propagating modes for an aluminum plate with stress-free boundaries and Poisson’s ratio . and are the normalized frequency and wavenumber respectively being the longitudinal wave velocity. The location of normalized central excitation frequency , normalized wavenumber and normalized value of the temporal bandwidth of the pulse in our case (see section 4) are schematically represented. The branches that cross the horizontal line at correspond to the propagating Lamb modes that could be excited at the corresponding normalized temporal frequency. qR identifies the region of the spectrum that corresponds to quasi-Rayleigh waves.

Quasi-Rayleigh waves result from particular combinations of Lamb modes in the so-called quasi-Rayleigh region of the Lamb spectrum where the associated branches of the S0 and A0 modes nearly overlap (see qR region with in figure 3). A quasi-Rayleigh harmonic guided wave is just the superposition of a symmetric S0 Lamb mode and an anti-symmetric A0 Lamb mode with the same temporal frequency within the qR range and with the same value of their amplitudes at the top surface of the plate. In this conditions, the harmonic A0 and S0 modes have practically the same propagating wavenumber and a common phase velocity (characterized by the common slope of the branches in the qR region) which is the Rayleigh phase velocity of the plate material. The Rayleigh equation

(8)

where and , shows that can be obtained as a function of the longitudinal and transversal phase velocities and , i.e. that only two of the three phase velocities , and are independent. If only harmonic quasi-Rayleigh waves of this type are efficiently excited by the ultrasonic source, the complex components of the complex displacement of the quasi-Rayleigh incident wavetrain (figure 2.a) can be written as the linear superposition

(9)

where each quasi-Rayleigh harmonic complex component of frequency can be discomposed as

(10)

being and , respectively, the corresponding harmonic complex components of the S0 and A0 modes with temporal frequency and a scalar function that specifies the weight of each harmonic quasi-Rayleigh displacement of frequency within the qR range spectrum of the quasi-Rayleigh pulse. We will assume that is peaked around a temporal frequency and takes non zero values over a bandwidth interval (figure 2.b), in such a way that each harmonic component of the pulse has a wavenumber that is very close to the central wavenumber of the train within the quasi-Rayleigh region qR of the Lamb spectrum (in our experiments we employ an excitation with a normalized central temporal frequency and central wavenumber of the order of 3 and 6 respectively -see below subsection 4 and figure 3). Provided that in the scattering problem the hole is located in the far-field region of the source it is justified, on the one hand, to assume that the phase and the amplitude of the incident field can be described as a quasi-cylindrical inhomogeneous wave that diverges from a line-source, parallel to the axis (which intersect the top surface of the plate at point with position vector , see figure 1.b) and, on the other hand, to neglect the contributions of non-propagating modes in such a way that expressions (9-10) give an adequate description of the components of the displacement for the a transient narrow-band quasi-Rayleigh incident field. In addition, as the normalized transversal distribution of the complex components of the displacement field of each harmonic quasi-Rayleigh contribution is a function only of its wavenumber and the thickness of the plate we can also assume, in a first approximation, that all the harmonic quasi-Rayleigh components of the pulse share the same normalized transversal distribution for their displacements (see Appendix B). In these conditions, the value of the complex components of the complex displacement of the quasi-Rayleigh wavetrain inside the volume of the plate can be written as

(11a)
(11b)

in terms of the complex out-of-plane component , being and the normalized transversal distributions of the complex amplitude of the in-plane and out-of-plane components of the S0 mode with frequency and wavenumber and and the corresponding transversal distributions for the A0 mode. Extracting the real part of the result of expression (11) we finally obtain the Cartesian components of the displacement of the incident field at any point inside the plate.

3 Materials and methods

In the following, we first describe the experimental set-up employed to generate and to detect quasi-Rayleigh wavetrains in plates (subsection 3.1), then we outline the foundations of the new FC(Gram) numerical method (subsection 3.2) and finally we present the discrete formulation of the direct scattering problem (subsection 3.3) and the procedure employed for obtaining an analytic approximation to the incident field (subsection 3.4).

3.1 Experimental methods

Several through-thickness holes, with nominal diameters in the range to , have been adequately prepared in aluminum plates with nominal thickness . In all of the cases the plates have been supported so that the constraints at their surface are minimized and plasticine has been used as acoustic absorber at the edges of the plates to minimize reflections of the incident and scattered waves that could disturb the measured acoustic fields inside the region of interest (ROI) (figure 4). The average thickness of the plates and the average diameter of the holes have been determined experimentally with a relative error lower than by averaging several direct measurements obtained with a caliper with resolution of 0.05 mm and a metalographic microscope with a translation stage with resolution of 0.01 mm. The longitudinal wave velocity in the plates has been measured by means of the classical pulse-echo method resulting .

Figure 4: Reference image for a through-thickness hole () of diameter showing the region of interest (ROI) employed in the images and profiles of numerical and experimental ultrasonic fields over the plane and the subzones of the ROI and used for quantitative assessment of the relative global error between theory and experiment (see section 4).
Figure 5: Experimental setup

The experimental system used to generate and to detect the elastic waves in this plates is schematically represented in figure 5. Quasi-Rayleigh transient wavetrains have been generated by means of the classical wedge method using a short tone-burst with central frequency and an integer number of cycles (between 5 and 10). The instantaneous out-of-plane acoustic field at the top surface of the plate associated to the wavetrain during propagation and scattering has been measured with our self-developed double-pulsed TV holography system trilloao03 (). As is common in TV holography techniques Doval2000 () we employ an image-hologram configuration, sensitive to the out-of-plane component of the displacement of the surface points, with the image sensor of a video camera as recording medium. There is not optical reconstruction of the recorded holograms but instead their intensity distribution is electronically processed to render the optical phase-difference map, which depends on the displacements of the specimen surface. A twin-cavity pulsed, injection seeded and frequency doubled Nd:YAG laser (Spectron SL404T), the core of the TV holography system, emits two laser pulses, with a duration of and a adjustable temporal inter-pulse delay , that are employed to record two correlograms in separate frames of a CCD camera (PCO Sensicam Double-Shutter). Each correlogram corresponds to the interference of the reference beam and the object beam scattered back by the plate surface. A processing procedure based on the spatial Fourier transform method has been applied to the correlograms trillomst03 (), which renders the so-called optical phase-change map

(12)

proportional to the instantaneous out-of-plane acoustic displacement field being the wavelength of the laser and the instant of emission of the first laser pulse. The usual scaling factor () between optical phase-change and displacement is doubled in equation (12) by selecting ( of the period associated to central excitation frequency , that is the minimum number of odd half-periods for which the camera can record the two correlograms in different frames).

With this scheme we have performed measurements in successive instants of the wavetrain propagation (, , , , , with ) to obtain a complete sequence of optical phase-change maps (we have selected with , i.e., a quarter of the period associated to the central excitation frequency of the wavetrain). Even though this sequence of maps represents by itself a useful means to assess the interaction of the quasi-Rayleigh wavetrain and the defect for every time instant, it has been processed as a whole employing the spatio-temporal 3D Fourier transform method trillonimes06 () to obtain a sequence of experimental complex displacement fields , from which the acoustic amplitude and the total acoustic phase can be retrieved independently and a reconstructed out-of-plane displacement is obtained with an improved signal-to-noise ratio. The sequence of complex displacement fields is the raw experimental data set for the comparison with numerical simulations and it is used for obtaining an analytic approximation to the incident complex field (see subsection 3.3), which is employed to define a non-homogeneous Dirichlet boundary condition in the discrete formulation of the direct scattering problem that is introduced in the following.

3.2 Numerical Method

In view of their simplicity and versatility, classical differentiation methods such as those based on use of finite differences (FD) or finite elements are often used as part of partial differential equation (PDE) discretization schemes. When applied to differential equations governing wave-like phenomena, however, approaches of this type do give rise to high numerical dispersion: phase errors accumulate over subsequent periods of a wave field, and thus FD and FEM solvers require increasing numbers of discretization points per wavelength to resolve the solution within a given accuracy tolerance as the size of the problem grows. As is well known, methods based on use of Fourier series do not suffer from this problem: a fixed number of discretization points per wavelength suffices for such methods to produce wave-like solutions with a fixed accuracy in arbitrarily large domains (i.e. domains that can contain an arbitrarily large number of wave periods).

Unfortunately, however, classical Fourier methods are only applicable to periodic functions—since Fourier expansion of a non-periodic functions gives rise to the well-known “ringing effect” (the Gibb’s phenomenon gocken ()): the jump that results in the periodic extension of an inherently non-periodic function leads to an “overshoot” in the Fourier representation at points of discontinuity of the periodic continuation of the given function. This inaccuracy does not subside as more terms are added and, in fact, it results in extremely slow pointwise convergence throughout the interior of the computational domain—in addition to the unmitigated error at domain boundaries.

A goal to extend the applicability of Fourier methods (together with its inherent excellent qualities, most notably dispersionless-ness and high-order accuracy) to general non-periodic configurations has lead to the development of the FC methods that form the basis of the numerical solvers utilized in this work. The Fourier continuation (FC) method produces a rapidly-convergent interpolating Fourier series representation of a given function on a region larger than the given physical domain. This is accomplished by relying on a “periodic extension” of the given function, that closely approximates the given function values in the original domain, but which is periodic (albeit in a slightly enlarged domain) and thus suitable for Fourier approximation. The procedure is illustrated in figure 6. In detail, given a function which is defined, without loss of generality, on the unit interval as

the FC method produces a periodic function defined on an extended interval,

which closely approximates on the original interval

Figure 6: Fourier continuation of the non-periodic function given by Triangles/circles and squares represent matching points and continuation points, respectively.

The proposed fully-discrete Fourier continuation algorithms proceed as follows: letting be the number of discretization points over the unit interval (yielding a uniform grid , ) together with point values of the function of interest, the Fourier continuation method produces a -periodic trigonometric polynomial of the form

(13)

that matches the given discrete values of : , and which, for smooth functions , rapidly approaches as the grows. Derivatives of the function can then be easily and accurately be computed through term-by-term differentiation, e.g.

(14)

In the simplest treatment brunohan (); boyd (); brunospringer (), the coefficients of (13) are found by invoking the Singular Value Decomposition (SVD) to solve the least squares system

(15)

Although straightforward, repeated applications of the SVD for time-dependent function values on each line of a three-dimensional mesh in a PDE is, in general, unacceptably costly. To circumvent this difficulty, an FFT accelerated method bruno2010high (); lyonbrunoII (); albin2011spectral () can be used. Relying on numbers of and of function values nearest to the left and right endpoints of the interval (e.g. ), the accelerated FC approach proceeds by projecting the known near-endpoint sets of function values onto certain basis of Gram discretely orthogonal polynomials. Each one of these orthogonal polynomials possesses a discrete continuation, and, thus, the near boundary values, which are given by certain linear combinations of the orthogonal-polynomial basis, can easily be continued as well. Once discrete periodic continuation functions have been determined, the corresponding Fourier coefficients can be obtained via an application of the FFT. The overall FC approximation procedure is demonstrated in Figure 6.

Figure 7: Plate-propagation test. Maximum numerical errors over all space and over one full temporal cycle (defined as the time required for any one crest to travel the length of the plate) are presented for a plane wave solution with increasing number of wavelengths for the 3D elastic wave equation with Neumann boundary conditions.
Figure 8: Close up of computational model, composed of six overlapping patches, for a through-hole on an elastic plate.

Using the FC method to evaluate function values and their derivatives, a time-domain PDE solver can be completed simply by incorporating an adequate ordinary differential equation (ODE) solver. Following amlani2016fc () we utilize the fourth-order Adams-Bashforth (AB4) algorithm to effect the timestepping. The errors produced by the FC procedure for a simple rectangular plate are illustrated in figure 7. For comparison errors resulting from use of a finite-difference method of second order are also included. The figure shows that the FC errors remain constant as the number of wavelengths increase, while the FD errors grow—in a classical manifestation of the dispersion error. Note that, additionally, even for an acoustically small plate (one wavelength along the length of the plate) a total of 150 FD discretization points per wavelength are necessary to achieve an error comparable to that resulting from the FC approach on the basis of only 20 FC points. Other comparisons in the figure provide more details in these regards. For a given accuracy tolerance, the reduced FC meshes lead to substantially reduced computing times in comparison with those required by the FD method. Figure 8 illustrates the overlapping-patch geometrical description used in conjunction with the FC approach. Full details concerning accuracy, geometry treatment and time-stepping methods used can be found in amlani2016fc ().

3.3 Computational specification of the experimental configuration

A discrete version of the IVBP problem specified by equations (5-7) has been developed in a prismatic computational domain with interior volume (figure 9) using a self-developed high-precision MATLAB pre-computation for the FC(Gram) basis and a self-developed C/C++ implementation of the full 3D FC(Gram)-based elasticity solver described in subsection 3.2. Lamé-Navier equation (5-6) has been employed in and stress-free boundary condition (7) has been applied over the top and bottom sides of the prismatic domain (which correspond to a portion of the top and bottom surfaces of the plate) and over the cylindrical boundary I1 (which corresponds to the hole border ). On top of that, the Dirichlet boundary conditions

(16a)
(16b)

have been imposed over the lateral exterior boundary of , non homogeneous ones for introducing the incident field at E1 and homogeneous ones in the remaining part E2. The thickness of the computational domain has been set equal to the nominal thickness of the plates (). The height and the width of the top and bottom sides of the computational domain have been adequately selected in order to avoid interference in the ROI region between direct scattering phenomena and waves reflected by the exterior boundary. Null initial conditions

(17a)
(17b)

have been assumed for both displacement and velocity in .

Figure 9: Top side of the prismatic computational domain, defined by two parallel planes at (the top and bottom sides of the plate) and the plane lateral boundaries E1 and E2, with a through-thickness void, defined by the lateral cylindrical boundary I1 (the hole border). The relative position of the ROI and the reference system employed are shown.

3.4 Analytic approximation to the incident field

For specifying the analytic approximation to the Cartesian components of the incident field used in the non-homogeneous Dirichlet boundary condition at E1, where , we use a two-stage procedure: first, we obtain a complex analytic approximation to the experimental measured out-of-plane complex component of the incident field at the top surface of the plate at E1, where . Then, from this complex field we reconstruct the Cartesian components of the real displacement of the incident field at any point over the whole section of the plate in E1 using the properties of quasi-Rayleigh wavetrains stated in subsection 2.2.

complex field amplitude phase
not used
Table 1: Models employed to obtain the analytic approximation to the out-of-plane component of the incident field at the top surface of the plate. The best-fit values of the parameters of the three models are obtained using a least squares procedure over three selected experimental profiles along the , and directions, respectively. The best-fit values and of the parameters and are used to obtain the central Rayleigh wavenumber and the central temporal frequency , respectively.

To obtain in the first stage, on the one hand, we use the part of the temporal record of the out-of-plane displacement at the top surface of the plate in E1 which corresponds to the incident field (see figure 10.a). On the other hand, from the experimental sequence (see f.i. figure 11), we select a particular frame at reference time instant in which the wavetrain is completely inside the ROI but without interaction with the hole. From the spatial distribution of the amplitude and phase in this particular frame we can assume that the phase and the amplitude of the incident field can be described as a quasi-cylindrical inhomogeneous wave that diverges from a line-source, parallel to the axis, which intersect the top surface of the plate at point (see figure 1.b). In this framework, we derive the complex analytic approximation by following an objective matching procedure in 6 steps Tesis_PRG (): 1) the position of the virtual source at the top surface of the plate is obtained as the center of curvature of the wavefronts of the incident field in the selected frame at ; 2) using over this frame we select the complex profile and from the best-fit value of the slope of a linear fit to its acoustic phase (the phase of the field in table 1) we estimate the central wavenumber of the train ; 3) over the temporal record we select a reference point with coordinates located roughly at the center of the part of the temporal record corresponding to the incident field (figure 10.a), obtaining the profiles and which characterize at the distribution of the incident field in the and directions, respectively; 4) it is assumed that the the out-of-plane component of the incident field at the top surface of the plate at E1 can be factorized as

(18)

where the models for and are congruent with the usual description of the phase of a cylindrical wave in the Fresnel approximation (table 1); 5) the best-fit values of the parameters of these models are found by matching to and to using a least squares procedure. In particular, the best-fit value (see table 1) divided by gives an estimate of the central temporal frequency of the wavetrain, that must be close to the nominal value of the selected excitation frequency. From and we obtain the Rayleigh velocity of the plate as that combined with allows us to determine any other relevant parameter that could be used to characterize the wave propagation in the plate (e.g. the phase velocity of transversal waves , the Lamé constants of the plate and , the Poisson’s ratio of the plate and the associated normalized spectrum of Lamb modes, etc.); 6) and are global renormalization factors in amplitude and phase that provide, if necessary, additional degrees of freedom to obtain a better global matching of the analytical approximation considering all the values of the complex amplitude of the incident field within the temporal record .

(a) (b)
Figure 10: Temporal record of the out-of-plane displacement at boundary E1 (): (a) experiment , (b) reconstructed incident field employing the analytic approximation . Point with coordinates and is marked in white. Only the rectangular subzone on the left in (a) has been employed for obtaining the analytic approximation to the non-homogeneous Dirichlet boundary condition in (b). Displacement in units of . in millimeters and in microseconds. The axis numbering corresponds to that of figure 4. Midgray level represents zero.

In the second stage we use expression (11) to reconstruct the components of the complex displacement of the incident field at any point over the whole section of the plate in E1 as

(19a)
(19b)

from which we derive the analytic approximation to the real displacement incident field as

(20)

Auxiliary real envelopes and (with pulse profiles with a smooth transition in amplitude from one to zero along and , respectively) have been employed in expression (20) to guarantee proper matching of the inhomogeneous boundary condition at E1 with the lateral homogeneous boundary condition at E2 and the zero initial values in the computational domain imposed by expression (17). Practical implementation and processing of numerical data based on this approach has been developed employing self-developed routines in Matlab.

4 Results and discussion

We have employed the previously described methodology to analyze the experimental and numerical sequences corresponding to the scattering of quasi-Rayleigh wavetrains by four through-thickness cylindrical holes with nominal diameters , , and . After obtaining the experimental sequence for each hole using the procedure described in subsection 3.1, we obtained the analytical approximation to the incident field and the propagation parameters that characterize the wave propagation in the plates using the procedure described in subsection 3.4. We found that the analytical approximation matches the experimental incident field in the corresponding temporal record in all cases, even though pixel to pixel agreement (figures 10.c and d) presents errors at some points in concordance with previously measured noise levels in our case, that are typically around of the mean signal level, both in amplitude and phase Tesis_PRG (). In addition, the values of the propagation parameters for each hole and their associated sample averages (table 2) have coherent and reasonable values: on the one hand, the measured central temporal frequency matches the nominal excitation frequency and the values of the Rayleigh phase velocity, the transversal phase velocity and the Poisson’s ratio are within the typical values reported for aluminum; on the other hand, the corresponding normalized values of the frequency and wavenumber of the wavetrains are localized in the qR region of the spectrum (see f.i. the point associated to the averaged normalized frequency an the associated normalized wavenumber in figure 3). Then, using the analytic approximation and the associated propagation parameters, we solved the IBVP for each hole following the procedure outlined in section 3.2 to obtain the corresponding numerical sequences . The typical computational time was of the order of 90 min using a computer with 512 cores for 40000 time-steps of ns and a discretization of the computational domain that leads to million of unknowns.

parameter expression unit sample mean
rad/m
(rad/s)
mm
s
m/s
equation (8) m/s
-
-
-
Table 2: Values of the propagation parameters obtained from the coefficients of the analytic approximation to the out-of-plane component of the incident field at the top surface of the plate. The dispersion around the sample mean is specified after the sign by the rounded-up value of the standard deviation of the mean.
Figure 11: Frame of the video sequence of the out-of-plane displacement field for the scattering of a transient narrow-band quasi-Rayleigh wave by a cylindrical hole with : reference field (experimental) on the left and problem field (numerical) on the right. Amplitude in units of . Lengths in millimeters. The axis numbering corresponds to that of figure 4. Full sequence in Video 1, MPEG, 1 MB, [video_1_figure_11.mp4].

To illustrate the obtained results in the comparison between experimental and numerical data we use representative sequences corresponding to experimental and numerical displacement fields for the cylindrical hole with nominal diameter (figure 11). A visual comparison of the sequences evidences a close agreement in the spatial distribution of the experimental and numerical values of the real displacement. This agreement is confirmed with a more careful visual comparison frame by frame as it is illustrated in figures 12 and 13 where, using a common grey level scale for the images and a common vertical scale for the profiles, the amplitude and phase of the experimental and numerical displacement are presented for two representative instants in the sequences: corresponding to the incident field before the interaction with the hole (figure 12), and corresponding to the scattered field after the interaction with the hole (figure 13). From these figures, apart from local differences associated to the fluctuations in the experimental field due to presence of speckle, we can conclude that there is a very good matching of the experimental and numerical values. This good matching can be found all along the sequence both for the incident field, which confirms that the employed approach for introducing the incident field using an analytic approximation works properly, and the scattered field, including the forward, lateral and backscattering areas in the image.

This agreement between experiment and numerics found in the visual comparison has been checked over different areas in the image in a more objective quantitative way using the relative rms global errors in L2 norm in amplitude

(21)

and in phase

(22)

being the number of pixels in the selected region and

(23)

the maximum values for amplitude and phase of the field over the ROI in the whole sequence. and were evaluated for all instants of the sequence in several sub-zones of interest (see figure 4) and the results (figure 14) show that the value of the relative rms global error is always lower than the typical value of the noise-to-signal ratio, in both amplitude and phase, in the experimental PTVH system (which is mainly associated to speckle fluctuations Tesis_PRG ()). In other words, a perfect match was obtained—up to a tolerance of the order of the experimental noise.

(a)
(b)
(c)
(d)
Amplitude(C1) Phase(C2)
Figure 12: Scattering of a transient narrow-band quasi-Rayleigh wave by a cylindrical hole with at : a) experimental field, (b) numerical field, (c) and (d) profiles along the lines in the and directions marked in white in (a) (thick) and (b) (thin). Amplitude in units of . Phase in radians. Lengths in millimeters. The axis numbering corresponds to that of figure 4.
(a)
(b)
(c)
(d)
Amplitude(C1) Phase(C2)
Figure 13: Scattering of a transient narrow-band quasi-Rayleigh wave by a cylindrical hole with at : a) experimental field, (b) numerical field, (c) and (d) profiles along the lines in the and directions marked in white in (a) (thick) and (b) (thin). Amplitude in units of . Phase in radians. Lengths in millimeters. The axis numbering corresponds to that of figure 4.
(a)
(b)
(c)
(d)
Amplitude(C1) Phase(C2)
Figure 14: Global relative total error in amplitude and phase between the problem field (numerical) and the reference field (experiment) for the complete sequence of the scattering of a transient narrow-band quasi-Rayleigh wave by cylindrical hole with for several regions within the ROI corresponding to those in figure 4: (a) , (b) , (c) and (d) . The arrows identify the time intervals in which there is an incident (i) or a scattered (s) acoustic field different from zero for each region.
Figure 15: Frame of a video sequence showing the components of the displacement vector in the scattering of a transient narrow-band quasi-Rayleigh wavetrain by a cylindrical hole with : (a) side view from the perspective of the negative axis, (b) side view from the perspective of the negative axis. The top side images show the out-of-plane component of the field in an area of centered at the hole within the ROI in figure 4 and the height of the side view is equal to the plate thickness mm. Axes orientation corresponds to those of figures 1 and 4. Amplitude in units of . Phase in radians. Full sequence in Video 2, MPEG, 2.5MB [video_2_figure_15.mp4].

To illustrate the potential of the FC(Gram) solver for quantitative characterization of temporal and spatial distribution of the acoustic field we present a full 3D sequence for the 12 mm hole that displays the temporal variation of all three components of the displacement vector in the interior volume of the plate (figure 15). In the backscattering area, the incident wave presents the characteristic spatial distribution of the amplitude of a quasi-Rayeligh wave, with significant values only in a layer of the order of one acoustic wavelength close to the top side of the plate—as clearly visible in the initial frames of the side view of the video corresponding to figure 15.a left. The scattered field, on the other hand, presents well-known features of scattering of elastic waves by defects ews133 (). Indeed, the reflected wave in the backscattering area is quasi-Rayleigh and it produces a standing wave pattern by interference with the incident field, as can be seen in the late frames of the aforementioned video. Further, the usual modal conversion to SH0 mode Dilligent () is evidenced by the large displacement that appears in a vicinity of the hole in figure 15 and in the time-interval in the associated video. The video also displays a wave perturbation that travels downward along the lateral side of the hole, which is subsequently scattered at the bottom, generating an upward perturbation which interferes with the downward perturbation and thus forms a standing wave pattern ews156 () in the lateral side of the hole (see side view of figure 15.b, center). At the same time, a scattered wave is generated which propagates close to the bottom side of the plate—which can be viewed in the final frames of the side view portions of the video for the and components and in figure 15.a left.

5 Conclusions

The scattering of transient narrow-band quasi-Rayleigh waves by through-thickness holes in aluminum plates in the high-frequency regime has been studied. Experimental sequences of scattering patterns measured with a self-developed full-field PTVH system have been compared with the corresponding simulated sequences derived by a state-of-the-art FC(Gram) elastodynamic solver. In addition to a broad agreement in the spatial distribution of the incident and scattered fields for every time instant, which is apparent by visual inspection, a detailed quantitative comparison between experiments and numerics has been presented. Consideration of the relative global errors observed show that a perfect match was obtained, both in amplitude and phase, up to a tolerance of the order of the experimental noise. It is suggested that the combined computation-and-measurement technique bears a significant potential for ultrasonic NDT applications in plate-like structures. Ongoing work focuses on generalization of these techniques to generally-shaped plate-like structures containing complex defects, such as slots, partly through-thickness holes, cracks, etc.

6 Acknowledgments

OB gratefully acknowledges support by NSF and AFOSR and DARPA through contracts DMS-1411876, DMS-1714169, FA9550-15-1-0043 and HR00111720035, and the NSSEFF Vannevar Bush Fellowship under contract number N00014-16-1-2808. The authors from the University of Vigo gratefully acknowledge support by the Spanish Ministerio de Economía y Competitividad and the European Comission (ERDF) in the context of the Plan Nacional de I+D+i (project number DPI2011-26163), by the Subdirección Xeral de Promoción Científica e Tecnolóxica universitaria da Xunta de Galicia in the context of the Plan Galego de I+D+i (project number 10PXIB303167PR) and supplementary co-funding from the Universidade de Vigo (contract number 09VIA07).

Appendix A Real and complex description of wave propagation

For any real displacement field described by equations (1-2) an associated complex field

(24)

with complex components

(25)

can be introduced in such a way that and , where stands for the real part of complex number and is the imaginary unit. The complex amplitude of each complex component is unambiguously determined by the acoustic amplitude and the acoustic phase of the corresponding real component and viceversa Born_Wolf (). Hence, there is a one-to-one correspondence between real and complex components and, as a consequence, between the real displacement field and the complex displacement field , which can be also written as where

(26)

is its associated complex amplitude.

Once the complex displacement is well-established, taking as a reference equations (3-4) we can introduce the components of the complex strain tensor as

(27)

and the complex components of the complex stress tensor in the linear elastic regime as

(28)

From the linearity of equations (3-4) and (