domain: math spi

text:
Philippe Destuynder
email: philippe.destuynder@cnam.fr

Caroline Fabre
email: caroline.fabre@u-psud.fr

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

Keywords: 2010 Mathematics Subject Classification. Primary: 58F15, 58F17; Secondary: 53C35 Non destructive testing, Wave equation, Transparent boundary conditions, Fracture mechanics, Inverse problems

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 [START_REF] Diot | Mesures de défauts par ultrasons laser dans des soudures d'alliage d'aluminium[END_REF]- [START_REF] Dumont-Fillon | Contrôle non destructif par les ondes de Love et Lamb[END_REF]- [START_REF] Galvagni | The reflection of guided waves from simple supports in pipes[END_REF]- [START_REF] Lowe | Characteristics of the reflection of Lamb waves from defects in plates and pipes[END_REF]- [START_REF] Marty | Modelling of ultrasonic guided wave field generated by piezoelectric transducers[END_REF]- [START_REF] Ribichini | Study and comparison of different EMAT configurations for SH wave inspection[END_REF]). 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. 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 the lateral boundaries Γ e and Γ s are artificial cuts on which transparent boundary conditions are applied in order to simulate the propagation on an infinite strip. The separation line between the two media is denoted by Γ i . A crack may exist somewhere on this line. Its two extremities are positioned at points A and B with abscissa x 1 = x A and x 1 = x B . 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 FIGURE 1. The open set used for the wave model 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:

                   ∂ 2 u ∂t 2 -div (c 2 ∇u) = f in Ω, ∂u ∂ν = 0 on Γ -∪ Γ + , and ∂u ∂t + c ∂u ∂ν = 0 on Γ e ∪ Γ s , u(x, 0) = u 0 (x), ∂u ∂t (x, 0) = u 1 (x) in Ω. (1)
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:

                 -div (c 2 ∇w n ) = λ n w n in Ω, 0 ≤ λ 1 ≤ . . . λ n ≤ λ n+1 . . . , ∂w n ∂ν = 0 on Γ + ∪ Γ -∪ Γ e ∪ Γ s , Ω w 2 n (x)dx = 1. (2)
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:

i=1,N α n (t)w n (x), (3)
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 [START_REF] Amara | On a finite element schem for plane crack problems[END_REF]. 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: f ∈ L 2 (]0, T [×Ω).
Then there is a unique solution u to (1) such that:

u ∈ C 0 ([0, T ]; H 1 (Ω)) ∩ C 1 ([0, T ]; L 2 (Ω)).
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:

∂u n ∂ν = 0,
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.

∂u ∂t ∈ C 0 ([0, T ]; H 1/2 (Γ e ))
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 c ∓ = c + , c ± = c -for x 2 < 0): q(x 1 , x 2 ) = ( c ∓ c ± )( c + -c - π(c 2 + + c 2 -) )Log( 1 + 2 cos(x 2 π)e -πx1 + e -2x1π 1 -2 cos(x 2 π)e -πx1 + e -2x1π ) (4)
A primitive of q with respect to x 1 is a obtained using the so-called Dilog function [START_REF] Zagler | The Dilog function[END_REF].

More precisely [START_REF] Ph | Singularities occuring in multimaterials with transparent boundary conditions[END_REF], let us set:

D(z) = n≥1 z n n 2 with z = e -π(x1-ix2) . (5)
The singular fonction S involved in our problem is:

S(x 1 , x 2 ) = Im(D(z) -D(-z))/2, and q(x 1 , x 2 ) = ( c + -c - π(c 2 + + c 2 -) )( c ∓ c ± ) ∂S ∂x 2 . ( 6
)
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. and H 1 (Ω), and that the right hand side satisfies f ∈ H 1 (]0, T [; L 2 (Ω), with support in Q, one has:

u = α(t)S(x) + u R (x, t), with u R (x, t) ∈ W 1,∞ (]0, T [; H 1 (Ω)) ∩ W 2,∞ (]0, T [; L 2 (Ω)).
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 Ω:

Function → Open set ↓ | ∂u ∂t | ∂u ∂x 1 | ∂u ∂x 2 U ⊂⊂ Ω ± | L ∞ (]0, T [; H 1 (Ω)) | L ∞ (]0, T [; H 1 (U )) | L ∞ (]0, T [; H 1 (U )) Q ⊂⊂ Γ e± | L ∞ (]0, T [; H 1 2 (Γ e )) | L ∞ (]0, T [; H 1 2 (Q)) | L ∞ (]0, T [; H 1 2 (Q)) P ⊂⊂ Γ i | L ∞ (]0, T [; H 1 2 (Γ i )) | L ∞ (]0, T [; H 1 2 (P )) | L ∞ (]0, T [; H -1 2 (Γ i ))
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). Along

Γ + v (x 2 > 0) or Γ - v (x 2 < 0) the function ∂u ∂x 2
belongs to the space:

L ∞ (]0, T [; H 3/2 (Γ + v ∪ Γ - v )) (it is discontinuous across Γ i ).
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 [START_REF] Destuynder | Numerical analysis of singularities appearing in bimaterials with transparent boundary conditions[END_REF]. This remark is meaningful in our discussion.

For instance, one has:

∂u ∂t ∈ L ∞ (]0, T [; H 1/2 (Γ e ∪ Γ s ))
and therefore:

c ∂u ∂ν |Γe∪Γs ∈ L ∞ (]0, T [; H 1/2 (Γ e ∪ Γ s )).
But c is discontinuous across Γ i and therefore one has:

∂u ∂ν |Γe∪Γs / ∈ L ∞ (]0, T [; H 1/2 (Γ e ∪ Γ s )),
but (assuming an additional regularity on the initial data and on the right hand side of (1)), one has (and even more):

∂u ∂ν |Γe∪Γs ∈ L ∞ (]0, T [; L 2 (Γ e ∪ Γ s )).
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:

