FEW REMARKS ON THE USE OF LOVE WAVES IN NON DESTRUCTIVE TESTING

. This paper concerns a theoretical study on the possibility of using Love waves for non destructive testing. A mathematical model is presented and analyzed. Several numerical tests are given in order to show the mechanical behaviour of this model.


1.
Introduction. The detection of cracks at the interface between two materials, for instance in pipes with a coating, is often treated using ultrasonic waves or Foucault currents propagating transversally to the interface. The detection is based on the analysis of the return signal which is different depending if the wave meets a crack or not. The possibility of using Love waves is a new idea suggested by many authors (for instance one can consult [8]- [10]- [11]- [13]- [15]- [18]). Such waves could enable one to explore larger area from a single shot. But the signal processing which should be implied is much more complicated. Therefore, a mathematical model can be useful in order to determine if there is a crack and where it is localized. One solution is to have a wave simulator model set on a finite part of the structure, but the boundary conditions should avoid artificial reflections. The goal of this paper is to make explicit these remarks and to suggest appropriate numerical transparent boundary conditions for this challenge. Nevertheless, this study is only a discussion on a mathematical model which could be used in order to improve the understanding of Love waves for NDT. Because the real problem in our example concerns cracks which are along the interface between two materials, this is the case that we consider in the following. But a large part of the analysis is still true for more general position of the crack. This will be discussed in a forthcoming paper. Concerning the details of the mathematical proofs and the application to an operational use, we refer the reader to the references. Let us consider a rectangular open set as the one drawn on figure 1. There are two subsets Ω − and Ω + corresponding to two different materials. The wave velocities are respectively c − and c + . For sake of clarity, it is assumed that c − < c + . The upper and lower boundaries -say Γ + and Γ − -are free edges and we set Γ 0 = Γ + ∪ Γ − . The lateral boundaries Γ e and Γ s are artificial cuts on which The global open set on which the model is set is Ω = Ω + ∪ Ω − ∪ {Γ i \ Γ f }. But the extension to several cracks could be discussed in the same framework without any new theoretical difficulty. A mechanical excitation is produced on a subdomain of Ω + or Ω − which is assumed to be remote from the crack tips. It can be an initial condition or a time dependent force denoted by f and such that: f (x, t) = g(t)η(x). The mathematical model that we discuss consists in finding u(x, t) solution for any t of the following system: The existence and uniqueness of a solution is not exactly the result contained in most books on partial differential equations, because of the dissipation term which appears on the boundaries Γ e and Γ s . The proof method consists in constructing a sequence of approximations of the model (1) using a Galerkin approximation with the variational formulation and the family of eigenfunctions w n , n ≥ 1 solution of: The first term w 1 is the constant function corresponding to λ 1 = 0. One can define a sequence of approximate solutions u n to the equation (1) using the space V N span by the N -first eigenvectors solution of (2). This is obtained from the variational formulation using test functions in V N . This solution can be written: and one can prove that it is bounded in the space L ∞ (]0, T [; H 1 (Ω)). By extracting a subsequence, one proves the existence and then the uniqueness of a solution to (1). The time regularity is obtained by taking the time derivative of the model. The functional space used for the space dependence of the solution is H 1 (Ω) and L 2 (Ω) for its time derivative. For instance, one has the following result: Theorem 1.1. Let us assume that u 0 ∈ H 1 (Ω), u 1 ∈ L 2 (Ω) and: Then there is a unique solution u to (1) such that: Remark 1. It is worth noting that the convergence of the series (3) can't be proved in the space H 2 (Ω). This can be justified by the fact that u n satisfies on the boundary of Ω the condition: which is incompatible with the boundary condition which should be satisfied by u at any time on Γ e ∪ Γ s .

