Vincent Dorval
email: vincent.dorval@cea.fr

Michel Darmon

Sylvain Chatillon

Larissa Fradkin

Simulation of the UT Inspection of Planar Defects Using a Generic GTD-Kirchhoff Approach

Keywords: Ultrasounds, Scattering, Kirchhoff approximation, GTD PACS: 43.20.El, 43.35.Bf, 43.35.Cg CRACK SCATTERING MODELS

The modeling of ultrasonic Non Destructive Evaluation often plays an important part in the assessment of detection capabilities or as a help to interpret experiments. The ultrasonic modeling tool of the CIVA platform uses semianalytical approximations for fast computations. Kirchhoff and GTD are two classical approximations for the modeling of echoes from plane-like defects such as cracks, and they aim at taking into account two different types of physical phenomena. The Kirchhoff approximation is mainly suitable to predict specular reflections from the flaw surface, whereas GTD is dedicated to the modeling of edge diffraction. As a consequence, these two approximations have distinct and complementary validity domains. Choosing between them requires expertise and is problematic in some inspection configurations. The Physical Theory of Diffraction (PTD) was developed based on both Kirchhoff and GTD in order to combine their advantages and overcome their limitations. The theoretical basis for PTD and its integration in the CIVA modeling approach are discussed in this communication. Several results that validate this newly developed model and illustrate its advantages are presented.

INTRODUCTION