E(t) = 1 2 Ω | ∂u ∂t | 2 (x, t)dx + 1 2 Ω c 2 |∇u| 2 (x, t)dx.
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[:

E(t) + t 0 Γe∪Γs c| ∂u ∂t | 2 (x 2 , t)dx 2 dt ≤ E(0) ∀t > 0. (7)
This proves two properties:

1. the function t → E(t) which is continuous, is decreasing when t → ∞. Hence there is a lower limit -say E ∞ -such that lim t→∞

E(t) = E ∞ ;
2. There exists a constant c independent on time such that:

∀t > 0, || ∂u ∂t || 0,2,Ω ≤ c, ||∇u|| 1,Ω ≤ c.
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:

                 u t0 (x, t) → t0→∞ u * in L ∞ (]0, T [; H 1 (Ω)) weakly, ∂u t0 ∂t (x, t) → t0→∞ ∂u * ∂t in L ∞ (]0, T [; L 2 (Ω)) weakly, ∂u t0 ∂t (x, t) → t0→∞ ∂u * ∂t in L 2 (]0, T [×{Γ e ∪ Γ s })) weakly. (8)
From the convergence of E(t) to E ∞ one can state that:

∀ε > 0, ∃t ε > 0, ∀t ≥ t ε : E(t) ≤ E ∞ + ε. (9)
Therefore:

∀t ≥ t ε , E ∞ ≤ lim inf t→∞ E(t) + ∞ tε Γe∪Γs c| ∂u ∂t (x 2 , t)| 2 dx 2 dt ≤ E(t 0 ) ≤ E ∞ + ε.
This implies that:

∀ε > 0, ∃t ε ∀t 1 ≥ t ε , t1 tε Γe∪Γs c| ∂u ∂t (x 2 , t)| 2 dx 2 dt ≤ ε,
and finally, because of the lower semi-continuity of a continuous convex function:

∂u * ∂t = 0 on ]0, T [×{Γ e ∪ Γ s }. ( 10
)
From the equations satisfied by u, one obtains that u * is solution of the limit model :

         ∂ 2 u * ∂t 2 -div(c 2 ∇u * ) = 0 in ]0, T [×Ω, ∂u * ∂ν = 0, on ]0, T [×Γ 0 , ∂u * ∂t = 0 and ∂u * ∂ν = 0 on ]0, T [×Γ e ∪ Γ s . (11)
The Holmgrem theorem [START_REF] Holmgrem | Über Systeme von linearen partiellen Differentialgleichungen, Öfversigt af Kongl[END_REF] 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:

∀t > 0, ∂ ∂t ( Ω u(x, t)dx 1 dx 2 ) + Γe∪Γs cu(x, t)dx 2 ) = Ω u 1 (x)dx 1 dx 2 + Γe∪Γs cu 0 (x)dx 2 .
If the initial conditions satisfy:

Ω u 1 + Γe∪Γs cu 0 = 0,
one deduces that when t → ∞ and because u * = constant:

Γe∪Γs cu * = 0 => u * = 0.
A similar analysis can be derived if f = 0 but with some conditions on this term.

3. The Fourier transform for the operational use. 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:

û(ω, x) = ∞ 0
e -iωt u(x, t)dt = u r + iu i , ( real and imaginary parts of û).

One has:

-ω 2 û -div(c 2 ∇û) = f in Ω, ∂ û ∂ν = 0, on Γ 0 , iωû + c ∂ û ∂ν = 0 on Γ e ∪ Γ s . (13)
By splitting these equations into the real and imaginary parts, one deduces that:

-ω 2 u r -div(c 2 ∇u r ) = f r in Ω, ∂u r ∂ν = 0, on Γ 0 , ωu i = c ∂u r ∂ν on Γ e ∪ Γ s . (14)
-ω 2 u i -div(c

2 ∇u i ) = f i in Ω, ∂u i ∂ν = 0, on Γ 0 , ωu r + c ∂u i ∂ν = 0 on Γ e ∪ Γ s . (15)
3.1. Property of the Fourier model. Let us consider û the Fourier transform of u defined at [START_REF] Galvagni | The reflection of guided waves from simple supports in pipes[END_REF]; (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:

                         ∀ũ = (u r , u i ), ṽ = (v r , v i ) ∈ H 1 (Ω) × H 1 (Ω), a(ũ, ṽ) = -ω 2 Ω u r v r + u i v i + Ω c 2 (∇u r .∇v r + ∇u i .∇v i ) + ω Γe∪Γs c(u r v i -u i v r ), l(ṽ) = Ω (f r v r + f i v i )(x, ω)dx. ( 16
)
One has:

a(ṽ, ṽ) = -ω 2 Ω (v 2 r + v 2 i ) + Ω c 2 (|∇v r | 2 + |∇v i | 2 ).
Therefore the bilinear form a(., .) is Garding-coercive [START_REF] Necas | Les méthodes directes en théorie des équations elliptiques[END_REF]. This enables one to apply the Fredholm alternative [START_REF] Raviart | Approximation des équations aux dérivées partielles[END_REF] ensuring that the model ( 14)- [START_REF] Marty | Modelling of ultrasonic guided wave field generated by piezoelectric transducers[END_REF] which is equivalent to find ũ ∈ H 1 (Ω) × H 1 (Ω) such that: ∀ṽ, a(ũ, ṽ) = l(ṽ),

has a unique solution as soon as ω = √ λ n , ∀n ∈ N. ✷ From ( 17), it can be deduced that if one uses Shannon wavelet [21], the solution û is zero excepted for the frequencies in the bandwidth concerned.

4.

Local behaviour of the solution (u r , u i ) near the crack tip. Let us assume that there is a crack on the interface Γ i as shown on figure 1. The two extremities are denoted by A and B and their abscissa on Γ i are respectively x 1 = a and x 1 = b. From classical results (see [?], [START_REF] Bui | Mécanique de la rupture fragile[END_REF]) one can write in the vicinity of each crack tip ((R A , θ A ), -respectively (r B , θ b )are the polar coordinates from A -respectively Bthe reference axis for the polar angle being the line Γ i as shown on figure 4):

                                        
On a neighbourhood V A of A, one can write: :

u(x, t) = K A (t)S A (x) + u RA , u RA ∈ L ∞ (]0, T [; H 2 (V A )), with S A (x) = √ r A c sin( θ A 2 ), c = c + if θ A > 0, c = c -if θ A < 0,
and on a neighbourhood V B of B :

u(x, t) = K B (t)S B (x) + u RB , u RB ∈ L ∞ (]0, T [; H 2 (V B )), with S B (x) = √ r B c sin( θ A 2 ), c = c + if θ B > 0, c = c -if θ B < 0. (18)
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:

û(x, ω) = KA (ω)S A (x) + ûRA , (19)
with:

KA (ω) = ∞ 0 e -iωt K A (t)dt and ûRA (x, ω) = ∞ 0 e -iωt u(x, t)dt.
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 so-called 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 [START_REF] Ph | Sur une interpretation mathématique de l'intégrale de Rice en mécanique de la rupture fragile[END_REF]).

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 :

         - Γe∪Γs [ω 2 (u 2 r + u 2 i )ν 1 + c 2 2 (| ∂u r ∂x 2 | 2 + | ∂u i ∂x 2 | 2 )]ν 1 - Ω (f r ∂u r ∂x 1 + f i ∂u i ∂x 1 ) = G(ω) = π 4 (| KB | 2 -| KA | 2 ) (20)

PHILIPPE DESTUYNDER AND CAROLINE FABRE

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 ( [START_REF] Ph | Sur une interpretation mathématique de l'intégrale de Rice en mécanique de la rupture fragile[END_REF]).

First of all we denote by ∂V A the boundary of the neighbourhood of A as shown on figure 4. Let us denote by KAr 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):

π 4 | KAr | 2 = ∂V A [- ω 2 2 u 2 r ν ′ 1 -c 2 ∂u r ∂ν ′ ∂u r ∂x 1 + c 2 2 |∇u r | 2 ν ′ 1 ]. ( 21
)
This formula enables one to compute accurately the stress intensity factor K Ar .

2. If there is no crack ( KAr = KBr = 0) one has always G = 0;

3. Let us notice that if:

G(ω) = 0, ∀ω ∈ [ω 1 , ω 2 ]
then (it is assumed that ω1 = ω 2 ):

| KA (ω)| = | KB (ω)|, ∀ω ∈ [ω 1 , ω 2 ]. 4. If G(ω) = 0 on [ω 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 [START_REF] Destuynder | On the detection of cracks in linear elasticity research report CNAM[END_REF] 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) [START_REF] Dobrowolski | Numerical Approximation of Elliptic Interface and Corner Problems[END_REF] -say S * A and S * B -which leads to the expressions:

K Ar = Ω [f r S * A + ω 2 u r S * B ]dx - Ω div(c 2 ∇S Ar )S * A , K Br = Ω [f r S * B + ω 2 u r S * B ]dx - Ω div(c 2 ∇S Br )S * B . (22)
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)):