2.
A more refined regularity analysis for u. Let us discuss locally what happens in the vicinity of points E (or/and F ) (see figure 2). On a close neighbourhood O of E (same is true for F ), the solution u can be locally split into its symmetrical and unsymmetrical parts with respect to the coordinate x 2 as shown on figure 2. They are respectively denoted by u s and u a . From a classical result, the function u s and its derivative with respect to x 2 are continous across Γ i . Therefore one can prove from the local equation satisfied by u s , that it is locally in the space H 2 (O).
The discontinuity of ∂u ∂x 1 is therefore taken into account by the unsymmetrical term ∂u a ∂x 1 as it is explained hereafter. Because: the function ∂u ∂x 1 is piecewise smooth as far as: which can't be discontinuous. Since the velocity c is discontinuous, ∂u ∂x 1 is also discontinuous. In fact, one can introduce the following function defined on O and which takes into account this discontinuity (c ∓ = c − , c ± = c + for x 2 > 0 and A primitive of q with respect to x 1 is a obtained using the so-called Dilog function [21]. More precisely [5], let us set: The singular fonction S involved in our problem is: From a simple computation one can check that this function is locally equivalent to Log(r) (r is the distance from E). Furthermore it can be truncated symmetrically around Γ i (the interface) such that its integral over Ω is zero. S and its derivative with respect to the coordinate x 2 have been plotted on figure 3.
Theorem 2.1. Let Q be a sub-open set strictly included in Ω + ∪ Ω − . In other words, Q doesn't cross the interface Γ i . Assuming that the initial data are such that u 0 ∈ H 2 (Ω) , with support in Q and u 1 ∈ H 1 (Ω), and that the right hand side satisfies f ∈ H 1 (]0, T [; L 2 (Ω)), one has: More generally, the following array using the same notations as above and the assumptions of theorem 2.1, describes the regularity of u in the different parts of Ω: With another respect, the time regularity can be derived up to any order assuming regularity on the initial data u 0 et u 1 and the right hand side f of the wave equation.
Remark 2. Let us consider a vertical line -say Γ v (corresponding to 0 < x 1 < L).
belongs to the space: Hence there is a finite limit for x 2 = 0 + and x 2 = 0 − . This property is no more true along the boundary Γ e ∪ Γ s because of the singularity S given at (6). This remark is meaningful in our discussion.
For instance, one has: and therefore: But c is discontinuous across Γ i and therefore one has: but (assuming an additional regularity on the initial data and on the right hand side of (1)), one has (and even more): 2.1. Effect of the transparent boundary condition on u. For sake of simplicity in this subsection, let us assume that f = 0. But this is just a facility, not a restriction. The mechanical energy of the structure is defined by: The mechanical energy is decreasing with respect to the time variable because there are only initial conditions and due to the transparent boundary condition which acts as a damping. Let us check it. One has, by multiplying the wave equation by ∂u ∂t and integrating over Ω×]0, t[: This proves two properties: There exists a constant c independent on time such that: Let us point out that the estimate (7) proves a somewhat hidden regularity concerning the term ∂u ∂t on the boundary Γ e ∪ Γ s . Let us introduce the sequence of functions u t0 (x, t) = u(x, t 0 + t). The previous estimate enables one to extract a subsequence which converges to a the term denoted by u * (with respect to the time t 0 and the uniqueness of the limit will prove that all the sequence converges) such that: From the convergence of E(t) to E ∞ one can state that: Therefore: This implies that: and finally, because of the lower semi-continuity of a continuous convex function: From the equations satisfied by u, one obtains that u * is solution of the limit model: The Holmgrem theorem [12] enables one to conclude that the function u * (in a first step one proves that the time derivative of u * is zero) is a constant with respect to space and time variables on ]0, T [×Ω. Such a constant doesn't generate wave. Hence when t → ∞ the waves disappears. One can get rid of this constant as follows. By integrating the wave equation (1) on ]0, t[×Ω, one obtains: If the initial conditions satisfy: one deduces that when t → ∞ and because u * = constant: A similar analysis can be derived if f = 0 but with some conditions on this term.
3. The Fourier transform for the operational use. In this paper, the identification of a crack is assumed to be performed a posteriori from measurements issued from the experiment. Hence the use of a Fourier transform of the wave equation is fully justified. Furthermore there is a mathematical equivalence between the time and the frequency domains. And the identification problem that we discuss is a mathematical one. For sake of simplicity in the notations, it is assumed in this section that the initial conditions on u are homogeneous and that u = 0 for t < 0. The only external force applied is f at the right hand side of the wave equation. Furthermore, it is assumed to decrease sufficiently fast when t → ∞ in order to make sense to the following. Let us set: e −iωt u(x, t)dt = u r + iu i , ( real and imaginary parts ofû). (12) One has: −ω 2û − div(c 2 ∇û) =f in Ω,

PHILIPPE DESTUYNDER AND CAROLINE FABRE
By splitting these equations into the real and imaginary parts, one deduces that: − 3.1. Property of the Fourier model. Let us considerû the Fourier transform of u defined at (12); (u is decreasing to zero when t → ∞ because of the transparent boundary condition which acts as a damping as shown in the previous section as soon as a compatibility condition is satisfied by the data): One can prove the following result.
Proof. Let us introduce the bilinear form a(..) defined on the space H 1 (Ω) × H 1 (Ω) and the linear form l(.) on H 1 (Ω) by: One has: Therefore the bilinear form a(., .) is Garding-coercive [16]. This enables one to apply the Fredholm alternative [17] ensuring that the model (14)- (15) which is equivalent to findũ ∈ H 1 (Ω) × H 1 (Ω) such that: has a unique solution as soon as ω = √ λ n , ∀n ∈ N. 2 From (17), it can be deduced that if one uses Shannon wavelet [14], the solution u is zero excepted for the frequencies in the bandwidth concerned.
On a neighbourhood V A of A, one can write: : and on a neighbourhood V B of B : The coefficients K A and K B are called the stress intensity factors and they depend on the time variable. Concerning the Fourier transform of u, the previous result becomes for instance around the extremity A: with:K The same result is also true for the extremity B of the crack tip.
5. An energy release rate for crack detection. Let us now introduce the socalled multiplier method for partial differential equations which in fact is nothing else but a domain derivative of the elastic energy with respect to a displacement of the domain Ω (see [7]). Several integrations by parts are used. Multiplying (15)-(14) by ∂u r ∂x 1 or ∂u i ∂x 1 and by integrating over Ω δ , one obtains the following relations (ν 1 is +1 on Γ s and −1 on Γ s : Let us discuss briefly this relation. This is an important point in the detection criterion of the crack located on the interface line Γ i . 1. There is a formula which gives the expression of the stress intensity factors with respect to the solutionû in a neighbourhood of the crack tip. Let us recall it ( [7]). First of all we denote by ∂V A the boundary of the neighbourhood of A as shown on figure 4. Let us denote byK Ar the stress intensity factor at point A for u r . If ν = (ν 1 , ν 2 ) is the unit normal to ∂V A inside V A , one has (a similar result is also true in the vicinity of point B): This formula enables one to compute accurately the stress intensity factor K Ar . 2. If there is no crack (K Ar =K Br = 0) one has always G = 0; 3. Let us notice that if: then (it is assumed that ω 1 = ω 2 ): , and if the function f is correctly chosen, one can prove that there is no crack. It looks like a controllability result. Let us give few ideas on the justification of this property. For sake of simplicity we only consider u r assuming for instance that u i = 0 and we restrict the discussion to the case c + = c − in order to avoid the two singularities appearing at E and F (see [6] for the complete case). In fact one can use the characterization of the stress intensity factor from the dual singular functions (one for each crack tip) [9] -say S * A and S * B -which leads to the expressions: In fact, one can use the eigenvector basis {w n } defined at (2) in order to represent the solution u r . Hence one has for any ω / ∈ Λ = {λ n } (see (2)): Thus, setting: one has:

NON DESTRUCTIVE TESTING AND LOVE WAVES 437
Let us consider that f r = χ [ω1,ω2] d(x). Let us introduce a function of the space L 2 (Ω) -say h(x, ω)-such that; If : : The expressions of S * A and S * B are known (for instance using the polar coordinates around the point A), one has [9]: Therefore, if d is chosen sucht that one can conclude that this is impossible for any position of A and B, one can claim that there is no crack. 5. Let us denote by G obs the quantity G but computed from an experimental measurement (or an estimation of it). This quantity is estimated from strain gauges set on the two boundaries Γ e and Γ s which give the terms required in the characterization of G obs . More precisely, the terms involved in G obs are: ∂u ∂x2 and u along these boundaries. Furthermore, because Love waves are suggested, and because they are exponentially decreasing inwards the domain Ω + , only a neighbourhood of [Γ e ∪ Γ s ] ∩ Ω − could be equipped. With another respect, for a given position of the crack (ie x A and x B given) one can compute the numerical expression of G(ω, x A , x B ) and compare it to G obs . Therefore, one can define an observability criterion of a crack by: And the localisation of the crack can be handle through the following optimisation problem for x C given between on ]0, L[ and (x A , x B ) ∈ R 2 being the control variables: The last term is introduced in order to select the solution for which the middle of the crack is the nearest of the point with abscissa x 1 = x C and ε is a small parameter (inducing a Tychonoff regularisation because of the a priori non uniqueness of the solution [20]). The continuity of CN D with respect to (x A , x B ) and its positiveness ensure the existence of a solution (but not the uniqueness). In fact, the main point is to detect the presence of a crack. Its exact position is a second step much more complicated from the computational point of view. It can performed numerically by using an adjoint state which leads to the expression of the gradient o the criterion to be minimized. This will be discussed in future works.
6. The numerical simulation in time. The Fourier transform can be computed directly or derived from a time simulation which is the choice we did. In fact we choose to solve numerically in time in order to have a full description of the mechanical phenomenon. Furthermore, in order to restrict the frequency range, we use a pseudo-Shannon wavelet setting:   7. The Fix method for the singularities at the points E and F . In order to improve a finite element scheme for solving a PDEs, Fix [19] has suggested in fracture mechanics modelling, to add the singular function S to the classical f.e.m. space V h spanned by first order polynomials (using here quadrangular elements). This leads to an approximate solution denoted by u h and such that (we restrict to the singularity on Γ e for instance, for sake of brevity): Our goal is to discuss this method in our framework. In particular, the error estimate between α(t) and α h (t) is analyzed in terms of error with respec to the mesh size say − h. Finally, this suggests to substract the term α h (t)S(x) to the approximate solution u h in order to try to eliminate the effect of this non-physical singularity in the computation of the criterion CN D(x A , x B ).  The variational formulation consists in finding u ∈ V 0 = {v ∈ H 1 (Ω) and such that: Let us denote by V h the approximation space spanned by the first degree Lagrange finite element with the classical assumptions. V h 0 is the subspace of functions v ∈ V h satisfying the additional condition Ω v = 0. Setting: The upgrade approximate model consists in finding u h ∈ V h 0s such that: Theorem 7.1. Let u ∈ V 0 and u h ∈ V h 0s be the solutions of (28) and (30). Let us set: The family of meshes is assumed to satisfy the classical regularity assumptions (see Ciarlet [4] or Raviart-Thomas [17]). Then there exists a constant c independent on the mesh size h such that: Indications on the proof following the method used in Amara-Destuynder-Djaoua [1]. Fom the variational formulation characterizing u at (28), one gets: Let us set (see figure 8): K h =]0, h[×] − h, h[, which are the two elements near the singular point E for instance. Furthermore, we introduce the projection -say P -from V 0 onto V h 0 with respect to the scalar product induced by the bilinear form a(., .). One has for instance (constant functions disappear in the expression of the bilinear form a(., .)): ∀v ∈ V h 0 , a(S − P S, v) = 0. Let us define by π the Lagrange interpolation operator from V ∩ C 0 (Ω) into V h . From the relations satisfied bu u and u h one has, using classical error estimates (see  and P.A. Raviart-J.M. Thomas [17] ): ≤ a(u − πu, u − πu) ≤ ch|u r | 2,2,Ω .
And by separating the term in V h 0 from the singular function S, one obtains: ∀v ∈ V h 0s , (α h − α)a(S, v) + a(u h r − u r , v) = 0, ∀v ∈ V h 0 , a(S − P S, v) = 0. Setting: v = S − P S,

PHILIPPE DESTUYNDER AND CAROLINE FABRE
one deduces that (with a new constant c): The computation of the previous infimum using a symbolic software (Maple) enables one to complete the proof , because S is known analytically. 2 Remark 3. The extension to the stationary case (ω = 0) is based on the same method as the one described quickly here. Once, α h known, the idea is to substract α h S from u h in order to avoid the effect of the singularity S in the computation of the detection criterion G defined at (20). The effect is to smooth the behaviour of the solution of the wave equation near the singular points.
Remark 4. This theoretical result suggests to use u h r in the CND criterion for detecting the crack. It appears has a way for eliminating the artifical singularity S even if the error estimate on α − α h is very coarse.
Another possibility is to project the finite element solution (without adding the singularity) orthogonally to the singularity S. This is done by setting for instance:

S.
We have plotted u h and R 0 u h on figure 10 the results obtained by this simple method in the computation of the term ∂u ∂x 2 . But up to now this is just an idea without mathematical foundations. Figure 10. Comparison between u h (left) and Ru h (right) inside Ω at the same time 8. Conclusion and future work. The non destructive testing is a challenge which requires efficient and reliable methods. The large amount of structures which should be tested requires new methods which enable to perform this duty faster than classical ones. The researches in this field are very active and a lot of problems are still to be solved from both the experimental aspect and the signal processing algorithms.
In this paper, we have only suggested few remarks based on mathematical modelling but which are to be continued in order to confirm the operational possibilities of the tool considered here.