The ultrasonic testing simulation module of the CIVA software offers ways of computing crack echoes, including semi-analytical models [START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for NDT[END_REF][START_REF] Darmon | [END_REF] and a hybrid model using finite element code (ATHENA code) [START_REF] Leymarie | A Semi-analytical FEM Hybrid Model for Ultrasonic Simulation of Complicated Wave Defect Interaction[END_REF]. In the CIVA 10 version, two semi-analytical models were available in the case of cracks: Kirchhoff and GTD [START_REF] Darmon | Recent modeling advances for ultrasonic TOFD inspections[END_REF]. They are based on different approximations and have complementary validity domains. Choosing between the two models required expertise. Moreover, in some inspection configurations, no model was entirely satisfactory. The PTD model, available in CIVA 11, aims at overcoming these limitations.

In this communication, the hypotheses of the Kirchhoff and GTD models are described in terms of Crack Opening Displacement. This formalism is then used to explain the principle of the more complete PTD model, which relies on both Kirchhoff and GTD. PTD simulated results are finally compared to both Kirchhoff and GTD simulations and in addition to finite element simulations.

Crack Opening Displacement

The Crack Opening Displacement (COD) is defined as the difference between the displacements on the two sides of a planar crack. It is a determining factor in the scattering of a planar crack. One way to make it appear is to develop the expression for the scattering signal given by Auld [START_REF] Auld | [END_REF] in the frequency domain:

𝛿Γ 𝑏𝑎 = 1 4𝑃 ∫ (𝑣 1 ⃗⃗⃗⃗ • 𝑇 2 -𝑣 2 ⃗⃗⃗⃗ • 𝑇 1 ) • 𝑛 ⃗ 𝑑𝑆 𝑆 𝐹 , (1)
where 𝛿Γ 𝑏𝑎 is the difference in the measured signal due to the presence of the crack, P is the excitation power, 𝑣 a particle velocity, 𝑇 a stress, 𝑆 𝐹 the flaw surface and 𝑛 ⃗ the local outward normal to that surface. The index 1 denotes the field obtained when the emitting probe emits and the flaw is absent. The index 2 denotes the field obtained when the receiving probe emits and the flaw is present.

We assume that the surface of a planar crack behaves as a free surface. Therefore, the presence of the flaw implies that the stress projected on the normal to the surface is zero. This leads to the disappearance of the first term of equation [START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for NDT[END_REF]. The COD appears in the remaining term when the contributions on the two sides on the crack are combined: the stress 𝑇 1 in the absence of the crack is continuous and the normal 𝑛 ⃗ has opposed directions on each sides. It is then possible to make the COD appear:

𝛿Γ 𝑏𝑎 = - 𝑖𝜔 4𝑃 ∫ 𝐶𝑂𝐷 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝑇 1 • 𝑛 ⃗ 𝑑𝑆 𝑆 + , (2)
S+ is one of the sides of the crack, arbitrary chosen, ω is the angular frequency and the COD is defined here as the difference between the particle displacement on the side S+ and on the other side at the same location.

The focus of this communication is the computation of crack echoes. The ultrasonic fields in the absence of the crack are known. In the applications presented here, they are calculated in CIVA software using a paraxial beam method [6]. In this framework, the modeling of crack scattering can be reduced to the calculation of the COD. In the next sections, the Kirchhoff and GTD models are described with an emphasis on the corresponding CODs.

The Kirchhoff approximation

The Kirchhoff approximation consists, for each point of the defect surface, in assuming that the surface is an infinite plane tangent to the actual surface. The ultrasonic field at this point is expressed as the sum of an incident plane wave and the plane waves reflected by the infinite plane. This approximation is also sometimes called the "tangent plane approximation". The ultrasonic field is supposed to be zero on the non-insonified side of the crack, therefore the COD coincides with the displacement on the insonified side.

Coefficients for different incident and diffracted waves and angles for the scattering of ultrasonic waves by cracks are given by Chapman [START_REF] Chapman | Ultrasonic scattering from smooth flat cracks: an elastodynamic Kirchhoff diffraction theory[END_REF]. The current CIVA implementation of Kirchhoff model computes the integral of equation ( 2) numerically based on the amplitudes of incident and reflected waves [START_REF] Dorval | A general Kirchhoff approximation for echo simulation in ultrasonic NDT[END_REF], instead of using expressions such as those of Chapman where the amplitudes of these waves are implicit. This approach has the advantage of being easily adapted to any material.

The Kirchhoff expression of the echo in an integral form over the flaw surface allows a lot of flexibility in practical applications, as the variations of the fields over the surface or complex flaw geometries can easily be taken into account.

The Kirchhoff approximation implies that the COD can be non-zero at the edge of a crack. This is non-physical, as the two sides of a crack join at the edge, therefore there cannot be any difference between their displacements at such locations. This leads to errors in the calculated diffraction when the contribution of the edge is significant. As a consequence, the Kirchhoff approximation will tend to be inaccurate in the calculation of edge diffraction.

Conversely, it will tend to give reliable results for specular echoes, where the contributions of the crack surface dominate compared to the contributions of the edges.

The Geometrical Theory of Diffraction

The Geometrical Theory of Diffraction (GTD) is based on an exact solution for the COD of a semi-infinite crack. As a consequence, it is complementary to the Kirchhoff approximation: it is appropriate for modeling edge diffraction, but not for specular echoes.

Notable early works on the application of GTD to elastodynamics have been done by Achenbach and Gautesen [9]. For the purpose of an interpretation in terms of COD, another way of obtaining coefficients [10] is more convenient. This later approach makes use of expressions similar to the integral of equation [START_REF] Darmon | [END_REF]. For incident plane waves, the amplitude of 𝑇 1 is an oscillating exponential and the integral can be treated as a spatial Fourier transform of the COD. This Fourier transform is calculated asymptotically using the stationary phase method. The asymptotic calculation is accurate for most wave directions, except for those near critical angle where it can cause errors and for those around specular and transmission directions.

The main limitation of the GTD coefficients is their divergence around specular and transmission directions, which is related to the semi-infinite crack assumption.

The Physical Theory of Diffraction

There are several ways of overcoming the limitations of the Kirchhoff and GTD models. One is to improve the GTD model in order to correct its results in the specular and diffraction regions. This leads to methods such as the Uniform Asymptotic Theory (UAT) of diffraction [11][START_REF] Achenbach | Ray methods for waves in elastic solids: with applications to scattering by cracks[END_REF][START_REF] Fradkin | [END_REF]. These methods constitute an improvement over GTD but they are inconvenient for handling complex geometries and taking into account the variations of the fields over the surface of the flaw.

The Physical Theory of Diffraction (PTD) is another way to overcome these limitations. It has been proposed under several forms by several authors: It was first developed in acoustics [START_REF] Ufimtsev | Fundamentals of the Physical Theory of Diffraction[END_REF], and later in elastodynamics [START_REF] Mcmaken | [END_REF]. The work presented here is based on another elastodynamic form proposed in [16]. It allows combining all the advantages of the Kirchhoff (suitable for specular and echoes, convenient for complex geometries and field variations) and GTD (suitable for diffraction echoes) models. The PTD echo of a 2D planar crack can be decomposed in the following way:

𝐸𝑐ℎ𝑜 𝑃𝑇𝐷 = 𝐷𝑖𝑓𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝐺𝑇𝐷 𝐸𝑑𝑔𝑒1 -𝐷𝑖𝑓𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝐾𝑖𝑟𝑐ℎℎ𝑜𝑓𝑓 𝐸𝑑𝑔𝑒1 + 𝐸𝑐ℎ𝑜 𝐾𝑖𝑟𝑐ℎℎ𝑜𝑓𝑓 + 𝐷𝑖𝑓𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝐺𝑇𝐷 𝐸𝑑𝑔𝑒2 -𝐷𝑖𝑓𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝐾𝑖𝑟𝑐ℎℎ𝑜𝑓𝑓 𝐸𝑑𝑔𝑒2 , (3)
where EchoKirchhoff and DiffractionGTD are the outputs of the Kirchhoff and GTD model. DiffractionKirchhoff is a correction term that ensures the consistency of the model. This term corresponds to the diffraction by a semi-infinite crack simulated in the Kirchhoff approximation: it can be obtained by applying the same asymptotic approach as for the GTD model, except the exact COD is replaced by the Kirchhoff approximated COD.

DiffractionKirchhoff has two notable properties. Around specular directions, it diverges in the same way as the GTD coefficient. For observation directions for which edge diffraction is predominant, it agrees with standard Kirchhoff coefficient. As a result, in expression (3), the two DiffractionKirchhoff terms will cancel EchoKirchhoff in diffraction directions and they will cancel the singularities of the two DiffractionGTD terms in specular directions. PTD therefore agrees with GTD for diffraction echoes and with Kirchhoff for specular echoes.

An interpretation in terms of COD clarifies the way PTD combines Kirchhoff and GTD models. As the echoes are proportional to the COD, the echo decomposition of equation (3) translates into a COD decomposition:

𝐶𝑂𝐷 𝑃𝑇𝐷 = 𝐶𝑂𝐷 𝐺𝑇𝐷 𝐸𝑑𝑔𝑒1 -𝐶𝑂𝐷 𝐾𝑖𝑟𝑐ℎℎ𝑜𝑓𝑓 𝐸𝑑𝑔𝑒1 + 𝐶𝑂𝐷 𝐾𝑖𝑟𝑐ℎℎ𝑜𝑓𝑓 + 𝐶𝑂𝐷 𝐺𝑇𝐷 𝐸𝑑𝑔𝑒2 -𝐶𝑂𝐷 𝐾𝑖𝑟𝑐ℎℎ𝑜𝑓𝑓 𝐸𝑑𝑔𝑒2 , ( 4
)
where the COD associated to an edge is the COD of a semi-infinite crack starting at the corresponding edge. The case of a 2D crack in normal incidence (figure 1) will be considered as an illustration. In this simplified case, the COD corresponding to the Kirchhoff approximation (figure 2(a)) is uniform along the crack. As expected with this approximation, it has an unrealistic discontinuity at the edges of the crack, where the COD should be zero. Figure 2(b) and 2(c) shows the COD corresponding to the Kirchhoff approximation for a semiinfinite crack starting at the left edge. Figure 2(d) and 2(e) show the exact solution for the COD corresponding to a semi-infinite crack starting at the left edge [START_REF] Borovikov | Uniform Stationary Phase Method[END_REF]: it illustrates the exact solution for the COD, which is implied but not explicitly given in the GTD model. Figure 2(f) is the COD corresponding to the sum of the 5 previous CODs, as expressed in equation ( 4). It illustrates how the corrective terms cancel the undesirable effects of each model: the discontinuity due to Kirchhoff disappears, and so do the semi-infinite contributions of GTD. Only small residual oscillations due to GTD remain.

For more complex crack geometries or for non-homogeneous fields, the treatment of the edges is similar as in the simple case: the crack is still assumed to be semi-infinite, whereas the field is assumed to be constant and to have the same value as at the edge everywhere. However, the surface integral of the Kirchhoff approximation is done by taking into account both the specific geometry of the crack and the spatial variations of the field.

In the CIVA implementation of this method, the surface and the edges of the crack are both meshed. The surface mesh is used for the computation the CODKirchhoff term in equation ( 4), and the edge mesh the DiffractionGTD-DiffractionKirchhoff terms. Figure 3 illustrates how defects can be meshed.

COMPARISON OF RESULTS OBTAINED WITH THE DIFFERENTS MODELS

In this section, simulated echoes obtained for a crack seen from different incident angles in a pulse-echo configuration are shown. They will illustrate how PTD behaves compared to Kirchhoff and to GTD for different directions, and compared to a finite elements result. Comparison of echoes obtained using the PTD model and a finite elements approach FIGURE 5. Variations of the echo amplitude in the pulse-echo configuration of Figure 4 for a wider range of incident angles.

Comparison echo calculated with the three models

The horizontal axis is the angle (in degrees) corresponding to the incidence on the crack, 90° being the normal incidence.

The results of figure 5 have been obtained by computing the echo in the configuration of figure 4 for a wider range of probe positions around the cylinder. Computations have been performed using the CIVA implementation of the PTD model and the CIVA-Athena finite element module [START_REF] Leymarie | A Semi-analytical FEM Hybrid Model for Ultrasonic Simulation of Complicated Wave Defect Interaction[END_REF]. A satisfactory agreement is obtained for all angles, regardless of whether the predominant phenomenon is edge diffraction or specular reflection. Significant differences, such as those of figure 4, would have been obtained for diffraction angles between PTD and Kirchhoff and for specular angles between PTD and GTD.

CONCLUSION

The PTD model combines the Kirchhoff and GTD models to cover the validity domain of both of them. The interpretation in terms of Crack Opening Displacement presented in this paper provides insight in the hypothesis of the model and could be used as a framework for new model developments.

PTD has been implemented in the CIVA 11 software as a way to improve the simulation of crack echoes. It makes the results more accurate in cases where Kirchhoff or GTD were not sufficient, and removes the need for the user to select the appropriate echo model. Limitations remain in cases that could not be treated by Kirchhoff or by GTD (namely, scattering at critical angles or echoes of very small cracks) and are of interest for additional developments.

FIGURE 1 .FIGURE 2 .

 12 FIGURE 1. Simplified case for the illustrations of the COD. The crack size has been taken as unity.

FIGURE 3 .

 3 FIGURE 3. Examples of crack meshes for the calculation of echoes in the PTD approximation. Red dots are the mesh points for the Kirchhoff model, blue crosses are the mesh points for the edge diffraction. (a): rectangular crack, (b): circular crack.

FIGURE 4 .

 4 FIGURE 4. Two measurement pulse-echo configurations and echo signals for a 5mm crack in a cylinder inspected using compressional waves, in a 2D configuration. (a): normal incidence configuration, (b): normal incidence echoes, (c): 45° incidence configuration, (d): 45° incidence echoes. For (b) and (d), the solid black lines are the PTD results, the red dashed lines (larges dashes) are the Kirchhoff results and the blue dashed lines (small dashes) are the GTD results. The echoes of figure 4 (b) correspond to specular reflection. As expected, PTD agrees with Kirchhoff and not with GTD. It should be noted that for practical use in CIVA, the divergence of the GTD coefficients is attenuated by preventing the coefficient to rise too much around specular directions. It explains why the GTD result is low. An unaltered GTD model would yield an echo with extremely high amplitude, and would still not agree with PTD. The echoes of figure 4 (d) correspond to edge diffraction. As expected, PTD agrees with GTD and not with Kirchhoff.


Anastassios Skarlatos
email: anastasios.skarlatos@cea.fr

Member, IEEE Theodoros Theodoulidis

Senior Member, IEEE Nikolaos Poulakis

A fast and robust semi-analytical approach for the calculation of coil transient eddy-current response above planar specimens

published or not. The documents may come

A fast and robust semi-analytical approach for the calculation of coil transient eddy-current response above planar specimens

Anastassios Skarlatos, Theodoros Theodoulidis, Nikolaos Poulakis

I. INTRODUCTION

Pulsed eddy-current testing (PECT) has attracted the interest of numerous researchers and non-destructive testing practitioners the last years for a number of reasons. The main motivation of this increased interest is the attempt of extending the sensitivity zone of traditional harmonic ECT techniques deeper in the material. Pulsed currents, having in principle broadband spectra, can respond to this requirement, without compromising the resolution for subsurface defects thanks to the higher frequency spectral components. From the point of view of laboratory implementation, PECT can be carried out with relatively simple equipment. This convenience is also enhanced by the latest developments of fast data acquisition cards nowadays.

From the theoretical point of view, this raised interest in the PECT method has motivated the development of fast and reliable theoretical models for the prediction of the measurement signals. Semi-analytical solutions are always attractive, when the problem under consideration is amenable to analytic formulation, due to very interesting computational times with respect to numerical solutions. There is a huge amount of published works on the semi-analytical modelling of harmonic ECT in canonical multilayer structures, stemming from the seminal articles of Dodd and Deeds [START_REF] Dodd | Analytical solutions to eddy current probe coil problems[END_REF], [START_REF] Dodd | Induction coils coaxial with an arbitrary number of cylindrical conductors[END_REF]. Taking advantage of this corpus of published works, transient solutions can be constructed by addressing the eddy-current problem under a monochromatic excitation using the above mentioned established harmonic solutions and resorting to an inverse Fourier or Laplace transform to obtain the corresponding time signals.

The main difference of the two approaches is the domain of the sampling frequencies, which in the case of the Laplace transform spans the half of the complex plane, whereas the Fourier frequencies lie on the imaginary axis.

The Fourier transform is perhaps the most straight forward approach for the calculation of the transient response due to the simplicity and robustness of the inverse transform operator, which in addition can take advantage of the well established Fast Fourier Transform (FFT) algorithms [START_REF] Li | Fast analytical modelling for pulsed eddy current evaluation[END_REF]- [START_REF] De Haan | Analytical expressions for transient induction voltage in a receiving coil due to a coaxial transmitting coil over a conducting plate[END_REF]. Even if FFT can significantly accelerate the calculation of the time signals, the increased number of frequencies required for the accurate representation of the usual transient signals met in PECT applications can inhibit the inherent excellent computational performances of the semi-analytical solutions. To overcome this problem, intelligent sampling of the frequency axis in combination with interpolation has been proposed in recent works [START_REF] Xie | Efficient numerical solver for simulation of pulsed eddy-current testing signals[END_REF]- [START_REF] Xie | Sizing of wall thinning defects using pulsed eddy current testing signals based on a hybrid inverse analysis method[END_REF]. An insurmountable difficulty inherit with the Fourier approach is the so-called Gibbs effect apparent at the vicinity of signal discontinuities, which is the case of step or pulse excitations, frequently used in PECT configurations.

Laplace transform on the other side remains the most efficient and elegant tool for the solution of initial value problems. The main challenge, however, of this approach comes from the lack of a generic and robust algorithm for the sampling of the complex plane and the thereupon linked inversion approach. As a result, the solution to the problem remains to a great extend ad hoc, with a number of works devising closed-form (whenever possible) inversion formulas [START_REF] Bowler | Pulsed eddy-current response to a conducting half-space[END_REF], [START_REF] Fu | Transient eddy-current driver pickup probe response due to a conductive plate[END_REF] or pole extraction using Padé approximation or decomposition using the generalised pencil of function (GPoF) approach [START_REF] Theodoulidis | Developments in calculating the transient eddy-current response from a conductive plate[END_REF], [START_REF] Skarlatos | A modal approach for the solution of the non-linear induction problem in ferromagnetic media[END_REF]. For the short-time transient signal calculation, which has the form of exponential damping, a number of specialised numerical algorithms such as the Stehfest method has also been proposed with excellent results [START_REF] Theodoulidis | Efficient calculation of transient eddy current response from multilayer cylindrical conductive media[END_REF].

A more direct approach, usually applied to mesh-based approaches, such as the finite element method (FEM) or the finite integration technique (FIT), relied in discretisation of the temporal derivatives and iterative solution using a timestepping scheme [START_REF] Clemens | Transient eddy-current calculation with the FI-method[END_REF]- [START_REF] Dutiné | Explicit time integration of transient eddy current problems[END_REF]. Stability issues impose the use of implicit time integration schemes for diffusion problems as the one examined here. Nevertheless, recent works have shown that quasi-explicit schemes can also be applied with the consequent benefits in computational performances [START_REF] Dutiné | Explicit time integration of transient eddy current problems[END_REF]. This latter approach has been barely considered, to the best of he authors' knowledge, in combination with modal schemes, and the reason lies with the right-hand term in the diffusion equation linked with the previous timestep(s). Such inhomogeneous terms may be resolved with the extension of the modal projection basis, as it was shown in the case of non-linear eddy-current problem [START_REF] Skarlatos | A modal approach for the solution of the non-linear induction problem in ferromagnetic media[END_REF], [START_REF] Skarlatos | Study of the non-linear eddy-current response in a ferromagnetic plate: Theoretical analysis for the 2D case[END_REF].

In this work, the ideas developed in [START_REF] Skarlatos | A modal approach for the solution of the non-linear induction problem in ferromagnetic media[END_REF], [START_REF] Skarlatos | Study of the non-linear eddy-current response in a ferromagnetic plate: Theoretical analysis for the 2D case[END_REF] are applied for the solution of the implicit time-stepping schemes obtained from the temporal discretisation of the eddy-current equation of state in the spectral domain for the solution of the PECT inspection of a conducting plate. The interest in combining the spectral (semi-analytical) approach with time-stepping relies in the robustness of the approach regardless the form of the excitation signal. It is shown that the state equation obtained upon time-discretisation bears the same form with the harmonic case with an effective frequency (linked with the discretisation rate), which is the same for all time-steps (provided that a uniform discretisation has been applied). This detail is important for the efficiency of the overall scheme since combined with the sparsity of the system matrix obtained upon the projection to the spectral basis (the problem modes) is translated to a problem factorisation. A time-domain version of the reciprocity theorem is then proposed in order to calculate the induced electromotive force (EMF) at the pick-up coil without the need of passing from the Fourier domain.

The paper is organised as follows. The considered problem is presented in section II. The mathematical formulation of the approach is developed in the next section (section III), followed by the calculation of receiving coil EMF in section IV. An analysis of the numerical performance of the method is presented in section V. Finally, the numerical results obtained by the proposed approach are compared against reference solutions in section VI.

II. PROBLEM DESCRIPTION

The problem which we shall consider is depicted in Fig. 1. An infinite conducting plate of thickness d is inspected by a couple of cylindrical coaxial coils, whose axes are perpendicular to the surface of the plate. The plate conductivity and magnetic permeability is equal to σ and µ = µ r µ 0 , respectively, with µ 0 standing for the permeability of the free space and µ r is the relative permeability. The coils are connected in driver-pickup mode, that is, current is supplied to the first one (driver coil), and the induced voltage at the second (pickup coil) is measured.

The excitation current has arbitrary temporal profile I = I(t), whose energy is principally distributed in the lowfrequency band between 0 and 100 kHz. Under this assumption, the quasi-static approximation for the Maxwell equations is valid, and we can proceed to the analysis in a similar way that we do for the typical harmonic eddy-current problem. In addition, the coil current is assumed to be weak enough in order that non-linear effects be negligible.

For the needs of the analysis, the three layers of the geometry (plate and air columns above and underneath it) will be referred to using the indices 1, 2 and 3, as shown in Fig. 1. The objective of the analysis is to calculate the induced voltage in the coil V (t).

z z = -d σ, µ 1 2 3 ρ = ρ L z 1p z 2p z 1d z 2d r 1p r 2p r 1d r 2d z = 0
It is convenient to artificially truncate the computational domain at a sufficiently remote from the coil distance, which will be denoted as ρ r L and apply a perfectly electric conductor (PEC) condition there. Since we are dealing with a diffusion problem, the field intensity is rapidly decreasing as we are moving far from the source, and hence it can be considered as negligible after a certain distance. This is the basic assumption adopted in the so called truncated region eigenfunction (TREE) approach, used in this work and which proves to result in a convenient simplification of the problem (with respect to the classical Dood and Deeds approach [START_REF] Dodd | Analytical solutions to eddy current probe coil problems[END_REF], [START_REF] Dodd | Induction coils coaxial with an arbitrary number of cylindrical conductors[END_REF]).

III. SOLUTION IN THE TIME DOMAIN USING TIME DESCRETISATION

A. Vector potential formulation in time domain

The mathematical analysis of the transient problem follows the A -Φ formulation for the quasi-static limit, which is the established approach for the treatment of axisymmetric problems. The magnetic vector potential A, and the electric scalar potential Φ are defined via the relations

B = ∇ × A (1)
and

E = - dA dt -∇Φ (2)
with E, B being the electric field and the magnetic flux density, respectively. Notice that for problems with axisymetric summetruy we can adopt a potential ansatz of the type A = Ae φ , where e φ is the unit vector along the φ direction, i.e. the magnetic vector potential lies in the azimuthal direction. Upon substitution to the Maxwell equations one obtains the diffusion equation for the potentials

∇ × µ -1 ∇ × A + σ dA dt = -σΦ + J ( 3
)
where J is the excitation current density, and σ, µ stand for the electrical conductivity and the magnetic permeability of the medium. The material will be considered as isotropic and linear throughout this work, i.e. the magnetic flux density B and the magnetic field H are related via a linear constitutive relation

B = µH. (4)
Assuming in addition a piecewise homogeneous material and using the Coulomb's gauge for the magnetic vector potential A, the magnetic permeability can be drawn outside the curl operator reducing (3) to the scalar equation

∇ 2 - 1 ρ 2 -µσ d dt A = -µJ, (5)
where we have taken into account the azimuthal orientation of the magnetic potential and the excitation current (J = Je φ ). We are interested in the source-free case J = 0, since all the layers except the internal and/or the external region of the tube wall do not contain excitation currents. The coil field in these latter two regions will be treated indirectly as we will see later.

Equation ( 5) is a parabolic equation, and hence its numerical integration should be carried out using an unconditionally stable scheme. Let us consider a homogeneous descretisation of the time axis t i = i∆t, with i = 0, 1, . . . and ∆t the discretisation step. Using an implicit Euler scheme for the time derivative, (5) reduces to the following update equation, where the solution at i is given by means of the solution at the previous time-step i -1

∇ 2 - 1 ρ 2 - µσ ∆t A i = - µσ ∆t A i-1 . (6)
The notation can be simplified by introducing the equivalent angular frequency ω s = -√ -1/∆t, with √ -1 being the imaginary unit, and (6) becomes

∇ 2 - 1 ρ 2 -k 2 A n = -f n . (7)
with k 2 = √ -1ω s µσ and f n = k 2 A n-1 . Equation ( 7) is formally equivalent with the inhomogeneous Helmholtz equation for frequency equal to the discretisation (sampling) frequency of the time integration scheme. Our problem thus reduces to the solution of the Helmholtz equation with a source term at each time-step, which is a very convenient observation since one can resort to tools developed for the harmonic regime in a planar medium [START_REF] Skarlatos | A modal approach for the solution of the non-linear induction problem in ferromagnetic media[END_REF], [START_REF] Skarlatos | Study of the non-linear eddy-current response in a ferromagnetic plate: Theoretical analysis for the 2D case[END_REF]- [START_REF] Theodoulidis | Eddy-current interaction of a long coil with a slot in a conductive plate[END_REF].

B. General solution of the inhomogeneous Helmholtz equation

For the solution of (7), we shall adopt the same strategy used in [START_REF] Skarlatos | A modal approach for the solution of the non-linear induction problem in ferromagnetic media[END_REF], [START_REF] Skarlatos | Study of the non-linear eddy-current response in a ferromagnetic plate: Theoretical analysis for the 2D case[END_REF]. The solution of the inhomogeneous equation is decomposed in the sum of a partial solution, which eliminates the right-hand-side term, and the solution of the homogeneous equation, which takes care of the current excitation and assures the fulfilment of the continuity relations across the geometry interfaces.

In the two air regions, above and underneath the plate both k and f i in the diffusion equation ( 6) vanish, which means that the problem there reduces to the magnetostatic formulation 1 . The solution for the magnetic potential in these two regions reads

A (1) i (ρ, z) = ∞ m=1 J 1 (κ m ρ) C (1) i;m e κmz + D (s) i;m e -κmz (8)
and

A (3) i (ρ, z) = ∞ m=1 C (3) i;m J 1 (κ m ρ) e -κmz (9)
respectively. Note that D (s)

im stand for the development coefficients for the coil field in the unbounded medium, and they are hence known. Their explicit expressions will be given at a later point. The κ m eigenvalues are determined by the zero tangential potential condition imposed by the PEC boundary at ρ L , and they are obtained by the roots of the Bessel functions

J 1 (κ m ρ L ) = 0, m = 1, 2, . . . , ∞. (10)
Inside the plate, we must address the inhomogeneous Helmholtz equation. The general solution can be expressed as the superposition of the homogeneous solution and a special solution, which satisfies the right-hand-side [START_REF] Skarlatos | Study of the non-linear eddy-current response in a ferromagnetic plate: Theoretical analysis for the 2D case[END_REF], [START_REF] Skarlatos | Calculation of the eddy-current flow around a cylindrical through-hole in a finite-thickness plate[END_REF]

A (2) i (ρ, z) = A (h) i (ρ, z) + A (s) i (ρ, z) (11)
The homogeneous solution reads

A (h) i (ρ, z) = ∞ m=1 J 1 (κ m ρ) C (2)
i;m e vmz + D

i;m e -vmz [START_REF] Skarlatos | A modal approach for the solution of the non-linear induction problem in ferromagnetic media[END_REF] with v 2 m = κ 2 m + k 2 . The modal basis used for the homogeneous solution belongs to the kernel of the Helmholtz operator (ker ∇ 2 -k 2 ) and hence must be complemented in the z direction for the development of the special solution, which must satisfy a nonzero right-hand side [START_REF] Skarlatos | Start vector generation for implicit Newmark time integration of the wave equation[END_REF], [START_REF] Skarlatos | Study of the non-linear eddy-current response in a ferromagnetic plate: Theoretical analysis for the 2D case[END_REF]. A suitable eigenbasis for this development is the Fourier basis. Since the solution does not satisfy any symmetry property along z the full Fourier basis must be used with period slightly larger than d (given that a Fourier series describes a periodic function). Such choice, however, would suffer from the Gibbs effect at the plate boundaries. To overcome this difficulty, the special solution can be split into two terms: one satisfying zero boundary conditions at the plate interfaces, and a second term accounting for the in general non-zero field value at those interfaces. This representation is schematically shown in Fig. 2. We can thus write

A (s) i (ρ, z) = ∞ m=1 J 1 (κ m ρ) ∞ n=1 [c i;mn sin(q n z) + g i;m (z)] (13)
with eigenvalues q n being obtained by the zero field condition at the plate interfaces, namely

q n = nπ d , n = 1, 2, . . . , ∞. (14)
We are free to chose the form of the g i;m (z) with the condition that it contains a sufficient number of degrees of freedom in order to account for the field values at z = 0 and z = -d. These conditions are easily met by any polynomial of degree greater than 0. Another possible choice is

g i;m (z) = 2 n=1 d i;mn cos(q n z) ( 15
)
with q n given by ( 14) and the coefficients d i;mn are determined by matching the boundary values for g i;m (z). In this work, this latter choice has been adopted since it leads to slightly simpler equations.

Fig. 2: Representation of an arbitrary function as a sum of a sine series and a residual term.

C. Application of the update equation: special solution calculation

Substituting (11) into (6), the first two terms vanish by construction (being solutions of the homogeneous equation) yielding

∞ m=1 J 1 (κ m ρ) ∞ n=1 v 2 m + q 2
n c i;mn sin(q n z)

+ 2 n=1 v 2 m + q 2 n d i;mn cos(q n z) = f i (ρ, z) . (16)
The d i;mn coefficients are easily obtained by observing (16) at the two plate interfaces (thus the first term vanishes) and weighting the two sides of the equation by the radial eigenfunctions J 1 (κ m ρ)

d i;mn = J 1 (κ m ρ) , f i (ρ, 0) ∓ f i (ρ, -d) 2 (v 2 m + q 2 n ) , n = 1, 2 (17)
with the minus sign being chosen for n = 1 and the plus for n = 2. The inner product is defined as

J 1 (κ m ρ) , R(ρ) := 2 ρ 2 L J 2 0 (κ m ρ 2 L ) ρ L 0 ρJ 1 (κ m ρ) R(ρ) dρ. ( 18
)
Recalling the definition of f i , which stands for the previous timestep solution, i.e.

f i = k 2 A (2)
i-1 and substituting the corresponding modal development, (17) reduces to

d i;mn = k 2 v 2 m + q 2 n d i-1;mn + 1 2 1 ∓ e -vmd C (2) i-1;m + 1 2 1 ∓ e vmd D (2) i-1;m , n = 1, 2. ( 19
)
Note that C i-1;m , D

i-1;m and d i-1;mn are known since they have been calculated in the previous step.

Having evaluated a i;m and b i;m , we can move this term to the right-hand-side of ( 16) and weight both terms with the Laplacian operator eigenfunctions J 1 (κ m ρ) sin(q n z) thus obtaining

c i,mn = 1 v 2 m + q 2 n J 1 (κ m ρ) sin(q n z) , f i (ρ, z) , n = 1, . . . , ∞ (20)
with the inner product along ρ being given by ( 18) and the inner product along z defined as

sin(q n z) , Z(z) := 2 d 0 -d sin(q n z) Z(z) dz. (21)
Recalling again the f i definition as the solution at the previous timestep, [START_REF] Auld | Eddy current probe response to open and closed surface flaws[END_REF] becomes

c i,mn = k 2 v 2 m + q 2 n c i-1;mn + C (2) i-1;m M mn + D (2) i-1;m N mn , n = 1, . . . , ∞ (22)
with the integrals

M mn N mn = 2 d 0 -d sin(q n z) e ±vmz dz = - 2 d q n v 2 m + q 2 n 1 -e ∓vmd (-1) n . (23)

D. Application of the continuity relations: homogeneous solution calculation

The development coefficients of the homogeneous solution

C (1) i;m , C (2) i;m , D (2)
i;m and D [START_REF] Li | Fast analytical modelling for pulsed eddy current evaluation[END_REF] i;m will be determined by the application of the continuity relations at the interfaces of the geometry. We have two conditions in our procession which the solution must satisfy: the continuity of the magnetic potential and its normal derivative. Hence at the z = 0 interface we have

A (1) ρ1 = A (2) ρ1 (24)
∂A (1) ∂z ρ1 = 1 µ r ∂A (2) ∂z ρ1 [START_REF] Cheo | A reciprocity theorems for electromagnetic fields with general time dependence[END_REF] and similarly for z = -d.

Substituting the expressions ( 8)-( 11) and weighting with the radial eigenfuction using [START_REF] Theodoulidis | Eddy Current Canonical Problems (with applications to nondestructive evaluation)[END_REF], we obtain for z 0

C (1) i;m + D (s) i;m = C (2) i;m + D (2) i;m + 2 n=1 d i;m1 (26) µ r κ m C (1) i;m -D (s) i;m = v m C (2) i;m -D (2) i;m + ∞ n=1 q n c i;mn . ( 27
)
In the same fashion the continuity relations at z = -d yield

C (3) i;m = C (2) i;m + D (2) i;m + 2 n=1 (-1) n d i;mn . ( 28
) µ r κ m C (3) i;m = v m C (2) i;m -D (2) i;m + ∞ n=1 (-1) n q n c i;mn . ( 29
)

E. Source field

In the previous development, the field of the source coil in the air, represented in the modal solution via the coefficients D (s) i;m has been considered as known. In order to calculate these coefficients, one has to recall that the coil field will satisfy the magnetostatic equations, since σ = 0 in air, and hence (5) reduces to

∇ 2 - 1 ρ 2 A 0 (ρ, z, t) = -µJ(ρ, z, t) . (30)
It becomes clear from [START_REF] Theodoulidis | Impedance of an induction coil at the opening of a borehole in a conductor[END_REF] that the solution A 0 follows instantaneously the current density value, which in its turn can be split in a spatial and a temporal part.

J(ρ, z, t) = ι(ρ, z) I(t) (31)
with

ι(ρ, z) = N d (r 2d -r 1d )(z 2d -z 1d ) , in coil. 0 elsewhere (32)
N d stands for the number of turns in the driving coil, and r 1d , r 2d , z 1d , z 2d its radial and axial dimensions, shown in Fig. 1. Following [START_REF] Theodoulidis | Eddy Current Canonical Problems (with applications to nondestructive evaluation)[END_REF], D

i;m are given from the relation

D (s) i;m = -4µ 0 ι(ρ, z) sinh κ m l 2 χ(κ m r 1d , κ m r 2d ) κ 5 n [ρ L J 1 (κ n ρ L )] 2 × J m (κ n ρ d ) e -κmz d I i ( 33
)
where r d = (r 2d + r 1d )/2 and z d = (z 2d + z 1d )/2 the coordinates of the coil centre and l = z 2d -z 1d its length. I i stands for the discretised current signal, i.e. I i = I(t i ).

IV. CALCULATION OF THE PICK-UP COIL VOLTAGE

For the driver-pick-up coil configuration that we examine here, an important quantity is the electromotive force (EMF) induced at the pick-up coil, for the reason that it constitutes the main experimental variable for this type of control. Its calculation is thus one of the major objectives of the mathematical modelling.

In harmonic regime, the respective quantity (or the equivalent variable of the coil mutual impedance) is calculated in an elegant way using the Auld's reciprocity theorem [START_REF] Auld | Theoretical characterization and comparison of resonantprobe microwave eddy-current testing with conventional low-frequency eddy-current methods[END_REF], [START_REF] Auld | Eddy current probe response to open and closed surface flaws[END_REF]. For transient measurements, this theorem can be readily extended using convolution integrals. The TD case has been treated principally in the antennas community for the calculation of receiver response to arbitrary signals [START_REF] Welch | Reciprocity theorems for electromagnetic fields whose time dependence is arbitrary[END_REF]- [START_REF] Shlivinski | Time-domain coupled response of antennas -reciprocity theorem approach[END_REF]. To the best of the authors' knowledge, this problem has been scarcely explored for low-frequency problems and therefore it will be considered herein in some detail.

The general form for the 2D problem with rotational symmetry reads

∆V (t) * I(t) = 1 µ 0 ∂V A ec * ∇ dA s dt - dA s dt * ∇A ec • dS (34)
It is obtained following the standard procedure found in almost all the classical textbooks. A detailed derivation is provided in the appendix for the sake of the text self-completeness. Interestingly, the convolution integral of (34) right-handside can be dropped by observing that the source field, being solution of the magnetostatic problem, instantly follows the time variations of the input current. Note that this is an approximation which holds only in the context of the low-frequency quasi-static formulation, which the eddy-current problem satisfies. By virtue of this property, the time derivatives of the source field will provide Dirac functions in case that the current excitation is a Heaviside function H(t). Equation (34) becomes

∆s(t) = 1 µ 0 ∂V [A ec (t) ∇A * s -A * s ∇A ec (t)] • dS (35)
where A * s is the source field for a unit constant current and ∆s is the EMF integral

s(t) = t 0 ∆V (τ ) dτ. (36)
The EMF expression for an arbitrary current profile I(t) can be easily obtained from the step response ∆s(t) by applying the Duhamel's integral [START_REF] Theodoulidis | Developments in calculating the transient eddy-current response from a conductive plate[END_REF] ∆V

I (t) = d dt t 0 ∆V H (τ ) I(t -τ ) dτ (37)
with the ∆V I and ∆V H denoting the signals for the arbitrary and step excitation, respectively. Noting that ∆V H is the time derivative of s(t) and carrying out integration by parts (37) becomes

∆V I (t) = - d dt t 0 s(τ ) İ(t -τ ) dτ. (38)
For a discrete signal, the time integral reduces to a simple trapezoidal rule. Furthermore, the derivative can be approxi-mated by a finite difference scheme yielding to the following relation

∆V Ii = - i k=0 s k İi-k + i-1 k=0 s k İi-k-1 = i-1 k=0 s k İi-k-1 -İi-k -s i İ0 . (39)
A case of special interest consists the excitation with a rectangular pulse since this is the most usual excitation signal used in PECT applications. In this case the EMF calculation can be readily derived noting that

∆V I (t) = d dt [s(t) -s(t -T )] (40)
where T stands for the pulse duration. Assuming that T = K∆t, the discrete form of (40) admits the particularly simple expression

∆V Ii = (s i -s i-1 -s i-K + s i-K-1 ) /∆t. (41)
We need now to calculate the surface integral in terms of the development coefficients C

i;m and D (s) i;m . The analysis is strictly identical with the harmonic case, and the result reads for the ith timestep

s i = - 2πρ 2 L µ 0 ∞ m=0 κ m J 1 (κ m ρ L ) C (1) i
;m D (s) m . (42)
where the timestep index i has been intentionally suppressed since the source field is constant for t > 0.

The absolute signal at the pickup coil can be easily obtained by recalling that for two coils in air

V 0 (t) = -M dI dt ( 43
)
where M is the mutual impedance between driver and pickup coil. The total signal will be the superposition of ∆V and the signal in air, i.e. V (t) = ∆V (t) -M İ(t).

V. NUMERICAL PERFORMANCE

A number of tuning parameters needs to be fixed in order to proceed to numerical computations. Rules of thump for the estimation of the domain truncation limit ρ L and the maximum number of radial modes N m taken into account for the solution have been proposed in a series of previous works (e.g. [START_REF] Theodoulidis | Eddy Current Canonical Problems (with applications to nondestructive evaluation)[END_REF], [START_REF] Theodoulidis | Impedance of an induction coil at the opening of a borehole in a conductor[END_REF], [START_REF] Skarlatos | Solution to the eddy-current induction problem in a conducting half-space with a vertical cylindrical borehole[END_REF]). The choice of ρ L is principally dictated by geometrical criteria such as the coil external radii. For the series truncation, general rules are more hard to find yet numerical experimentation proves that a number between 100 and 150 modes would be sufficient for the most of the cases.

As far as the number of axial modes q l is concerned, there is again no general rule than can be followed. However, the small thickness of the specimens considered in PECT applications indicate that truncating the Fourier series to no more than 15-20 terms, would hardly affect the accuracy of the solution. The arguments resemble the discussion given in [START_REF] Skarlatos | A modal approach for the solution of the non-linear induction problem in ferromagnetic media[END_REF], [START_REF] Skarlatos | Solution to the eddy-current induction problem in a conducting half-space with a vertical cylindrical borehole[END_REF] concerning the series truncation for the TM potential (W b second order vector potential) in planar specimens with a borehole. Final judge for a safe estimation of these parameters should be the numerical experimentation.

An important speed-up of the time-stepping approach with respect to spectral decomposition using Fourier or Laplace transform, consists in the fact that the sampling frequency appearing in the Helmholtz equation [START_REF]Development of a very fast simulator for pulsed eddy current testing signals of local wall thinning[END_REF] is the same for all iterations, meaning that the system can be factorised using LU decomposition once thus significantly accelerating the overall solution. It is noted that the equation system ( 27),( 29) is diagonal along the φ direction, a fact that facilitates the factorisation. A similar decomposition has been exploited also in a hybrid FIT-modal scheme leading to important speedup factors [START_REF] Skarlatos | A mixed spatial-spectral eddy-current formulation for pieces with one symmetry axis[END_REF]. As a figure of merit, all the numerical computations carried out in the context of this work were of the order of a few seconds for a complete excitation cycle of 100 timesteps, using a modern PC (Intel(R) Core(TM) i7-8850H CPU @ 2.60 GHz).

VI. RESULTS

The proposed approach has been applied for the calculation of the PECT response in the inspection scenario described by Fig. 1. The plate conductivity and relative magnetic permeability are σ = 5 MS/m and µ r = 150, respectively, which correspond more or less the parameters of a typical construction steel. The plate thickness is 1 mm. The inner and outer radius of the driver coil is r 1d = 1 mm and r 2d = 2.65 mm, respectively, its length is z 2d -z 1d = 2 mm and its number of turns is N d 336. The corresponding dimensions for the pickup coil are r 1p = 2 mm and r 2d = 5 mm and z 2p -z 1p = 2 mm, and it is wound with N p =700 turns.

The first computation concerns the magnetic field calculation at a representative point in underneath the plate (ρ o = 0 mm, z o = -3d/2). Two waveforms for the input current have been considered: a rectangular pulse and a pulse with exponential relaxation. The latter case is representative of a realistic current signal provided by a square voltage generator. The pulse duration is taken of 1 ms duration with a 40% duty cycle (dc) for both waveforms, and the pulse amplitude reaches 1 A. The time constant for the exponential pulse is 0.05 ms. The signal waveform is depicted in Fig. 3 The field transients for both inputs are shown in Fig. 4. As reference solution are taken numerical simulations using the FIT method. The FIT solution is also calculated in TD using a timestepping scheme and the difference with the herein presented approach consists in the discretisation of the spatial part using a computational grid [START_REF] Clemens | Transient eddy-current calculation with the FI-method[END_REF], [START_REF] Skarlatos | Start vector generation for implicit Newmark time integration of the wave equation[END_REF].

The induced EMF signals computed with the semianalytical approach for the two signals of Fig. 3 are compared against the corresponding FIT results in Fig. 5. Since the rectangular pulse is a highly idealised situation, the exponential curve described by the second signal is taken as a characteristic current profile when the coil is fed by a typical voltage source [START_REF] Bowler | Pulsed eddy-current response to a conducting half-space[END_REF]. Note that the EMF response is calculated with FIT via direct integration of the magnetic vector potential where the integration is carried out on the cross-section of the pickup coil, with S p = (z 2p -z 1p )(r 2p -r 1p ) and N p /S p standing for the windings density.

V = -2π N p S p d dt pickup A(r) dS (44)
A final comparison is realised between the semi-analytical results for the EMF in the pickup coil calculated using the timestepping approach and the TD reciprocity theorem and the corresponding results obtained using the Laplace approach described in [START_REF] Theodoulidis | Efficient calculation of transient eddy current response from multilayer cylindrical conductive media[END_REF], where the Stehfest algorithm is used for the calculation of the inverse Laplace transform. The calculations have been carried out for step excitation since the Stehfest algorithm does not work for non-monotonous excitations. The comparison is shown in Fig. 6 All results show a very good agreement verifying the accuracy of the TD approach.

VII. CONCLUSION

In this work the TREE semi-analytical approach has been combined with time-stepping integration for the calculation of the eddy-current transient response under arbitrary signals. The main advantage of the time-stepping approach with respect to Fourier or Laplace inversion is the robustness of the solution (Gibbs effect or instabilities due to non-monotonicity of the input). To the drawbacks of the approach can be counted the mathematicam complexity of the scheme (one has to deal with the inhomogeneous Helmholtz equation) and the potentially large number of time samples one has to consider for certain excitations. Laplace-based approaches are much more flexible, providing the possibility of a highly inhomogeneous sampling of the time-axis. Hence, the herein presented approach can be considered as complementary to the Laplace solutions for situations where smooth transients in a reasonable time interval need to be calculated.

The TD version of the reciprocity theorem provides an elegant way of calculation the receiver EMF. Note that the approach presented in section IV is compatible with Laplace calculations of the step response, making the algorithm also applicable beyond the scope of the time-stepping solution of this work.

APPENDIX

In this appendix, we shall derive the TD version of the reciprocity theorem for the special case of the eddy-current problem. The derivation will follow in general Felsen and Marcuvitz [START_REF] Felsen | Radiation and Scattering of Waves[END_REF]. Further references can be sought mostly in the wave propagation literature such as [START_REF] Welch | Reciprocity theorems for electromagnetic fields whose time dependence is arbitrary[END_REF]- [START_REF] Shlivinski | Time-domain coupled response of antennas -reciprocity theorem approach[END_REF]. The list does not claim exhaustiveness.

We start the analysis by writing the Maxwell equations in the free space for an arbitrary current source vector J ∇ × E s = -dB s dt (45)

∇ × H s = J (46)
where the index s stands for source. Introducing the conductive medium, and subtracting the free-space solution from the one in the presence of the conductor, we obtain for the field difference in the source region (the air region where source

Fig. 1 :

 1 Fig. 1: Infinite conducting plate inspected via a pair of coaxial cylindrical coils connected in driver(d) -pickup (p) mode. The geometry consists a three-layer medium.

Fig. 3 :

 3 Fig. 3: Excitation signals used for the solution validation: 40% dc square pulse and pulse with exponential relaxation (chargedischarge mode).

Fig. 4 :

 4 Fig. 4: Semi-analytical vs. numerical (FIT) solution for (a) B ρ and (b) B z field components as a function of time at the selected observation point under the specimen. The transient signals are shown for both current excitations of Fig. 3.

Strictly speaking, the field is not stationary but its value follows the excitation variations instantaneously, that is without delay due to diffusion. In this sense the problem formulation is isomorphic with the magnetostatic equations and hence it will be considered as such.

with E ec = E -E s and similarly for the magnetic field and magnetic induction variations, i.e. H ec and B ec , respectively. We apply the transformation t → t -t, which performs a temporal reflection and a shift of the time axis, to (47),(48) resulting to the adjoint set of equations

with the adjoint electric and magnetic fields E † = E † (r, t) and H † = H † (r, t), which are linked with the original solution of the untransformed equation E and H, via the relations

The magnetic induction B † is defined accordingly. Following the usual procedure, we multiply (45),( 46) with H † ec and E † ec and (47),(48) with H s and E s , respectively, use of the vector identity

one gets after some trivial manipulations

and multiplying with the permeability of the free-space µ 0

Integrating (54) throughout the air region and applying the Gauss theorem in order to get rid of the divergence operator, we obtain

where dS stands for the outwards pointing normal to the boundary ∂V vector.

Using (1),( 2), the integrand of the boundary integral can be written successively

and recalling that in the 2D problem the magnetic potential is lying along the azimuthal direction A = Ae φ , the last relation reduces to

Integrating (55) in time, taken into account the previous relation for the boundary integral and rearranging the terms we arrive at the equation

Making the assumption that the the system initially is in the zero state and that the excitation current starts from zero value the first term of the right-hand-side vanishes since A s (r, 0) = 0 and A † ec (r, t ) = A ec (-r, 0) = 0. Furthermore, the left-hand-side is equal to the time convolution of the EMF at the pick-up coil with the excitation current ∆V (t ) * I(t ), yielding the relation

with (• * •) standing for the time convolution integral. Equation (57) provides a very elegant way of calculating the EMF of the receiving coil, in direct analogy with the classical Auld's reciprocity formula [START_REF] Auld | Theoretical characterization and comparison of resonantprobe microwave eddy-current testing with conventional low-frequency eddy-current methods[END_REF], [START_REF] Auld | Eddy current probe response to open and closed surface flaws[END_REF] for the impedance variation calculations in harmonic problems.


Michel Darmon

Adrien Ferrand

Vincent Dorval

Sylvain Chatillon

Sebastien Lonne

Sébastien Lonné

Recent Modelling Advances for Ultrasonic TOFD Inspections

Keywords: TOFD simulation, edge diffraction, lateral wave, SDH PACS: 43.20.Gp -43.35.Yb -43.35.Cg -43.35.Pt -43.35.Zc

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

TOFD (Time Of Flight Diffraction) is a widespread NDT ultrasonic technique used for locating and sizing cracks by their diffraction echoes. This technique relies on an arrangement of two probes of opposite beam directions. Figure 1 shows some typical phenomena arising for such an arrangement over a planar specimen containing an embedded crack-like flaw:

• Diffraction echoes (2-2') are scattered by the edges of the flaw, and caught by the receiving probe • A lateral wave (1) is transmitted, which follows the specimen profile • A backwall echo (3) is observed, which corresponds to the specular reflection of the beam over the outer wall of the component.

From the B-Scan image of Figure 1 one clearly distinguishes the top and bottom diffraction echoes of parabolic shapes, as well as the lateral and backwall echoes. Those echoes are conventionally used to size the defect: from the knowledge of backwall and lateral waves, beam propagation and wave speed in the material, one may straightforward estimate the locations of the edges of the flaw.

FIGURE 1. Illustration of TOFD inspection technique over a planar specimen and related echoes for an embedded crack

Depending on the size of the flaw and its position, shadowing effects over the backwall echo may arise, as the beam will be partly reflected by the flaw and lead to a lower backwall echo. This is notably the case for the configuration of a backwall breaking flaw depicted in Figure 2.

FIGURE 2. Illustration of TOFD inspection technique over a planar specimen and related echoes for a backwall breaking crack

Simulation is helpful for evaluating the performances of inspection techniques [1]. Several developments have been carried out in the past to provide, TOFD modelling in the CIVA software. Simulation of TOFD inspection henceforth involves the necessary simulation of complex phenomena due to the 3D geometries of flaws and structures. Consequently several improvements have been newly integrated in the TOFD simulation tools.

LATERAL WAVES

So as to simulate lateral waves, two different developments have been devised. The first one is devoted to planar components: it was previously added to CIVA and recently improved. The second one is newly developed and dedicated to lateral waves on irregular surfaces; the validation of its complete process has to be successfully finished.

Lateral Waves on Planar Components

Concerning the simulation of lateral waves on planar components, the main recent improvement has consisted in the replacement of the previous 2D model [START_REF] Cerveny | Theory of seismic head waves[END_REF] by a 3D one. Such a simulation tool follows a procedure depicted in Figure 3: 1) Discretization of transmitting and receiving surfaces into elementary surfaces surrounding mesh points 𝑃 𝑖 and 𝑃 𝑗

2) Between two points 𝑃 𝑖 and 𝑃 𝑗 , respectively located on the emitting and receiving surfaces: Calculating the 3D path of the lateral head wave (HW). Application of asymptotic ray theory (Cerveny's model [START_REF] Cerveny | Theory of seismic head waves[END_REF]) for calculating each ray amplitude. 3) Frequency loop and summation over the transmitter and receiver surfaces. The lateral wave is also called head wave since it is the first arrival wave received during a TOFD (Time Of Flight Diffraction) inspection. For a planar surface, a head wave is generated at critical incidence at the interface between the material and the emitter, propagates in the material along the surface, and is radiated at critical angle as a bulk wave in the wedge, reaching the receiver. For planar surfaces, the head wave is seen as the result of a critical refraction phenomenon of the incident wave on the interface between the two media. Previously the head wave amplitude was computed by summing contributions from only 2D paths aligned along the axis connecting the two probes centers [START_REF] Cerveny | Theory of seismic head waves[END_REF]. In the framework of the 3D new model for lateral waves on planar components, the paths of the lateral wave are calculated in 3D; consequently the contribution of some new rays skewed with respect the probe centers axis is truly taken into account. Numerical and experimental validations are in progress.

Shadowing by a near surface crack

In TOFD inspections, a surface breaking flaw prevents the lateral wave to be transmitted from one probe to the other (Figure 4 left). That's why head waves can sometimes be used for surface breaking cracks detection. On the other hand, when the crack is located near the surface and separated from it by a distance called ligament, complex interferences between the head wave and the top edge diffracted wave occur. There is not anymore interruption of the lateral wave but there is still an effect of the ligament on the lateral wave amplitude.

To model this ligament effect, an empirical model which has been integrated in CIVA11 release has been developed to simulate, in a simplified manner, the interferences between the head wave and the top edge diffraction. An amplitude law for the head wave amplitude has been adjusted using CIVA/ATHENA [START_REF] Mahaut | Modeling of Corner Echo in UT Inspection Combining Bulk and Head Waves Effect[END_REF] finite elements calculations as reference:

𝐴 𝐻𝑊 = 𝐴 𝐻𝑊0 ( 𝑙 𝑙𝑖𝑔𝑎𝑚𝑒𝑛𝑡 𝑙 𝑐 ) 2
where 𝐴 𝐻𝑊0 is the simulated lateral head wave without flaw, 𝑙 𝑙𝑖𝑔𝑎𝑚𝑒𝑛𝑡 is the ligament and 𝑙 𝑐 = 3√𝜆 where λ is the wavelength.

Thanks to this empirical adjustment of the amplitude law, this simplified model leads for several metallic materials to results similar to finite elements ones except for ligaments close to the critical ligament 𝑙 𝑐 .

Examples of simulations of the lateral wave using this empirical model are given hereafter. Analytical Civa simulations have been performed for TOFD inspections on a plane specimen using two identical probes with a 70mm Probe Center Space (PCS) emitting P60° waves at 2MHz. The specimen contains a flaw of fixed bottom edge location but varying ligament. For 6 mm ligament, a total transmission of the lateral wave is observed and there is no flaw influence on the head wave. When decreasing the ligament, a decrease of the HW amplitude is predicted, the HW becoming less predominant than the bottom edge diffraction echo. When the flaw is a surface breaking one, a total interruption of the lateral wave is modelled and the bottom edge diffraction echo is the only observed one.

The simulation provides thus a progressive extinction of the lateral wave with a decreasing ligament.

Lateral Waves on Irregular Surfaces

TOFD and other pitch-catch configurations are also classically used on non-planar interfaces in NDE (inspections of cylinders, nozzles [3]…).

The previously developed model for head waves simulation is based on the asymptotic ray theory [START_REF] Cerveny | Theory of seismic head waves[END_REF] which has also been used to improve the modeling of flaw corner echoes [START_REF] Mahaut | Modeling of Corner Echo in UT Inspection Combining Bulk and Head Waves Effect[END_REF]. Nevertheless, this existing head wave simulation can only be applied for a planar specimen composed of an isotropic medium.

In [START_REF] Ferrand | [END_REF], complete interpretation and simulation of the head wave propagation mechanisms for irregular surfaces are proposed. This approach leads to the development of a generic algorithm of ray tracing between interface points (called Generic Interface Ray Tracing/GIRT) which is able to compute the travel path of head waves in specimens inspected by TOFD technique. Based on an adaptation to the NDE domain of seismic ray tracing algorithms, GIRT is a new method for solving the two-points ray tracing problem, which is the calculation of the true ray between two defined points: the source and observation ones.

The GIRT algorithm is adapted to find not only the head wave ray path (first arrival) but also all complex later arrival waves including mode or nature conversion by adding constraints on the searched ray path. Secondly, GIRT can take into account all types and natures of waves propagation occurring in NDE. Indeed the propagation simulation is available for both P and SV waves and accounts for modes conversions. Surface waves as Rayleigh waves and head waves are taken into account as well as diffractions from surfaces in any kind of waves. Finally, GIRT approach is extended to flaw scattering modeling since the flaws contained in the specimen are also meshed by the GIRT algorithm in order to model rays diffracted from flaws. Consequently the GIRT algorithm is generic for application to NDE as it can deal with any ultrasonic wave propagating near irregular surfaces or flaws of any geometry (CAD defined).

By modeling the complex wave interactions with the surfaces thanks to this ray approach, GIRT is thus able to correctly predict the head wave time of flight and to interpret the complex waves propagation at the vicinity of an irregular interface. Indeed, comparisons between GIRT and FEM simulations of the head wave fronts near the receiver surface have been previously presented in [START_REF] Ferrand | [END_REF]. A good agreement is obtained between GIRT and FEM results, which validates the GIRT algorithm. For very corrugated surfaces, a bulk propagation mechanism has to be taken into account for the head wave modeling. GIRT has been integrated in CIVA11 release to predict the valid time of flight for lateral waves on irregular surfaces.

A 2D model for simulating precisely the lateral waves amplitude on irregular surfaces has then been devised; it has not been yet integrated in CIVA, its validation being in progress. For the moment, this model can deal with cylindrical (Figure 6) or scouring irregularities (Figure 7). The computation of head wave amplitude is derived from ray-based models.

The head wave ray path is determined by GIRT between source and observation points located over the probes surfaces. The ray path is divided into elementary paths and the application of an amplitude ray model for each interaction corresponding to each elementary path is carried out. Two specific new ray models have been developed: one to simulate the creeping ray along a cylindrical part and one to model a grazing ray along a planar part (Figure 7). Indeed, these two elementary paths are involved in the head wave path on scoured surfaces.

The developed models and their validations will be described in a future publication. Encouraging results have been obtained for simulating 2D lateral waves on cylindrical or scoured surfaces.

CRACK DIFFRACTION

In order to predict TOFD inspection performance in realistic configurations, the simulation needs to account for arbitrarily oriented flaws with respect to the probes orientation. In addition to tilt angle (flaws not necessarily perpendicular to a planar backwall), skew or vertical disorientation angles may apply, for instance using inspection over pipes or nozzles. Moreover, realistic flaws, may have a complex contour shape (not necessarily rectangular aperture), so developments were also carried out to account for CAD contour planar flaws and in addition branched or multifaceted flaws. Edge diffraction echoes can be modelled identically by GTD or new PTD [6] [7] models. Two different strategies have been investigated concerning the calculation of the GTD coefficients. A so-called "projected 2D" option, which is based onto the projection of the incoming and scattered wave vectors over the plane which contains the normal of the flaw edge, was previously available. It has been removed and replaced in CIVA11 release by a "pure" 3D GTD code, (metre une reference). The advantages of the 3D codes are a better prediction for strong 3D crack misorientation and the use of a new TV-TH decomposition: continuous limit when turning the 3D problem into 2D one. An example of experimental validation of the GTD model implemented in CIVA is performed hereafter. The reference block, of 30mm thickness, is composed of steel and incorporates 11 artificial notches of extension 15mm and heights varying from 0.5mm to 15mm (Figure 8). The material homogeneity was experimentally checked and the specimen material (steel) was modelled as isotropic, homogeneous and the attenuation was ignored. The longitudinal (L) and transversal (T) velocities were estimated using successive backwall echoes (VL = 5900m.s -1 and VT = 3230m.s -1 ), the density was 7.8g.cm -3 .

The two probes arranged in TOFD mode were moved in two directions over the notches with a scanning step of 0.2 mm and an increment step of 1mm. The probes used were single element conventional round probes with a centre frequency of 5MHz and a diameter of 6.35mm. Concerning the material properties of the wedges, the longitudinal and transversal velocities are respectively VL = 2730m.s -1 and VT = 1340m.s -1 and the density is 1.18g.cm -3 .

TOFD inspection of the back-wall breaking notches was performed for probe centre space "PCS" (Figure 9) varying from 35mm to 100mm, allowing various incidence angles of the L waves on the top edge and various positions of the "L axis crossing point" (point where the L axis of the transmitter and receiver intersect) relative to that of the top edge (Figure 9). Figure 10 shows the experimental and simulated amplitudes of the top edge echo obtained for the TOFD inspection of the 15mm height notch performed with the probes generating longitudinal waves at 60° in the specimen and for various PCS. Three different simulations are shown: one using GTD model and two simulations using the finite elements model CIVA/ATHENA [START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for NDT[END_REF] both for a null or 0.2mm notch aperture. The reference for the amplitude is the L direct echo of the SDH Ø2mm at 20mm depth obtained in TOFD mode with a PCS of 70mm.

When the PCS increases, the "L axis crossing point", located above the top edge for the PCS of 35mm, moves deeper, reaches the top edge for PCS = 50 mm and then drops below it. The amplitude of the top edge echo reaches its maximum when the crossing point is close to the top edge.

We have a good agreement between experiment and GTD or CIVA-ATHENA predictions whatever the notch aperture. In this case, we observe almost no effect of the notch aperture. Indeed we can see on the right of Figure 10 that, whatever the PCS, the incident angle θ is less than 130°, and influence of the notch aperture on the top edge echo amplitude occurs only for θ angles greater than 130°. Differences between GTD and CIVA-ATHENA appear for such angles owing to an additional specular echo coming from the horizontal part of the notch top. Contrary to CIVA-ATHENA, the GTD model does not compute this specular echo.

The same good agreement between experiments and the three simulations (with GTD and CIVA-ATHENA with quasi null or 0.2mm aperture) was obtained for the notches of 10mm and 5mm.

CALIBRATION DEFECTS

Until CIVA10 release, the response of a side drilled hole (SDH) which is a classical calibration flaw was only modelled using the Kirchhoff approximation: it takes into account specular reflection on a cavity but not the creeping waves which propagate around the cavity circumference (Figure 13). Therefore, an exact analytical solution for the scattering from a cylindrical cavity, based on the Separation Of Variables (SOV) method [START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for NDT[END_REF], has been implemented in CIVA 10 to improve simulation of the SDH response [START_REF] Darmon | [END_REF].

For validation purposes, TOFD inspections of a planar block containing SDHs were performed by varying the PCS of 50mm to 80mm. When the PCS increases, L theoretical axes of the two sensors intersect above and below the TG of different diameters located at 20mm deep. The crossing point is at the SDH depth when the PCS is about 70mm. These axes and the angles of incidence on the TG are shown in Figure 11. The SOV model predicts the existence of the two SDH echoes (specular reflection and creeping wave) but their relative amplitudes are not always perfectly correct: in the case of small PCS, 50mm or 60mm) for which the intersection point of axes is above the SDH, the second echo is overestimated by CIVA. When the depth of the intersection point is on the SDH the amplitudes of two echoes (which are similar) are very well predicted. Nevertheless, the maximal amplitudes obtained with CIVA10 for this Ø2mm SDH versus the PCS are always in good agreement with those measured even for PCS 50mm or 60mm when the amplitude of the second echo is not correct (Figure 12) because its amplitude, smaller than the first echo, does not really affect the measurement of the maximum amplitude of the SDH echo. For more important PCS, there is more field amplitude on the SDH bottom; moreover, when increasing the PCS as shown in Figure 13, the path for creeping waves is shorter and their attenuation is smaller. Consequently the creeping waves amplitude increases with the PCS.

Improvement due to creeping waves modelling (not accounted in Kirchhoff model) is henceforth obtained using the CIVA SOV model.

CONCLUSION

Recent developments have been carried out to enhance the abilities of ultrasonic TOFD inspection simulation. New SOV exact analytical model enables to improve simulation of a typical calibration flaw, a SDH, by accounting for creeping waves. Lateral waves on planar components are henceforth simulated using a 3D model and accounting their shadowing by a near surface crack. A 2D model of lateral waves on irregular surfaces has also been devised for cylindrical or scoured irregularities. Intensive experimental validations of all these improvements are in progress. Simulation of lateral waves on irregular surfaces could be further improved to deal with 3D configurations, wedgelike irregularities, impedance interfaces and anisotropic media. So as to simulate edge diffraction, integration of 3D GTD model and of the generic PTD has been performed leading to successful validations. Works are in progress to improve edge diffraction simulation near critical angles.

FIGURE 3 .

 3 FIGURE 3. Procedure for simulating lateral waves on planar components using contact probes.

FIGURE 4 .

 4 FIGURE 4. TOFD inspections on a planar component: at left, case of a surface breaking crack; at right, case of a near surface crack.

FIGURE 5 .

 5 FIGURE 5. CIVA Simulations of TOFD inspections on a planar component containing a flaw of fixed bottom edge location but varying ligament (from 6 to 0 mm).

FIGURE 6 .

 6 FIGURE 6. Cylindrical irregularity of radius a.

FIGURE 7

 7 FIGURE 7 Amplitude calculation along a ray path calculated by GIRT over a scoured irregularity. Red dashed blocks: creeping ray along cylindrical parts. Blue dashed block: grazing ray along the planar part.

FIGURE 8 :

 8 FIGURE 8: Steel specimen, 30 mm thickness, with back-wall breaking artificial notches of height 0.5mm, 1mm, 1.5mm, 2mm, 3mm, 4mm, 5mm, 7.5mm, 10mm, 12.5mm and 15mm (15 mm extension).

FIGURE 9 :

 9 FIGURE 9: CIVA schematics of TOFD configurations with different PCSs, showing the variation of the L waves incident angle on the top edge of the notch (θ) and of the L axis crossing point position relative to that of the top edge with the PCS.

FIGURE 10 :

 10 FIGURE 10: L60° TOFD inspection, notch of 15mm height. Left) Compared experimental and simulated amplitudes of the L top edge echo versus the PCS. Three different simulations shown: one using CIVA GTD model and two simulations using the finite elements model CIVA/ATHENA both for a null or 0.2mm notch aperture. Right) Inspection configurations for different PCS showing the relative positions of the notch and the probes.

FIGURE 11 :

 11 FIGURE 11: SDH at 20mm depth, ray path display and value of the direct ray incidence angle on the SDH. TOFD inspection, contact probe Ø6.35mm, 5MHz, P60°, various PCS. Experimental Ascans, measured on a Ø2mm SDH located at 20mm depth are compared to simulated ones, using SOV model, for different PCS.The SOV model predicts the existence of the two SDH echoes (specular reflection and creeping wave) but their relative amplitudes are not always perfectly correct: in the case of small PCS, 50mm or 60mm) for which the intersection point of axes is above the SDH, the second echo is overestimated by CIVA. When the depth of the intersection point is on the SDH the amplitudes of two echoes (which are similar) are very well predicted. Nevertheless, the maximal amplitudes obtained with CIVA10 for this Ø2mm SDH versus the PCS are always in good agreement with those measured even for PCS 50mm or 60mm when the amplitude of the second echo is not correct (Figure12) because its amplitude, smaller than the first echo, does not really affect the measurement of the maximum amplitude of the SDH echo. For more important PCS, there is more field amplitude on the SDH bottom; moreover, when increasing the PCS as shown in Figure13, the path for creeping waves is shorter and their attenuation is smaller. Consequently the creeping waves amplitude increases with the PCS.Improvement due to creeping waves modelling (not accounted in Kirchhoff model) is henceforth obtained using the CIVA SOV model.

FIGURE 12 :

 12 FIGURE 12: Comparison of measured (black) and simulated with CIVA Ascans from a 2mm diameter SDH for different diameters PCS. TOFD inspection with Ø6.35mm probe, 5MHz, P60°.

FIGURE 13 :

 13 FIGURE 13: Specular path (in green) and creeping wave path (in blue) around the side drilled hole which is shorter for great PCS than for the little corresponding to a weaker energy loss by irradiation along this path.


Vincent Dorval
email: vincent.dorval@cea.fr

Nicolas Leymarie

Sylvain Chatillon

Improved Semi-Analytical Simulation of UT Inspections Using a Ray-Based Decomposition of the Incident Fields

and research institutions in France or abroad, or from public or private research centers.

Improved semi-analytical simulation of UT inspections using a ray-based decomposition of the incident fields

INTRODUCTION

The simulation of ultrasonic Non Destructive Evaluation has a wide range of industrial applications. It can, for example, help the design of inspection procedures or the interpretation of measurements. Simulation methods can rely either on semi-analytical methods, numerical methods (such as finite elements), or on a combination of both. Semi-analytical methods apply approximations to perform faster computations than numerical ones, which can come at the cost of a loss of accuracy in their results.

Plane-wave scattering coefficients are used by many semi-analytical methods, including those described by Thompson and Gray [1] and by Schmerr and Song [START_REF] Schmerr | Ultrasonic Measurement Models[END_REF], as well as those implemented in the Civa software [START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for in Ultrasonic Wave Propagation in Non Homogeneous Media[END_REF]. Planewave and quasi plane-wave approximations of the ultrasonic fields provide the input values for the calculations of these coefficients. The quasi plane-wave approximation allows accounting for some effects of the probe aperture in the spectral content of the echo. This can be done through analytically calculated field coefficients for simple geometries [1]. The Civa approach is meant to be generic and able to handle complex geometries, therefore its field calculations rely on a semi-analytical paraxial ray method. A strict plane wave approximation is applied to these fields during the echo computation and some information is lost about aperture and spectral content. Such an approximation can be more or less accurate, depending on the situation. It typically yields satisfactory results when applied in far field and in the focal axis of the probe, but can produce significant errors when these conditions are not met. The simulation method presented in this communication is being developed with the aim of correcting these errors.

EXPRESSION OF ULTRASONIC ECHOES

Ultrasonic echoes can be expressed using Auld's reciprocity theorem [START_REF] Auld | [END_REF]. With some notation modifications and the introduction of the amplitude Sinput of the input signal, their expression can be adapted into:

𝑆 𝑒𝑐ℎ𝑜 = 𝑆 𝑖𝑛𝑝𝑢𝑡 ∫ (𝒖 𝐸𝑚 • 𝑻 𝑅𝑒𝐹 -𝒖 𝑅𝑒𝐹 • 𝑻𝒏 𝐸𝑚 )𝑑Σ Σ 𝐹 . ( 1
)
Sinput is related to the electric signal that enters the acquisition chain and also accounts for all the filtering effects not related to ultrasonic phenomena. It is considered as an input parameter of the model. u and T refer respectively to the displacement vectors of ultrasonic fields and to their stress tensors projected on the normal of the defect. The indexes Em and ReF correspond to two configurations: Em is the ultrasonic field radiated by the emitting probe when the flaw is absent. ReF is the ultrasonic field radiated by the receiving probe when it is acting as an emitter and when the flaw is present. ΣF is the surface of the flaw. Secho, Sinput and ultrasonic fields all depend on the frequency.

Expression ( 1) is general and independent of any approximations. In practice, its evaluation generally requires various degrees of approximation, depending on the approach. In many approaches the expression is approximated in ways that allows replacing the integrand by plane-wave interaction coefficients [1][START_REF] Schmerr | Ultrasonic Measurement Models[END_REF][START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for in Ultrasonic Wave Propagation in Non Homogeneous Media[END_REF]. These coefficients have been obtained analytically for several approximations and types of flaws. Several of them are listed by Schmerr and Song [START_REF] Schmerr | Ultrasonic Measurement Models[END_REF].

Based on this, a possible way to compute an echo can be expressed:

𝑆 𝑒𝑐ℎ𝑜 𝑝𝑙𝑎𝑛𝑒_𝑤𝑎𝑣𝑒 =
The integral has been replaced by a summation over the surface of the flaw on elements of surface ΔΣ centered around positions x. The integrand has been replaced by the product of plane wave amplitudes A and of a scattering coefficient CF. The interaction with the flaw is handled by CF, the amplitudes A of the fields are in the absence of defect. The directions and polarizations d and p correspond to approximations of the field as plane waves at the location x.

In this formalism, the amplitudes A also contain phase information concerning the incident and received waves. In the strictest sense of the plane-wave approximation, the phase varies with angular frequency ω as 𝑒𝑥𝑝(𝑖𝜔𝜏) (or its complex conjugate, depending on conventions), expressing only a time-of-flight τ. In the less quasi plane waveapproximation, the frequency dependence of A can account for some effects related to probe aperture.

The semi-analytical echo computation of the Civa software relies on a strict plane-wave approximation, where the waves at a location are described by only their amplitude and time-of-flight. This simplification allows a gain in computation time: both these quantities can be interpolated, which allows doing the costly field computations on a coarse grid, extracting quantities, and then interpolating them on a finer grid. However, it comes at a cost in accuracy.

RAY BASED APPROACH

A decomposition of the field as a sum of rays is readily available in the Civa software, as the ultrasonic field computations of the Civa software rely on them. They can handle complex probe geometry and arbitrary numbers of interfaces with various shapes. To perform echo computations in the plane-wave approximation, the contributions of the rays are summed and the resulting field is approximated. A straightforward way to improve the accuracy of that approach would be to treat each ray as a different plane wave instead of approximating their sum as one. It could be written as follows: (

𝑆 𝑒𝑐ℎ𝑜 𝑟𝑎𝑦𝑠 = 𝑆 𝑖𝑛𝑝𝑢𝑡 ∑ ∑ ∑ 𝐴 𝑟𝑎𝑦
) 3
The sums up to NEm x and NRe x are over the ensembles of rays found between the current point x of the defect and the surface of the emitting and receiving probes. This formulations implies that NEm•NRe coefficients CF need to be computed for each point x. This can induce significant increase in computation times compared to equation ( 2) that only required one coefficient for each point. Keeping computation times low is an important concern here and specific adjustments have been made to this end.

In general, parts of the computation of the interaction coefficient can be performed independently of emission, reception, or both. This allows limiting the operations that are performed NEm•NRe times. Different ways to make the computation faster can be available, depending on the model considered. Examples are given below for the Kirchhoff model for cracks and for an exact model for cylindrical defects. They give a general idea of various ways to perform fast computation. Similar developments have been made for other models, such as the GTD and Born models. Additionally to that, the implementation of the model uses parallel computing over the points x of the defect.

Example: the Kirchhoff Approximation for Cracks

The Kirchhoff approximation is also known as the tangent-plane approximation: the field at a point of the surface of a defect is assumed to be equal to the sum of the incident wave and of the waves that would be reflected on a tangent infinite plane. Chapman [START_REF] Chapman | Ultrasonic scattering from smooth flat cracks: an elastodynamic Kirchhoff diffraction theory[END_REF] gives plane-wave coefficients for the entire surface of defects of canonical shapes. It's also possible to perform a discrete sum over the surface, following the forms of equation ( 2) or [START_REF] Darmon | Modelling of scattering of ultrasounds by flaws for in Ultrasonic Wave Propagation in Non Homogeneous Media[END_REF], in order to deal with defects of any shape. This approach has also the advantage of for variations of the field at the surface, and it is followed by the Civa computation. In the plane-wave approximation, the implementation corresponds to the following equation, whose terms can be identified to those of both ( 1) and ( 2):

𝑆 𝑒𝑐ℎ𝑜 𝑝𝑙𝑎𝑛𝑒_𝑤𝑎𝑣𝑒 = -𝑆 𝑖𝑛𝑝𝑢𝑡 ∑ 𝒖 𝑅𝑒𝐹 𝑃𝑙𝑎𝑛𝑒𝑊𝑎𝑣𝑒𝐴𝑛𝑑𝐼𝑛𝑡𝑒𝑟𝑎𝑐𝑡𝑖𝑜𝑛_𝒙 • 𝑇 𝐸𝑚 𝑃𝑙𝑎𝑛𝑒𝑊𝑎𝑣𝑒_𝒙 𝒙∈𝑃𝑜𝑖𝑛𝑡𝑠 ∆Σ. (4)
The term that includes stress in presence of the flaw that appears in (1) vanishes here as the crack is considered as a free surface with no normal stress.

The proposed ray-based version of (4) can be written:

𝑆 𝑒 𝑝𝑙𝑎𝑛𝑒_𝑤𝑎𝑣𝑒 = -𝑆 𝑖𝑛𝑝𝑢𝑡 ∑ (∑ 𝒖 𝑟𝑎𝑦 𝑅𝑒 𝑗 𝑃𝑙𝑎𝑛𝑒𝑊𝑎𝑣𝑒𝐴𝑛𝑑𝐼𝑛𝑡𝑒𝑟𝑎𝑐𝑡𝑖𝑜𝑛_𝒙 𝑁 𝑅𝑒 𝒙 𝑗=1 ) • (∑ 𝑇 𝑟𝑎𝑦 𝐸𝑚 𝑖 𝑃𝑙𝑎𝑛𝑒𝑊𝑎𝑣𝑒_𝒙 𝑁 𝐸𝑚 𝒙 𝑖=1 ) 𝒙∈𝑃𝑜𝑖𝑛𝑡𝑠 ∆Σ. (5)
The brackets around the sums emphasize something important: the emitted and received part can be computed and summed independently from each other. Therefore, the interaction needs only to be computed NRe times, not NEm•NRe times. This can lead to significant gains in computation times.

In the implementation, the results are given the time domain.

Example: Exact Model for Cylindrical Flaws

The problem of the scattering of plane elastic waves by cylinders can be solved exactly for isotropic materials using a separation of variables. For example, Lopez-Sanchez et al. [START_REF] Lopez-Sanchez | [END_REF] give expressions for voids in two dimensions. In Civa, two dimensions versions of the coefficients are implemented for both voids and inclusions. They are used for computations both in two and three dimensions, the latter case being handled by projections and discretization along the axis of the cylindrical defects.

Due to the separation of variable, the scattering coefficient C can be written in the following general form (in the two dimensions case):

𝐶(𝜔, 𝑅, 𝜃) = ∑ 𝐹 𝑛 (𝜔𝑅)𝐺 𝑛 (𝜃) ∞ 𝑛=0 . ( 6
)
In practice, the infinite sum is truncated for calculations. Fn is a function that does not depend on the scattering angle θ and Gn is a function that does not depend on the angular frequency ω or the flaw radius R. For a given echo calculation, the radius and the range of frequencies of interest are constants. Therefore, the required values of Fn can be precomputed. This accelerates significantly the evaluation of equation ( 3), as the many evaluations of C only require the evaluations of the Gn coefficients and their combination with the corresponding precomputed Fn.

COMPUTATION RESULTS

The ray-based model has been implemented in a development version of the Civa software. Computations were run using three approaches: the existing Civa model based on the plane-wave approximation, the new ray-based model, and the finite element Civa-Athena hybrid approach [START_REF] Mahaut | Simulation of complex ultrasonic NDT cases using coupled analytical-numerical method : the Mohycan Project[END_REF]. The latter relies on the Athena finite element model for phenomena that occur in a small zone surrounding the defect, while the propagation between the probes and that zone is handled by the same paraxial ray model as the other two models. Finite elements are used here as a reference model for the computation of the interaction between waves and defect. As this approach is only available in 2D, computations with all three models have been performed under 2D assumptions to allow for meaningful comparisons.

The inspected defects are cylindrical in all the cases discussed below, and handled by the exact model described above in the two semi-analytical approaches. This model has the advantage of yielding exact interaction coefficients, which removes one possible source of difference between semi-analytical and finite elements model. It helps isolating the effects of the ray-based description of the fields versus the plane-wave one.

The cylindrical defects are voids with a 2mm diameter. They are located in blocks of steel with longitudinal and transverse wave velocity of respectively 5900 m/s and 3230 m/s, and a density of 7.8 g/cm3. The probes have a diameter of 12.7 mm and a center frequency of 2 MHz.

Comparison of the Results

In the first application case, the defect is located at a depth of 15 mm in the focal axis. The probe moves from a distance of 100 mm to 1 mm of the block. The configuration is illustrated in figure 1(a). The signals obtained by the three methods (not here) over all probe positions have similar shapes, but their amplitudes vary and are plotted in figure 1(b). The ray-based method is overall in better agreement with finite elements that the plane wave approximation. The improvement is less than 1 dB in this example but is clearly noticeable on the right side of the plot, which corresponds to short probe-defect distances. A possible explanation for that result is fact that the plane-wave approximation only simulates backscattering in that situation, whereas the ray-based method accounts for the wider range of incident and scattering direction involved. This range is widest when the probe is closest to the defect.

The second case is similar to the first, except that the probe remains at a distance of 10 mm from the block and moves parallel to the surface to a maximum distance of 25mm relative to the vertical of the defect, as illustrated in figure 2(a). shows that all three models predict signal echoes of similar shapes when the defect is along the focal axis of the probe, though the plane-wave result is distinguishable from the others and has a slightly lower amplitude (this is more apparent on figure 1(b), where the result at -10 mm correspond to the same probe position as here).

A clear disagreement appears on figure 2(c) where the defect is no longer located along the focal axis. The planewave model, which relies on only identifying only a unique plane wave associated with a unique time of flight, is unable to account for the wide range of times-of-flight corresponding to the surface of the probe seen from that angle. To the credit of that model, the maximum amplitude of its signal is not too far from amplitudes of the two others: even in this unfavorable case the plane-wave model can be sufficient if the objective of the calculation is only to have an estimation of the maximum amplitude of the echo. The ray-based model, however, remains in good agreement with finite-elements even for signal shapes thanks to its more complete description of the ultrasonic field.

The next case (figure 3) is a Time Of Flight Diffraction (TOFD) configuration where wedges have been used to emit and receive L waves at a 45 degrees angle. The impact points of the probes are separated by 50 mm and the SDH is at a depth of 15 mm. The results are presented for the L wave echoes at two probe positions. In this example, the SDH is out of the focal zone for the two probe positions shown. The ray-based method is therefore expected to provide an improvement compared to the plane-wave approximation in both. It appears that this improvement is more significant in figure 3(c), where the signals are more distorted by probe diffraction, than in figure 3(b). As in the previous case, the plane-wave approximation gives a reasonable approximation of the maximum echo amplitudes but the signal shapes are significantly improved by the ray-based method.

In the example illustrated figure 4(a), the surface of the specimen contains a 3mm drop in order to test the ability of the model to take into account the effect of beam-splitting due to a complex geometry. Using the upper section of the surface as a reference, the probe is located 10mm above and the SDH is at a depth of 15mm. The result of figure 4(b) and (c) shows this case is, as expected, very unfavorable to the plane-wave approximation. The ray-based model, on the other hand, manages to match closely finite element results. For all three models, the propagation through the irregular surface is computed by the same ray method: more complex phenomena may occur in reality for this particular case.

No comparisons with experimental results are presented in this communication, but some have been studied in various configurations. Examples are given in [START_REF] Toullelan | Results of the 2015 UT Modeling Benchmark Obtained with Models Implemented in CIVA[END_REF]. They lead to a similar conclusion as here: the ray-based method tends to be more accurate than the plane-wave approximation.

Computation Times

The previous results show that the new ray-based model allows improving on the plane-wave approximation by being in better agreement with finite elements. However, finite elements are expected to be the most reliable approach of the three in terms of accuracy. The main advantage of semi-analytical methods compared to finite elements, from a user point of view, is their speed. Therefore, the computation times of the new ray-based method have to remain low for it to be viable.

Table 1 summarizes the computation times corresponding to the examples presented in figures 1 to 4. In all of them, the signals have been computed for 101 probe positions. For all these configurations, the ray-based method is slower than the plane-wave approximation, but remains several orders of magnitude faster than finite elements.

CONCLUSION

The ray-based method presented here allows overcoming limitations of the plane-wave approximation currently used in the semi-analytical ultrasonic echo computations of the Civa software. Examples of results illustrate some of the situations where the ray-based method is significantly better than the plane-wave approximation, due to the

FIGURE 1 .

 1 FIGURE 1. (a) Simulated setup for the illustration of near-field effects; (b) Maximum amplitude of the echo from the sidedrilled hole as a function of probe distance.

FIGURE 2 .

 2 FIGURE 2. (a) Probe positions for the illustration of the effect of probe diffraction; (b) Signals simulated when the probe is above the SDH; (c) Signals simulated when the probe is shifted 25mm away.

Figure 2 (

 2 Figure 2(b)shows that all three models predict signal echoes of similar shapes when the defect is along the focal axis of the probe, though the plane-wave result is distinguishable from the others and has a slightly lower amplitude (this is more apparent on figure1(b), where the result at -10 mm correspond to the same probe position as here).A clear disagreement appears on figure 2(c) where the defect is no longer located along the focal axis. The planewave model, which relies on only identifying only a unique plane wave associated with a unique time of flight, is unable to account for the wide range of times-of-flight corresponding to the surface of the probe seen from that angle. To the credit of that model, the maximum amplitude of its signal is not too far from amplitudes of the two others: even in this unfavorable case the plane-wave model can be sufficient if the objective of the calculation is only to have an estimation of the maximum amplitude of the echo. The ray-based model, however, remains in good agreement with finite-elements even for signal shapes thanks to its more complete description of the ultrasonic field.The next case (figure3) is a Time Of Flight Diffraction (TOFD) configuration where wedges have been used to emit and receive L waves at a 45 degrees angle. The impact points of the probes are separated by 50 mm and the SDH is at a depth of 15 mm. The results are presented for the L wave echoes at two probe positions.

FIGURE 3 .

 3 FIGURE 3. (a) Probe positions for the illustration of TOFD simulations out of the focal zone; (b) L-wave signals simulated when probe are positioned symmetrically around the SDH; (c) Signals simulated when one probe is above the SDH.

FIGURE 4 .

 4 FIGURE 4. (a) Probe positions for the illustration of TOFD simulations out of the focal zone; (b) Signals simulated when the probe is above the SDH; (c) B-scans (time/position representations) for a 10mm horizontal scanning around that position.

TABLE 1 .

 1 Computation times in seconds.

		Example 1	Example 2	Example 3	Example 4
	Finite elements	786	1320	508	683
	Plane-wave approximation	0.8	0.8	0.9	2
	Ray-based method	1.2	4.7	4.4	8.6

position of the defect relative to the focal zone or to the part geometry. The most noticeable improvement is the shape of the signals. The maximum amplitudes are also generally more accurate with the ray-based model.

The improvement due to the new method comes at a cost in computation times, though that cost has been limited through optimization and parallelization. Ongoing work aims at further reducing computation times and at making the ray-based method available in a future release of Civa.