u r (x, ω) = n≥1 ξ n (ω)w n (x) with ξ n (ω) = Ω f r w n λ n -ω 2 .
Thus, setting:

δ A = 1 - Ω div(c 2 ∇S A )S * A , δ B = 1 - Ω div(c 2 ∇S B )S * B ,
one has:

K Ar (ω)=δ A n≥1 λ n Ω f r w n Ω S * A w n λ n -ω 2 K Br (ω)=δ B n≥1 λ n Ω f r w n Ω S * B w n λ n -ω 2 .
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;

h(x, ω) = n≥1 [ λ n λ n -ω 2 Ω d(x)w n (x)dx]w n (x).
For instance

h(x, 0) = d(x).
If :

∀ω ∈ [ω 1 , ω 2 ], |K Ar | = |K Br |, one has, assuming that f (x, ω) = χ [ω1,ω2] d(x), one can deduce: Ω h(ω)S * A = ± Ω h(ω)S * B = 0.
The expressions of S * A and S * B are known (for instance using the polar coordinate around the point A one has [START_REF] Dobrowolski | Numerical Approximation of Elliptic Interface and Corner Problems[END_REF]:

S * A (r A , θ A ) = 1 c √ r A sin( θ A 2 ).
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). 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:

CN D(x A , x B ) = ω2 ω1 |G(ω, x A , x B ) -G obs |) 2 dω. (23)
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:

min 0≤x A ≤x C ≤x B ≤L CN D(x A , x B ) + ε[(x A -x C ) 2 + (x B -x C ) 2 ]. (24)
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 Tikhonov regularisation because of the a priori non uniqueness of the solution [START_REF] Tikhonov | [END_REF]). 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 usinf 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 [START_REF] Strang | Analysis of the Finite Elements Method[END_REF] 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:

f (x, t) = d(x)[ sin(at) πt - sin(bt) πt ], a > b, (25)
u h (x, t) = α h (t)S(x) + u h r (x, t), u h r (x, t) ∈ V h . ( 26
)
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 ).

For sake of simplicity, we discuss the case of the following model (the functions f and g are given):

         -div(c 2 ∇u) = f in Ω assuming that: Ω f + Γe∪Γi gv = 0, ∂u ∂ν = 0 on Γ 1 , ∂u ∂ν = g on Γ e ∪ Γ s , Ω u = 0 (27)
The variational formulation consists in finding u ∈ V 0 = {v ∈ H 1 (Ω) and

Ω v = 0}
such that:

∀v ∈ V 0 , Ω c 2 ∇u.∇v = Ω f v + Γe∪Γs gv. (28)
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:

V h s = {R × S} ⊕ V h , and V h 0s = {v ∈ V h s , Ω v = 0}. (29)
The upgrade approximate model consists in finding u h ∈ V h 0s such that:

∀v ∈ V h 0s , Ω c 2 ∇u h .∇v = Ω f v + Γe∪Γs gv. (30)
Theorem 7.1. Let u ∈ V 0 and u h ∈ V h 0s be the solutions of ( 28) and (30). Let us set:

u = αS + u r , u h = α h S + u rh where u r ∈ H 2 (Ω) ∩ V 0 , u rh ∈ V h 0 .
The family of meshes is assumed to satisfy the classical regularity assumptions (see Ciarlet [START_REF] Ciarlet | The finite element mehod for elliptic problems[END_REF] or Raviart-Thomas [START_REF] Raviart | Approximation des équations aux dérivées partielles[END_REF]). Then there exists a constant c independent on the mesh size h such that:

             ||u -u h || 1,2,Ω ≤ ch|u r | 2,2,Ω and |α -α h | ≤ c Log(h) |u r | 2,2,Ω .
Indications on the proof following the method used in Amara-Destuynder-Djaoua [START_REF] Amara | On a finite element schem for plane crack problems[END_REF]. 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 Strang-Fix [START_REF] Strang | Analysis of the Finite Elements Method[END_REF] and P.A. Raviart-J.M. Thomas [START_REF] Raviart | Approximation des équations aux dérivées partielles[END_REF] ):

∀v ∈ V h 0s , a(u h -u, v) = 0 => a(u h -u, u h -u) ≤ 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, one deduces that (with a new constant c):

                 |α h -α| = ch
a(S -P S, S -P S)

≤ ch 1 inf z=ax1+bx2+cx1x2 K h c 2 | ∂(S -z) ∂x 2 | 2
.

The computation of the previous infimum using a symbolic software (Maple) enables one to complete the proof , because S is known analytically. ✷ 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 [START_REF] Tikhonov | [END_REF]. 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:

R 0 u h = u h -Ω Su h Ω S 2 S or R 1 u h = u h -Ω c 2 ∇u h .∇S Ω c 2 |∇S| 2 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. 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.

FIGURE 3 .

 3 FIGURE 3. The singular function S on left and ∂S/∂x 2 on right

Theorem 3 . 1 .

 31 Let us denote by (λ n , w n ) the eigenvalues and eigenfunctions solution of the model[START_REF] Brezis | Analyse fonctionnelle[END_REF]. Let us set Λ = { √ λ n , n ∈ N * }. For any ω ∈ R \ Λ, the coupled system (15) and (14) has a unique solution.

FIGURE 4 .

 4 FIGURE 4. Local polar coordinates

FIGURE 5 .

 5 FIGURE 5. Right hand side (left); solution with c + = 2c -(right)In order to compare the influence of the soft media we have drawn on figure6two cases of ratio between c + and c -. One can see the trapping effect in the softest media (Love waves) which increases with the ratio c + /c -compared to the figure 5 on the right.

FIGURE 6 .

 6 FIGURE 6. Left : c = c - Right: c + = 4c -

FIGURE 7 .FIGURE 8 .

 78 FIGURE 7. Without crack Solution of (1): low frequencies (left), ∂u/∂x 2 (right); f is given at (25): a = 50 b = 10. Comparing with figure 8 one can see that the energy is less localized in the softest media than with high frequencies.

FIGURE 9 .

 9 FIGURE 9. With a crack on the interface between the two materials The left figure represents the wave equation solution. The right one is ∂u/∂x 2 on Γ e and on Γ s . Same conditions as on figure8but with a crack between x 1 = .4L and x 1 = .6L. One can see the influence of the crack on this quantity and therefore on the criterion CND defined at (23)

FIGURE 10 .time 8 .

 108 FIGURE 10. Comparison betwwen u h (left) and Ru h (right) inside Ω at the same time

Acknowledgments. This work is developed in the context of the French-English GDR 2501: Wave propagation in complex media for quantitative and non destructive evaluation. The support of this research group from both scientific exchanges and experimental point of view has been particularly useful in our study.


Samuel Shannon

Victor Peron

Zohar Yosibash

Singular asymptotic solution along an elliptical edge for the Laplace equation in 3-D

Keywords: 3-D singularities, elliptical singular edge, edge flux intensity functions

An explicit asymptotic solution to the elasticity system in a three-dimensional domain in the vicinity of an elliptical crack front, or for an elliptical sharp V-notch is still unavailable. Towards its derivation we first consider the explicit asymptotic solutions of the Laplace equation in the vicinity of an elliptical singular edge in a three-dimensional domain. Both homogeneous Dirichlet and Neumann boundary conditions on the surfaces intersecting at the elliptical edge are considered. The dual singular solution is also provided to be used in a future study to extract the edges flux intensity functions by the quasi-dual-function method. We show that just as for the circular edge case, the solution in the vicinity of an elliptical edge is composed of three series, with eigenfunctions being functions of two coordinates.

Introduction

Asymptotic solutions of the elasticity system or the Laplace equation in the vicinity of singular points over two dimensional domains have been investigated for over half a century. For the Email address: zohary@bgu.ac.il (Zohar YOSIBASH) Preprint submitted to ??? November 7, 2014 Laplace equation these are described by an infinite series of the eigenpairs (α k , φ k (ϕ)) and their coefficients named flux intensity factors A k (FIFs) as follows:

τ (ρ, ϕ) = k A k ρ α k φ k (ϕ). (1)
In three-dimensional domains such as polyhedra, both vertex and edge singularities exist, see [START_REF] Dauge | Elliptic boundary value problems in corner domains -smoothness and asymptotics of solutions[END_REF][START_REF] Petersdorff | Decompositions in edge and corner singularities for the solution of the Dirichlet problem of the Laplacian in a polyhedron[END_REF]. For straight edges an explicit representation of the singular solutions along the edge was provided in [START_REF] Hartranft | The use of eigenfunction expansions in the general solution of three-dimensional crack problems[END_REF][START_REF] Costabel | A quasidual function method for extracting edge stress intensity functions[END_REF][START_REF] Omer | Edge flux intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF][START_REF] Yosibash | Edge stress intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF]. It is given by two infinite series containing the 2-D series and another series of the derivative of the FIFs A k which are functions along the straight singular edge as follows:

τ (ρ, ϕ, z) = k ∞ ℓ=0,2,4,6••• ∂ ℓ z A k (z)ρ α k +ℓ φ k,ℓ (ϕ). (2)
It was shown in [START_REF] Yosibash | Circular edge singularities for the Laplace equation and the elasticity system in 3-D domains[END_REF] that for three-dimensional domains with a circular singular edges, the solution of the Laplace equation is given by three series:

τ (ρ, ϕ, θ) = k ∞ ℓ=0,2,4,6••• ∂ ℓ θ A k (θ)ρ α k ρ R ℓ ∞ i=0 ρ R i φ k,ℓ,i (ϕ), (3)
where θ is the angle surround the circular edge, ρ is the normal distance from the circular edge and R is the radius of the circular edge. A k (θ) are the Edge Flux Intensity Functions (EFIFs), Sadowsky and Sternberg [START_REF] Sadowsky | Stress concentration around a triaxial ellipsoidal cavity[END_REF] considered the stress around a 3-D ellipsoidal cavity, but these ellipsoids do not contain singular edges. Green and Sneddon were the first (to the best of our knowledge) to analyze an elliptical crack in an infinite body under uniform tension [START_REF] Green | The distribution of stress in the neighbourhood of a flat elliptical crack in an elastic solid[END_REF]. They used the solutions in [START_REF] Sadowsky | Stress concentration around a triaxial ellipsoidal cavity[END_REF] to compute the normal displacement to the crack plane. Irwin [START_REF] Irwin | Crack-extension force for a part-through crack in a plate[END_REF] derived the pointwise edge stress intensity factor K I based on [START_REF] Green | The distribution of stress in the neighbourhood of a flat elliptical crack in an elastic solid[END_REF] for an elliptical crack in an infinite body under uniform tension. Kassir and Sih computed mode II and mode III edges stress intensity functions K II and K III for an elliptical crack in an infinite body under uniform shear (see [START_REF] Kassir | Three-dimensional stress distribution around an elliptical crack under arbitrary loadings[END_REF]). These solutions have been used in many papers for computing the edge stress intensity functions by the superposition method ( [START_REF] Shah | Streess intensity factor for an elliptical crack under arbitrary normal loading[END_REF][START_REF] Shah | Streess intensity factor for an elliptical crack approaching the surface of a semi-infinite solid[END_REF]), the weight function method (for example see [START_REF] Nishioka | The first-order variation of the displacement field due to geometrical changes in an elliptical crack[END_REF][START_REF] Singh | Closed form solution for an annular elliptic crack around an elliptic rigid inclusion in an infinite solid[END_REF]) and by Fourier approximation of the boundary condition or of the geometry (see [START_REF] Vijayakumar | An embedded elliptical crack, in an infinite solid, subject to arbitrary crack-face tractions[END_REF][START_REF] Roy | An elliptic crack in an elastic half-space[END_REF][START_REF] Singh | Closed form solution for an annular elliptic crack around an elliptic rigid inclusion in an infinite solid[END_REF][START_REF] Livieri | An analysis of three-dimensional planar embedded cracks subjected to uniform tensile stress[END_REF]). Leblond and Torlai provided the mechanism for the point-wise derivation of the elastic solution up to second order for a general curved crack ( [START_REF] Leblond | The stress-field near the front of an arbitrarily shaped crack in a 3-dimensional elastic body[END_REF]). Much attention has been devoted to the computation of the first coefficient of the asymptotic series in elasticity, namely the pointwise evaluation of the mode I stress intensity function along an elliptical crack. In [START_REF] Atroshchenko | Stress intensity factors for elliptical and semi-elliptical cracks subjected to an arbitrary mode I loading[END_REF] the exact analytical pointwise stress intensity function has been obtained for an elliptical crack embedded in an infinite elastic body and subjected to an arbitrary applied normal stress (Mode I). However, the general asymptotic solution of the Laplace equation and the elasticity system in the vicinity of an elliptical singular edge is yet unknown. As a result there are no methods that may provide the flux/stress intensity factors along the edge as functions of the coordinate along it. Towards this goal we herein consider first the general expression of the asymptotic solution to the Laplace equation in a domain that contains an elliptical singular edge.

α
For the Laplace equation we show that as for the circular edge case, the solution in the vicinity of an elliptical edge is composed of three series, but here the eigenfunctions φ k,ℓ,i are functions of the coordinate ϕ and of the coordinate γ along the elliptical edge. Specific examples for an elliptical crack with homogeneous Dirichlet or Neumann BCs are provided. The asymptotic solution and especially the dual solution is required to extract the edge flux intensity function (EFIFs) using the quasi-dual function method as in [START_REF] Costabel | A quasidual function method for extracting edge stress intensity functions[END_REF][START_REF] Yosibash | Circular edge singularities for the Laplace equation and the elasticity system in 3-D domains[END_REF]. Therefore we explicit provide in Appendix Appendix A the first terms of the dual solution. This is the first step for the computation of the solution to the elasticity system in the vicinity of elliptical singular edges which is of major engineering importance because most surface cracks are semi-elliptical.

The Laplace equation in the vicinity of an elliptic singular edge

Consider a singular edge of an elliptic shape in the x 1 -x 2 plane. We use a standard orthogonal curvilinear elliptic coordinate system (β, γ) related to the Cartesian coordinates by:

x 1 = a cosh(β) cos(γ) (4)
x 2 = a sinh(β) sin(γ) (5)
where β ≥ 0 , γ ∈ [0, 2π] , and ±a are the ellipse foci. A planar ellipse is determined by a fixed β = β 0 and a given a . The outer normal vector n at any point along such an ellipse (determined by γ ) is:

n = √ 2 cosh (2β 0 ) -cos(2γ) (sinh(β 0 ) cos(γ)x 1 , cosh (β 0 ) sin(γ)x 2 , 0) (6)
Remark 1. The ellipse reduces to a circle of a radius R by substituting a = R cosh(β 0 ) or a = R sinh(β 0 ) at the limit β 0 → ∞ in (4)- [START_REF] Hartranft | The use of eigenfunction expansions in the general solution of three-dimensional crack problems[END_REF].

We locate the two orthogonal coordinates (ρ, ϕ) initiating at the elliptical edge in the plane containing the normal n as shown in Fig. 1. The connection between the coordinates ρ, γ, ϕ and the Cartesian coordinates x 1 , x 2 , x 3 is:

x 1 = a cosh(β 0 ) + ρ cos(ϕ) √ 2 sinh(β 0 ) cosh(2β 0 ) -cos(2γ) cos(γ) (7)
x 2 = a sinh(β 0 ) + ρ cos(ϕ) √ 2 cosh(β 0 ) cosh(2β 0 ) -cos(2γ) sin(γ) (8)
x 3 = ρ sin(ϕ) (9)
and the scale factors (also known as the Lamé factors) are:

h ρ = 1 (10)
h ϕ = ρ (11)
h γ = a q(γ) g(γ) + ρ a f (ϕ) (12)
with:

f (ϕ) = sinh (2β 0 ) cos(ϕ) (13)
q(γ) = cosh (2β 0 ) -cos(2γ) (14)
g(γ) = q 3 (γ)/2 (15)
The Laplace equation in the new coordinate system (ρ, ϕ, γ) located on the ellipse edge (see Figure 1 ) is:

∆ {ρ,ϕ,γ} = ∂ ρρ + 1 ρ ∂ ρ + 1 ρ 2 ∂ ϕ,ϕ (16)
+ 1 a sinh (2β 0 ) g(γ) + ρ a f (ϕ) cos(ϕ)∂ ρ - 1 ρ sin(ϕ)∂ ϕ - 1 a 2 1 g(γ) + ρ a f (ϕ) 2 q(γ) sin(2γ)∂ γ + ρ a 3 1 g(γ) + ρ a f (ϕ) 3 q(γ) sin(2γ)3f (ϕ)∂ γ + 1 a 2 q(γ) g(γ) + ρ a f (ϕ) 2 ∂ γγ Remark 2.
For the limit case of a circular crack a = R cosh(β 0 ) or a = R sinh(β 0 ) and β 0 → ∞ :

lim β 0 → ∞ 1 a 1 g(γ) + ρ a f (ϕ) sinh(2β 0 ) = 1 R 1 1 + ρ R cos(ϕ) lim β 0 → ∞ 1 a 2 1 g(γ) + ρ a f (ϕ) 2 q(γ) sin(2γ) = 0 lim β 0 → ∞ ρ a 3 1 g(γ) + ρ a f (ϕ) 3 3q(γ) sin(2γ)f (ϕ) = 0 lim β 0 → ∞ 1 a 2 q(γ) g(γ) + ρ a f (ϕ) 2 = 1 R 2 1 1 + ρ R cos(ϕ) 2
and (16) reduces to the Laplace equation along a circular edge as presented in [START_REF] Yosibash | Circular edge singularities for the Laplace equation and the elasticity system in 3-D domains[END_REF]).

Multiplying ∆τ = 0 by 1 + ρ a g(γ) f (ϕ)

3 ρ 2 = 0 , (16) becomes:

1 + ρf (ϕ) ag(γ) 3 ρ 2 × ∆ = M 0 + ρ sinh (2β 0 ) ag(γ) (3 cos ϕM 0 + ρM 01 ) (17)
+ ρ sinh (2β 0 ) ag(γ) 2 3 cos 2 ϕM 0 + 2 cos ϕM 01 + 1 sinh (2β 0 ) 2 -q(γ) sin(2γ)∂γ + q(γ) 2 ∂γγ + ρ sinh (2β 0 ) ag(γ) 3 cos 3 ϕM 0 + cos 2 ϕM 01 + cos ϕ sinh (2β 0 ) 2 2q(γ) sin(2γ)∂γ + q(γ) 2 ∂γγ
with:

M 0 = (ρ∂ ρ ) 2 + ∂ ϕϕ (18)
M 1 = cos(ϕ)ρ∂ ρ -sin(ϕ)∂ ϕ (19)
Following the Laplace solution for the circular edge (see [START_REF] Yosibash | Circular edge singularities for the Laplace equation and the elasticity system in 3-D domains[END_REF]), we consider a solution of (17) of the form:

τ = k ∞ ℓ=0 ∂ ℓ γ A k (γ)ρ α k ρ q(γ) a g(γ) ℓ ∞ i=0 ρ sinh (2β 0 ) a g(γ) i φ k,ℓ,i (ϕ, γ) (20)
with either homogeneous Dirichlet or Neumann BCs on Γ 1 and Γ 2 (see Fig. 2): Remark 3. For the limit case of circular crack:

lim β 0 → ∞ ρ sinh (2β 0 ) a g(γ) i = ρ R i lim β 0 → ∞ ρ q(γ) a g(γ) ℓ = ρ R ℓ and (20) reduces to (3).
The homogeneous Dirichlet boundary conditions on ϕ 1 and ϕ 2 in view of (20) are:

φ k,ℓ,i (ϕ 1 , γ) = φ k,ℓ,i (ϕ 2 , γ) = 0 on Γ 1 ∪ Γ 2 ∀ γ ∈ [0, 2π], (21)
∂ ϕ φ k,ℓ,i (ϕ 1 , γ) = ∂ ϕ φ k,ℓ,i (ϕ 2 , γ) = 0 on Γ 1 ∪ Γ 2 ∀ γ ∈ [0, 2π]. (22)
In Figure 3 Substituting [START_REF] Yosibash | Circular edge singularities for the Laplace equation and the elasticity system in 3-D domains[END_REF] in ( 17) one obtains:

0 = A k (γ)ρ α k m 0 φ k,0,0 (23)
+ ρ sinh (2β 0 ) a g(γ) m 0 φ k,0,1 + (3 cos(ϕ) m 0 + m 01 ) φ k,0,0 + ρ sinh (2β 0 ) a g(γ) 2 m 0 φ k,0,2 + (3 cos(ϕ) m 0 + m 01 ) φ k,0,1 + 3 cos 2 (ϕ) m 0 + 2 cos(ϕ) m 01 φ k,0,0 + q(γ) sinh 2 (2β 0 ) (-sin(2γ)∂γ + q(γ)∂γγ ) φ k,0,0 + ρ sinh (2β 0 ) a g(γ) 3 m 0 φ k,0,3 + (3 cos(ϕ) m 0 + m 01 ) φ k,0,2 + 3 cos 2 (ϕ)m 0 + 2 cos(ϕ)m 01 φ k,0,1 + cos 3 (ϕ)m 0 + cos 2 (ϕ) m 01 φ k,0,0 + q(γ) sinh (2β 0 ) 2 g(γ) (-sin(2γ)∂γ + q(γ)∂γγ ) φ k,0,1 g(γ) + cos(ϕ) (2 sin(2γ)∂γ + q(γ)∂γγ ) φ k,0,0 + • • • + A ′ k (γ)ρ α k ρ q(γ) a g(γ) m 0 φ k,1,0 + ρ sinh (2β 0 ) a g(γ) m 0 φ k,1,1 + (3 cos(ϕ)m 0 + m 01 ) φ k,1,0 + 1 sinh (2β 0 ) (-sin(2γ) + 2q(γ)∂γ ) φ k,0,0 + ρ sinh (2β 0 ) a g(γ) 2 m 0 φ k,1,2 + (3 cos(ϕ)m 0 + m 01 ) φ k,1,1 + 3 cos 2 (ϕ) m 0 + 2 cos(ϕ)m 01 φ k,1,0 + 1 sinh (2β 0 ) -sin(2γ)φ k,0,1 + 2q(γ)g(γ)∂γ φ k,0,1 g(γ) + 2 cos(ϕ) (sin(2γ) + q(γ)∂γ ) φ k,0,0 + g(γ) sinh (2β 0 ) 2 (-sin(2γ)∂γ + q(γ)∂γγ ) q(γ)φ k,1,0 g(γ) + • • • (24) + A ′′ k (γ)ρ α k ρ q(γ) a g(γ) 2 m 0 φ k,2,0 + φ k,0,0 + ρ sinh (2β 0 ) a g(γ) m 0 φ k,2,1 + (3 cos(ϕ) m 0 + m 01 ) φ k,2,0 + φ k,0,1 + cos(ϕ)φ k,0,0 + 1 sinh (2β 0 ) -sin(2γ)φ k,1,0 + 2g(γ)∂γ q(γ)φ k,1,0 g(γ) + • • • + A ′′′ k (γ)ρ α k ρ q(γ) a g(γ) 3 m 0 φ k,3,0 + φ k,1,0 + • • • + • • •
where m 0 and m 1 are defined by:

m 0 φ k,l,i ∆ = (α + l + i) 2 φ k,l,i + ∂ ϕϕ φ k,l,i (25)
m 01 φ k,l,i ∆ = (α + l + i) cos(ϕ)φ k,l,i -sin(ϕ)∂ ϕ φ k,l,i (26)
Equation ( 23) has to hold true for any (ρ/a) i and for any ∂ ℓ γ A k (γ) , resulting in the following recursive set of ODEs for the determination of φ k,ℓ,i (ϕ, γ) :

m 0 φ k,0,0 = 0 (27)
m 0 φ k,0,1 = -(3f (ϕ)m 0 + m 1 ) φ k,0,0 (28)
m 0 φ k,0,2 = -(3 cos(ϕ)m 0 + m 01 ) φ k,0,1 -cos(ϕ) (3 cos(ϕ)m 0 + 2m 01 ) φ k,0,0 (29)
+ q(γ) sinh 2 (2β 0 ) (sin(2γ)∂γ -q(γ)∂γγ ) φ k,0,0 m 0 φ k,0,3 = -(3 cos(ϕ)m 0 + m 01 ) φ k,0,2 -cos(ϕ) (3 cos(ϕ)m 0 + 2m 01 ) φ k,0,1 -cos 2 (ϕ) (cos(ϕ)m 0 + m 01 ) φ k,0,0 (30)
q(γ) sinh 2 (2β 0 ) g(γ) (-sin(2γ)∂γ + q(γ)∂γγ ) φ k,0,1 g(γ) + cos(ϕ) (2 sin(2γ)∂γ + q(γ)∂γγ ) φ k,0,0

m 0 φ k,1,0 = 0 (31) m 0 φ k,1,1 = -(3 cos(ϕ)m 0 + m 01 ) φ k,1,0 - 1 sinh (2β 0 ) (-sin(2γ) + 2q(γ)∂γ ) φ k,0,0 (32)
m 0 φ k,1,2 = -(3 cos(ϕ)m 0 + m 01 ) φ k,1,1 -cos(ϕ) (3 cos(ϕ)m 0 + 2m 01 ) φ k,1,0 (33)
- 1 sinh (2β 0 ) -sin(2γ)φ k,0,1 + 2q(γ)g(γ)∂γ φ k,0,1 g(γ) + 2 cos(ϕ) (sin(2γ) + q(γ)∂γ ) φ k,0,0 - g(γ) sinh 2 (2β 0 ) (-sin(2γ)∂γ + q(γ)∂γγ ) q(γ)φ k,1,0 g(γ) m 0 φ k,2,0 = -φ k,0,0 (34)
m 0 φ k,2,1 = -(3 cos(ϕ)m 0 + m 01 ) φ k,2,0 -φ k,0,1 + cos(ϕ)φ k,0,0 (35)
- 1 sinh (2β 0 ) -sin(2γ)φ k,1,0 + 2g(γ)∂γ q(γ)φ k,1,0 g(γ) m 0 φ k,3,0 = -φ k,1,0 (36)
α k = 0, 1/2, 1, 3/2, 2, • • • .
Notice that equations ( 27)-(30) depend only on ϕ , so:

φ k,0,0 (ϕ, γ) = φ k,0,0 (ϕ) (37)
φ k,0,1 (ϕ, γ) = φ k,0,1 (ϕ) (38)
φ k,0,2 (ϕ, γ) = φ k,0,2 (ϕ) (39)
For simplicity, we can choose the homogeneous solution of ( 31) and (36) as zero, so:

φ k,1,0 (ϕ, γ) = 0 (40) φ k,3,0 (ϕ, γ) = 0 (41)
and ( 27)-(36) reduce to:

m 0 φ k,0,0 = 0 (42) m 0 φ k,0,1 = -(3 cos(ϕ)m 0 + m 01 ) φ k,0,0 (43)
cos 2 (ϕ) (cos(ϕ)m 0 + m 01 ) φ k,0,0 + q(γ)

sinh 2 (2β 0 ) g(γ) (sin(2γ)∂ γ -q(γ)∂ γγ ) φ k,0,1 g(γ) m 0 φ k,1,1 = sin(2γ) sinh (2β 0 ) φ k,0,0 (46)
m 0 φ k,1,2 = -(3 cos(ϕ)m 0 + m 01 ) φ k,1,1 (47)
+ 1 sinh (2β 0 ) sin(2γ) (φ k,0,1 -2 cos(ϕ)φ k,0,0 ) -2q(γ)g(γ)∂ γ φ k,0,1 g(γ) m 0 φ k,2,0 = -φ k,0,0 (48)
m 0 φ k,2,1 = -(3 cos(ϕ) m 0 + m 01 ) φ k,2,0 -φ k,0,1 + cos(ϕ)φ k,0,0 (49)

A specific example: An elliptical Crack with homogeneous Dirichlet BCs

For elliptical cracks with homogeneous Dirichlet BCs the eigenvalues are α k = 1/2, 1, 3/2, 2, • • • . Equations ( 42)-(48) can be solved with the homogeneous Dirichlet BCs (21) for the different α k to find φ k,ℓ,i (ϕ, γ) . For example, the singular part of the solution (i.e. for α 1 = 1/2 ) up to order of (ρ) 3 is: 42)-( 48) can be solved with homogeneous Neumann BCs (22) for the different α k to find φ k,ℓ,i (ϕ, γ) . For example, the singular part of the solution (i.e. for α 1 = 1/2 ) up to order of (ρ) 3 is:

τ = A 1 (γ)ρ 1/2 cos ϕ 2 - ρ sinh (2β 0 ) a g(γ) 1
) sinh (2β 0 ) 2 cos ϕ 2 + • • • +A ′ 1 (γ)ρ 1/2 ρ q(γ) a g(γ) sin(2γ) sinh (2β 0 ) ρ sinh (2β 0 ) a g(γ) 1 6 cos ϕ 2 - ρ sinh (2β 0 ) a g(γ) 2 3 8 cos ϕ 2 + 4 15 cos 3ϕ 2 + • • • A ′′ 1 (γ)ρ 1/2 ρ q(γ) a g(γ) 2 - 1 6 cos ϕ 2 + ρ sinh (2β 0 ) a g(γ) 1
τ = A 1 (γ)ρ 1/2 sin ϕ 2 + ρ sinh (2β 0 ) a g(γ) 1

Summary and Conclusions

We provided the machinery to develop an explicit asymptotic solution for the Laplace equation in the vicinity of an elliptical crack or sharp V-notch in a three dimensional domain. Specifically, for an elliptical crack with homogeneous Dirichlet or Neumann boundary conditions we provided the asymptotic series up to an order of ρ 3 . We have also shown that at the limit when the ellipse turns into a circle the asymptotic solution simplifies to the one obtained for a penny shaped crack in [START_REF] Yosibash | Circular edge singularities for the Laplace equation and the elasticity system in 3-D domains[END_REF]. The asymptotic dual solution (associated with the negative eigenvalue -1/2 ) has been also derived and will serve for the extraction of the flux intensity functions using the quasi dual function method.

The machinery and solution for the Laplace equation serves as a cornerstone for the asymptotic solution for the elasticity system, and attempts are underway to obtain an asymptotic solution to the elasticity system in the vicinity of an elliptical and semi-elliptic crack.

  k are the 2-D eigenvalues and φ k,ℓ,i (ϕ) are the 2-D eigenfunctions and their shadow functions. This expansion contains the 3-D straight edge series and another series associated to the curvature of the singular edge. Explicit asymptotic solutions in the vicinity of singular elliptical edges in a 3-D domain are not yet available neither for the elasticity system nor for the Laplace equation. The elasticity asymptotic solution has same characteristics as the Laplace asymptotic solution and towards its derivation we here consider first the Laplace equation in the vicinity of an elliptical singular edge in a three-dimensional domain. Thereafter same techniques may be used to derive the asymptotic solution for the elasticity system.

Figure 1 :

 1 Figure 1: An ellipse in the x-y plane and the notations.

Figure 2 :

 2 Figure 2: A slice of the ellipse singular edge.

Figure 3 :

 3 Figure 3: The surfaces that intersect at the elliptic singular edge on which homogeneous boundary conditions are prescribed. Left -The entire ellipse, Right -Zoom in the right part ( γ = 0 ).

  with the BCs (21) or (22) is the well known eigenvalues problem in a 2-D Vnotched domain. For example, for a crack with homogeneous Neumann BCs the eigenvalues are:

m

  0 φ k,0,2 = -(3 cos(ϕ)m 0 + m 01 ) φ k,0,1 -cos(ϕ) (3 cos(ϕ) m 0 + 2m 01 ) φ k,0,0 (44) m 0 φ k,0,3 = -(3 cos(ϕ)m 0 + m 01 ) φ k,0,2 -cos(ϕ) (3 cos(ϕ) m 0 + 2m 01 ) φ k,0,1

  An elliptical Crack with homogeneous Neumann BCs For elliptical cracks with homogeneous Neumann BCs the eigenvalues are α k = 0, 1/2, 1, 3/2, 2, • • • . Equations (

	8	cos	ϕ 2	+	7 60	cos	3ϕ 2	+ • • •
	+O ρ 4							(50)
	2.2. A specific example:

Acknowledgements

S. Shannon and Z. Yosibash gratefully acknowledge the support by the Israel Science Foundation (grant No. 593/14).

Appendix A. The dual solution for an elliptical crack

Here we provide the first dual solution (associated with the negative eigenvalues α k < 0 ) for an elliptical crack with homogenoues Dirichlet or Neumann BCs on the crack faces.

Appendix A.1. An elliptical crack with homogenous Dirichlet BCs

For an elliptical crack with homogenous Dirichlet BCs the first dual eigenvalue is: α = -1/2 . Substituting α = -1/2 in equations ( 42)-(48) with homogeneous Dirichlet BCs (21) one obtains the first terms of the dual solution K (α 1 =-1 2 ) :

. An elliptical crack with homogeneous Neumann BCs

For an elliptical crack with homogenous Neumann BCs the first dual eigenvalue is: α = -1/2 . Substituting α = -1/2 in equations ( 42)-( 48) with the homogeneous BCs (22) one obtains the first terms of the dual solution K (α 1 =-1 2 ) :