Marc Bonnet
email: bonnet@lms.polytechnique.fr

Bojan B Guzina

Sounding of finite solid bodies by way of topological derivative

Keywords: Topological derivative, cavity identification, elastodynamics, adjoint fiend method, inverse scattering, shape sensitivity

This paper is concerned with an application of the concept of topological derivative to elastic-wave imaging of finite solid bodies containing cavities. Building on the approach originally proposed in the (elastostatic) theory of shape optimization, the topological derivative, which quantifies the sensitivity of a featured cost functional due to the creation of an infinitesimal hole in the cavity-free (reference) body, is used as a void indicator through an assembly of sampling points where it attains negative values. The computation of topological derivative is shown to involve an elastodynamic solution to a set of supplementary boundary-value problems for the reference body, which are here formulated as boundary integral equations. For a comprehensive treatment of the subject, formulas for topological sensitivity are obtained using three alternative methodologies, namely i) direct differentiation approach, ii) adjoint field method, and iii) limiting form of the shape sensitivity analysis. The competing techniques are further shown to lead to distinct computational procedures. Methodologies (i) and (ii) are implemented within a BEM-based platform and validated against an analytical solution. A set of numerical results is included to illustrate the utility of topological derivative for 3D elastic-wave sounding of solid bodies; an approach that may perform best when used as a pre-conditioning tool for more accurate, gradient-based imaging algorithms. Despite the fact that the formulation and results presented in this investigation are established on the basis of a boundary integral solution, the proposed methodology is readily applicable to other computational platforms such as the finite element and finite difference techniques.

INTRODUCTION

In recent times, rapid advances in sensor and signal processing technologies have exposed the need for efficient 3D elastodynamic solutions to inverse scattering problems. Similar to their electromagnetic counterpart, elastic waves find applications as a sensing tool in diverse areas of engineering such as medical imaging, seismic surveys, and nondestructive material testing. In general, inverse elastodynamic solutions are derived from either of the three computational cornerstones that include i) (far-field) ray theory [1], ii) finite-difference approximation of the wave equation [START_REF] Sheriff | Exploration Seismology[END_REF], and iii) boundary integral equation (BIE) formulations [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF]. In the context of three-dimensional imaging, these approaches bear substantial computational cost associated with solving the forward elastic-scattering problem. This precludes the use of global search techniques such as genetic algorithms which entail a large number of forward simulations. To mitigate the problem, gradient-based optimization algorithms have been proposed as a computationally-tractable alternative to solving inverse scattering problems, especially when ameliorated by the analytical shape sensitivity estimates [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF][START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF]. Unfortunately, the latter class of solutions necessitate a reliable prior information about the location, topology and geometry of a hidden scatterer (e.g. defect) for satisfactory performance. As a result, their stand-alone use may not be sufficient for 3D elastodynamic imaging which intrinsically involves minimization in a highdimensional parametric space.

Driven by the foregoing considerations, the focus of this study is a robust, yet computationally efficient approach for preliminary elastic-wave imaging of 3D solid bodies (containing cavities) based on the concept of topological derivative. With reference to a generic cost functional J of the body shape, its topological derivative, T (x o ), synthesizes the sensitivity of J with respect to the creation of an infinitesimal cavity at a prescribed location x o inside the reference, i.e. cavity-free counterpart of the probed body. The concept of topological derivative first appeared in [START_REF] Eschenauer | Bubble method for topology and shape optimization of structures[END_REF] and [START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Kopchpositionierungkrieterien[END_REF] in the context of topological optimization of mechanical structures (see [8] for a description of the subject), wherein the spatial distribution of T (x o ) was used in an iterative procedure as a criterion for the removal of "excess" material through regions where T < 0. Recently, its rigorous mathematical formulation has been established within the framework of elastostatic problems and Laplace equation [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF][START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF]. Beyond its direct application to the topology and shape optimization of structures, however, the topological derivative is also expected to facilitate the minimization-based solution of inverse scattering problems by providing, through a first-order approximation, a rational basis for selecting the number, location, and geometry of hidden defects, all of which are essential for establishing a reliable initial "guess" (see also [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF][START_REF] Jackowska-Strumillo | Topological optimization and inverse problems[END_REF] for the developments in 2D elastostatics). In the previous study [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], this idea was considered in the context of inverse elastic scattering pertaining to semi-infinite and infinite domains, where the availability of suitable fundamental solutions made it possible to establish explicit expressions for T (x o ).

The present study is a part of ongoing research by the authors on the identification of cavities in a finite 3D body from non-invasive elastodynamic measurements. The topological derivative is hence established in connection with a generic cost function expressed in the form of an integral over the body's external boundary and whose kernel involves the solution of the forward elastic scattering problem. This format constitutes a natural generalization of the least-squares misfit function, assumed in the ensuing numerical study. Owing to the fact that the fundamental solutions for a finite body of arbitrary shape are generally unavailable, it is no longer possible to obtain explicit expressions for T (x o ). Instead, the computation of T (x o ) is shown to entail the computation of the solutions of a set of supplementary boundary-value problems for the reference (cavity-free) body. In an effort to provide a comprehensive treatment of the subject, formulas for topological derivative are obtained using: i) the direct differentiation approach, ii) adjoint field approach, and iii) shape sensitivity analysis. The point of departure for the first two methods is an asymptotic expansion of the featured cost functional with respect to the creation of an infinitesimal hole in the otherwise intact reference solid. This in turn requires the knowledge of the leading contribution of the scattered field for a vanishing cavity. By means of a boundary integral analysis, the latter solution is found to be governed by a pair of separable problems: a) an elastostatic exterior problem for the scaled cavity in an elastic free space, and b) an elastodynamic interior problem for the reference body involving point sources acting at x o . With such result, the first approach yields a formula for topological derivative involving the solutions to problems a) and b), in a fashion similar to the direct differentiation approach for shape sensitivity analyses. The second approach leads to a distinct formula, involving the solutions to problem a) and an adjoint field defined on the cavity-free reference body. The latter formula for topological derivative is also obtained via the third approach, this time by taking a zero-size limit of the suitable shape-sensitivity expression for a finite cavity (e.g. [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF]). One should mention that the foregoing three methodologies have been employed elsewhere, e.g. in [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] for inverse scattering in semi-infinite media, and in [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF] and [START_REF] Éa | The shape and topological optimization connection[END_REF] for static optimization problems, as a tool for obtaining T . To the authors' knowledge, however, their parallel development and comparison for three-dimensional elastodynamics and a common computational platform have eluded the previous studies. A set of numerical results for a basic 3D configuration is included to examine the utility of topological derivative for preliminary elastic-wave imaging of solid bodies. Notwithstanding the fact that the formulation and results presented in this study are established on the basis of a BIE solution, the proposed methodology is readily applicable to other computational platforms such as the finite element method. The present study emphasizes the development of topological sensitivity as a preparatory tool for more accurate, gradient-based solution algorithms. In such framework, the computational cost entailed by the topological sensitivity calculation is of the order of one forward elastodynamic solution, and therefore modest compared to that of the subsequent minimization of the output least squares functional. Another possibility, not pursued here, consists of defining an imaging algorithm entirely based on the topological derivative, in which the preliminary computation of T would be followed by an iterative procedure, each step involving matter removal based on the current distribution of T and re-computation of T for the new domain.

PRELIMINARIES

Let Ω true denote a finite elastic body bounded by the external surface S and characterized by the shear modulus µ, Poisson's ratio ν and mass density ρ. An unknown cavity (or a set thereof) B true bounded by the closed surface(s) Γ true is embedded in Ω true , so that Ω true = Ω \ B true where Ω denotes the reference, i.e. cavity-free counterpart of Ω true . On applying a steady-state traction p on S with angular frequency ω, an elastodynamic state u true arises which satisfies the field equations and boundary conditions

∇•(C : ∇u true ) = -ρω 2 u true (ξ ∈ Ω true ), t true = p (ξ ∈ S N ), t true = 0 (ξ ∈ Γ true ), u true = 0 (ξ ∈ S D ).
(1) where ξ denotes the position vector. Here, the external surface S is divided into complementary subsets S N and S D supporting prescribed tractions and displacements, respectively; t true ≡ σ true • n = (C : ∇u true )•n denotes the traction vector associated with the displacement u true through Hooke's law, and

C = 2µ ν 1 -2ν I 2 ⊗ I 2 + I sym 4
is the isotropic elasticity tensor with I 2 and I sym 4

symbolizing the second-order and symmetric fourthorder identity tensors, respectively. For simplicity, it is assumed that ω is not an eigenfrequency of any of the boundary-value problems appearing in the ensuing developments.

For the inverse scattering problem of interest, where the location, topology and geometry of B true (or equivalently Γ true ) is being sought, the trace of u true on S is assumed to be available over the measurement region S obs ⊂ S N . In what follows, these measurements will be denoted by u obs , so that u obs (ξ, ω) = u true (ξ, ω), ξ ∈ S obs .

To facilitate the ensuing developments, let u c denote the solution to the forward problem for a given excitation p and a trial cavity B c bounded by Γ. Accordingly, the forward solution u c is defined over Ω c = Ω \ Bc , and is governed by the equations

∇•(C : ∇u c ) = -ρω 2 u c (ξ ∈ Ω c ), t c = p (ξ ∈ S N ), t c = 0 (ξ ∈ Γ). u c = 0 (ξ ∈ S D ).
(

) 2
where t c is the traction vector associated with u c .

To establish a rational framework for solving the inverse problem, a misfit function is set up in order to minimize the difference between u obs and u c . The weighted least-squares misfit function, commonly used for such purpose, is defined here as

J (Ω c ; p) = 1 2 S obs (u c -u obs (ξ)) • W (ξ) • (u c -u obs (ξ)) dΓ ξ (3)
where W (ξ) is a 3×3 matrix-valued weighting function, assumed to be symmetric and positive definite (the simplest choice being W (ξ) = I 2 ), while overbar denotes complex conjugation. In general, (3) can be considered as a special case of the class of cost functions permitting the format

J (Ω c ; p) = S obs ϕ(u c (ξ), ξ) dΓ ξ (4)
where

ϕ(w, ξ) = 1 2 (w(ξ)-u obs (ξ)) • W (ξ) • (w(ξ)-u obs (ξ)) (5)

Governing integral equation for the forward problem

Let ûk i (ξ, x, ω) and tk i (ξ, x, ω; n) = σk ij (ξ, x, ω)n j (ξ) denote respectively the i-th Cartesian components of the elastodynamic fundamental solution, i.e. of the displacement and traction vectors at ξ due to a unit (time-harmonic) point force acting at x in the k-th direction in an infinite homogeneous elastic solid characterized by the shear modulus µ, Poisson's ratio ν and mass density ρ. This fundamental solution can be decomposed into a singular part and a regular (residual) component via

ûk i (ξ, x, ω) = ûk i (ξ, x) + [û k i (ξ, x, ω)] 2 , tk i (ξ, x, ω; n) = tk i (ξ, x; n) + [ tk i (ξ, x, ω; n)] 2 , (6)
where ûk i (ξ, x) and tk i (ξ, x; n), which constitute the singular part, are given by the elastostatic fundamental solution for an infinite solid (i.e. the Kelvin solution). With such definitions, the forward problem (2) can be formulated in terms of a boundary integral equation (BIE):

S∪Γ u c i (ξ, ω) -u c i (x, ω) tk i (ξ, x; n) dΓ ξ + S∪Γ u c i (ξ, ω)[ tk i (ξ, x, ω; n)] 2 dΓ ξ = S t c i (ξ)û k i (ξ, x, ω) dΓ ξ (x ∈ Γ ∪ S) (7)
with u c = 0 on S D and t c = p on S N .

TOPOLOGICAL DERIVATIVE

To aid the gradient-based minimization of (4) that is often used as a tool for identifying B true on the basis of u obs , of interest here is the development of topological derivative for the cost functionals J (Ω; p) of form (4), which would facilitate rational selection of the necessary initial "guess" in terms of the location, topology and geometry of B true . To this end, let B ε (x o ) = x o + εB, where B ⊂ R 3 is a fixed bounded open set with boundary S and volume |B| containing the origin, define the region of space occupied by a cavity of (small) size ε > 0 containing a fixed sampling point x o . Following [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF][START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF], one is in particular interested in the asymptotic behavior of J (Ω ε ; p) for infinitesimal ε > 0,

where

Ω ε = Ω \ B ε (x o ), and B ε (x o ) is the closure of B ε (x o ).
With reference to this limiting behavior, the topological derivative T (x o , ω) of the cost functional J (Ω; p) at x o for a cavity-free body is defined through the expansion:

J (Ω ε ; p) = J (Ω; p) + ε 3 |B| T (x o , ω) + o(ε 3 ) (ε ≪ Diameter(Ω), B ε (x o ) ⊂ Ω) (8)
One may note that this definition is not restricted to spherical infinitesimal cavities (for which B is the unit ball, S the unit sphere and |B| = 4π/3). In general, the value T (x o , ω) is expected to depend on the shape of B. Also, the postulated O(ε 3 ) asymptotic behavior of J (Ω ε ; p) -J (Ω; p) is chosen for simplicity because it is known to hold for traction-free cavities in the 3-D elastostatic case [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF][START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF]. The analysis to follow will corroborate this choice, and would allow to find this behavior were it left unspecified in definition (8).

With reference to (8), the evaluation of J (Ω ε ; p) requires the knowledge of the elastodynamic solution u ε to the forward problem (2) with B c replaced by B ε ≡ B ε (x o ). To this end, it is convenient to decompose the total displacement field u ε as

u ε = u + ũε (9)
where ũε denotes the scattered field, and u is the free field defined as the response of the void-free (reference) solid Ω due to given excitation (i.e. boundary traction) p, so that

∇•(C : ∇u) = -ρω 2 u (ξ ∈ Ω), t = p (ξ ∈ S N ), u = 0 (ξ ∈ S D ), ( 10
)
where t is the traction vector associated with u. At this point, it is important to emphasize that the notions of void-free reference configuration and the assocated free field, postulated for clarity reasons, do not restrict the applicability of topological sensitivity to more general geometric configurations. For example, the techniques developed in this study are equally applicable to reference solids containing pre-existing cavities and, in particular, to iterative imaging algorithms where the voids identified in a previous iteration are used to update the reference configuration for the next step. By analogy to [START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Kopchpositionierungkrieterien[END_REF], the free field is governed by the integral equation

u k (x, ω) + S u i (ξ, ω) -u i (x, ω) tk i (ξ, x; n) dΓ ξ + S u i (ξ, ω)[ tk i (ξ, x, ω; n)] 2 dΓ ξ = S t i (ξ)û k i (ξ, x, ω) dΓ ξ (x ∈ S) (11)
with u = 0 on S D and t = p on S N . For infinitesimal ε the scattered field is expected to vanish, i.e.

lim ε→0 |ũ ε (x)| = 0 (x ∈ Ω c ) (12)
whereas the free-field, by its definition [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF], does not depend on ε. On the basis of ( 4) and ( 12), one may expand J (Ω ε ; p) with respect to ũε as

J (Ω ε ; p) = S obs ϕ(u ε (ξ), ξ) dΓ ξ = S obs ϕ(u(ξ), ξ) + Re ∂ϕ ∂u u(ξ), ξ • ũε (ξ) + o(|ũ ε (ξ)|) dΓ ξ = J (Ω; p) + S obs Re ∂ϕ ∂u u(ξ), ξ • ũε (ξ) dΓ ξ + o( ũε ) (13)
where

∂ϕ ∂w ≡ ∂ϕ ∂w R -i ∂ϕ ∂w I w R = Re(w) , w I = Im(w) (14)
By means of ( 8) and ( 13), the topological derivative of J (Ω; p) can be recast as:

T (x o , ω) = lim ε→0 1 ε 3 |B| S obs Re ∂ϕ ∂u u(ξ), ξ • ũε (ξ) dΓ ξ . (15)
To estimate the leading perturbation term

S obs Re ∂ϕ ∂u u(ξ), ξ • ũε (ξ) dΓ ξ
one can resort to either i) a direct approach, ii) an adjoint field method, or iii) limiting form of the shape sensitivity analysis. These alternative methodologies for computing the topological sensitivity are presented in Sections 4, 5 and 6, respectively.

DIRECT APPROACH

The direct approach for the computation of topological derivative proposed in this study represents a generalization of the methodology developed in [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] for the inverse scattering problems in infinite and semi-infinite media. It entails a substitution of the leading contribution of ũε into (15) and seeking the limit of the resulting expression as ε → 0. It is therefore essential for this purpose to find the leading asymptotic behavior of ũε as ε → 0.

Leading contribution of ũε as ε → 0

From (2) with u c replaced with u ε , ( 9) and [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF], it follows that the scattered field ũε solves the boundary-value problem:

∇•(C : ∇ũ ε ) = -ρω 2 ũε (ξ ∈ Ω ε ), tε = 0 (ξ ∈ S N ), ũε = 0 (ξ ∈ S D ), tε = -σ•n (ξ ∈ Γ ε ) (16)
where σ = C : ∇u is the stress tensor associated with the free field [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF], Γ ε is the boundary of B ǫ , and n is the normal on S ∪ Γ ε outward to Ω ε . As a result, ũε is governed by a boundary integral equation analogous to [START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Kopchpositionierungkrieterien[END_REF] but with i) Γ superseded by Γ ε , and ii) applied loading being given by the surface traction -σ•n distributed over Γ ε instead of p acting on S N . To investigate the limiting form of this integral equation as ε → 0, it is convenient to distinguish between the two distinct cases, namely x ∈ S and x ∈ Γ ε (x being the collocation point), so that the boundary traces of (ũ ε , tε ) on S and Γ ε can be defined in terms of a pair of integral equations:

ũε k (x, ω) + Γε ũε i (ξ, ω) -ũε i (x, ω) tk i (ξ, x; n) dΓ ξ + Γε ũε i (ξ, ω)[ tk i (ξ, x, ω; n)] 2 dΓ ξ + S ũε i (ξ, ω) tk i (ξ, x, ω; n) dΓ ξ - S D tε i (ξ, ω)û k i (ξ, x, ω) dΓ ξ = - Γε t i (ξ, ω)û k i (ξ, x, ω) dΓ ξ (x ∈ Γ ε ) (17)
and

S ũε i (ξ, ω) -ũε i (x, ω) tk i (ξ, x; n) dΓ ξ + S ũε i (ξ, ω)[ tk i (ξ, x, ω; n)] 2 dΓ ξ + Γε ũε i (ξ, ω) tk i (ξ, x, ω; n) dΓ ξ = S D tε i (ξ, ω)û k i (ξ, x, ω) dΓ ξ - Γε t i (ξ, ω)û k i (ξ, x, ω) dΓ ξ (x ∈ S) (18)
where t i = σ ij n j is the traction associated with the free field. One may note the absence of the free term (ũ ε k (x, ω)) from BIE [START_REF] Sokolowski | Introduction to shape optimization[END_REF] owing to the fact that S is an external boundary of Ω ε (e.g. [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF]). To find the equations governing the leading contribution of (ũ ε , tε ) as ε → 0, one has to establish the asymptotic form of the system ( 17) and [START_REF] Sokolowski | Introduction to shape optimization[END_REF]. For this purpose, scaled coordinates

x = (x -x o )/ε and ξ = (ξ -x o )/ε (where x, ξ ∈ S ) are introduced when x ∈ Γ ε or ξ ∈ Γ ε , respectively. In particular, one has dΓ ξ = ε 2 dϑξ (ξ ∈ Γ ε , ξ ∈ S ) (19)
where dϑ denotes a differential element of the surface S bounding B ⊂ R 3 .

Preliminary estimate of (ũ ε , tε ) on S. In what follows, the traces of ũε on S N and Γ ε and of tε on S D will be denoted respectively as ũε S N , ũε Γε and tε S D . In view of [START_REF] Jackowska-Strumillo | Topological optimization and inverse problems[END_REF], it is natural to assume tentatively that ũε Γε = O(ε d ), where the exponent d > 0 is to be found. Based on this assumption, (18) can be used to estimate the behavior of ũε S N and tε S D as ε → 0. For such purpose, one may note that

Γε ũε i (ξ, ω) tk i (ξ, x, ω; n) dΓ ξ = ε 2 σk ij (x o , x, ω) S ũε i (x o + ε ξ, ω)n j ( ξ) dϑξ + o(ε d+2 ) (20)
and

Γε t i (ξ, ω)û k i (ξ, x, ω) dΓ ξ = ûk i (x o , x, ω)σ ij (x o , ω) Γε n j (ξ) dΓ ξ =0 + ûk i (x o , x, ω)σ ij (x o , ω) ,ℓ Γε n j (ξ)(ξ ℓ -x 0 ℓ ) dΓ ξ =-δ jℓ ε 3 |B| + o(ε 3 ) = ρω 2 ûk i (x o , x, ω)u i (x o , ω) -ûk i,j (x o , x, ω)σ ij (x o , ω) ε 3 |B| + o(ε 3 ) (x ∈ S) (21)
On the other hand, the remaining terms in [START_REF] Sokolowski | Introduction to shape optimization[END_REF] have the form

L S ũε S N , tε S D
where the linear integral operator L S clearly does not depend on ε. On the basis of this result, [START_REF] Dominguez | Boundary elements in dynamics[END_REF], and ( 21), integral equation ( 18) can be recast as

L S ũε S N , tε S D = O ε min(d+2, 3) as ε → 0 (22)
It follows that ũε S N and tε S D (i.e. the traces of ũε on S N and tε on S D ) are O(ε min(d+2, 3) ), i.e. that one can represent the leading contributions of ũε S N and tε S D as ε → 0 in terms of auxiliary functions V (ξ) and T (ξ) which are independent of ε, so that:

ũε S N (ξ) = ε min(d + 2, 3) |B| V (ξ) + o(1) , |V | = O(1) (ξ ∈ S N ) tε S D (ξ) = ε min(d + 2, 3) |B| T (ξ) + o(1) , |T | = O(1) (ξ ∈ S D ) (23)
Asymptotic form of equation [START_REF] Wheeler | Some theorems in classical elastodynamics[END_REF]. One may note that (17) involves integrals over S and Γ ε . By virtue of (23a), the integrals over S in [START_REF] Wheeler | Some theorems in classical elastodynamics[END_REF] are readily seen to be of order 3) .

S ũε i (ξ, ω) tk i (ξ, x, ω; n) dΓ ξ = O ε min(d+2, 3) S tε i (ξ, ω)û k i (ξ, x, ω) dΓ ξ = O ε min(d+2,
(

) 24
For the integrals over Γ ε in [START_REF] Wheeler | Some theorems in classical elastodynamics[END_REF], on the other hand, one has x ∈ Γ ε and ξ ∈ Γ ε . In this case, one can show that

ûk i (ξ, x, ω) = 1 ε ûk i ( ξ, x) + O(1) tk i (ξ, x, ω; n) = 1 ε 2 tk i ( ξ, x; n) + O(1) , x, ξ ∈ Γ ε (25)
and

[û k i (x o + ε ξ, x o + εx, ω)] 2 = O(1) [ tk i (x o + ε ξ, x o + εx, ω; n)] 2 = O(1) , x, ξ ∈ Γ ε (26)
where ûk i ( ξ, x) and tk i ( ξ, x; n) are respectively the displacement and traction given by the Kelvin fundamental solution. Now, again assuming that ũε

Γε = O(ε d ) with d > 0, one has Γε ũε i (ξ, ω)[ tk i (ξ, x, ω; n)] 2 dΓ ξ = O(ε d+2 ) (27)
by virtue of ( 19) and (26), and

Γε t i (ξ, ω)û k i (ξ, x, ω) dΓ ξ = εσ ij (x o , ω) S ûk i ( ξ, x)n j ( ξ) dϑξ + O(ε 2 ) (28)
by virtue of ( 19) and (25). On introducing the scaled coordinates x and ξ into the remaining integrals featured in [START_REF] Wheeler | Some theorems in classical elastodynamics[END_REF] and taking into account estimates (24), ( 27) and (28), equation ( 17) is readily seen to take the form

ũε k (x o + εx) + S ũε i (x o + ε ξ) -ũε i (x o + εx) tk i ( ξ, x; n) dϑξ = -εσ ij (x o , ω) S ûk i ( ξ, x)n j ( ξ) dϑξ + O(ε 2 ) , x ∈ S (29)
On the basis of (29), integral equation ( 17) is therefore of the form

L S ũε Γε = O(ε) as ε → 0
where the linear integral operator L S is O(1). It immediately follows that d = 1, i.e. that one can seek the solution ũε k (x o + εx) of (29) in the form ũε

Γε (ξ) = εU ( ξ) + o(ε) (ξ ∈ Γ ε ) (30)
where the auxiliary field U ( ξ) does not depend on ε and solves the counterpart of (29)

U k (x) + S U i ( ξ) -U i (x) tk i ( ξ, x; n) dϑξ = -σ ij (x o , ω) S ûk i ( ξ, x; n)n j ( ξ) dϑξ (x ∈ S ) (31)
over the normalized cavity surface S . Upon inspection, integral equation (31) turns out to be associated with the exterior elastostatic problem:

∇ξ •(C : ∇ξU ) = 0 ( ξ ∈ R 3 \ B), T = -σ(x o , ω)•n ( ξ ∈ S ). ( 32
)
where T = (C : ∇ξU )•n is the traction vector associated with U , and n is the normal on S outward to R 3 \ B. In fact, since the prescribed boundary traction T is defined in terms of a constant stress tensor, the solution to (32) can be conveniently recast as

U ( ξ) ≡ σ kℓ (x o , ω)U kℓ ( ξ) (33)
in terms of the individual solutions U kℓ = U ℓk to six canonical problems:

∇ξ •(C : ∇ξU kℓ ) = 0 ( ξ ∈ R 3 \ B), (C : ∇ξU kℓ )•n = - 1 2 n•(e k ⊗ e ℓ + e ℓ ⊗ e k ) ( ξ ∈ S ). (34)
which are independent of x o and ε.

Asymptotic form of equation [START_REF] Sokolowski | Introduction to shape optimization[END_REF]. By virtue of ( 23), the result (30), i.e. d = 1, implies that

ũε S N (ξ) = ε 3 |B| V (ξ) + o(1) , |V | = O(1) (ξ ∈ S N ) tε S D (ξ) = ε 3 |B| T (ξ) + o(1) , |T | = O(1) (ξ ∈ S D ) (35)
On the basis of ( 20), ( 21), (30), (33) and (35), the contributions of integrals over Γ ε in equation [START_REF] Sokolowski | Introduction to shape optimization[END_REF] reduce to

Γε ũε i (ξ, ω) tk i (ξ, x, ω; n) dΓ ξ + Γε t i (ξ, ω)û k i (ξ, x, ω) dΓ ξ = ε 3 |B| ρω 2 ûk (x o , x, ω)•u(x o , ω) -σk (x o , x, ω) : A : σ(x o , ω) + o(ε 3 ) (36)
where the constant fourth-order tensor A is defined by

A ijkℓ = 1 2µ I ijkℓ - ν 1 + ν δ ij δ kℓ - 1 |B| S U kℓ i ( ξ)n j ( ξ) dϑξ (37)
It should be noted that the tensor A and the associated canonical problems (34) are the same as in the previous study [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] dealing with infinite and semi-infinite media. For an arbitrarily shaped infinitesimal cavity, these canonical problems should in general be solved numerically. This is a modest computational task, since in fact one only needs to solve six elementary static problems which do not depend on x o . For the particular case of a spherical infinitesimal cavity, for which B is the unit ball, S is the unit sphere and |B| = 4π/3, problems (34) have an analytical solution, from which the following closed-form expression of A is obtained [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF]:

A = 3(1-ν) 2µ(7-5ν) 5 I sym 4 - 1+5ν 2(1+ν) I 2 ⊗ I 2 (38)
Recapitulation: governing problem for the scattered field in the limit as ε → 0. On the basis of the foregoing analysis, the leading contributions of the scattered field in the limit as ε → 0 are obtained as follows. The scattered displacement on the surface of the vanishing cavity, ũε Γε , is of the form (30) with U given by (33) in terms of the six solutions, U kℓ , to the canonical exterior elastostatic problems (34). As noted earlier, U kℓ depend solely on the void shape (as defined by S or B), and therefore need to be computed only once irrespective of the set of grid points x o or angular frequencies ω being considered. The constant fourth-order tensor A can then be determined on the basis of (37).

With the foregoing result (i.e. with ũε Γε pre-computed), the scattered field (ũ ε S N , tε S D ) on S is of the form (35), where (V , T ) are governed by the limiting form of the integral equation ( 18) as ε → 0. On the basis of (36), the latter BIE is obtained as

S V i (ξ, ω) -V i (x, ω) tk i (ξ, x; n) dΓ ξ + S V i (ξ, ω)[ tk i (ξ, x, ω; n)] 2 dΓ ξ - S T i (ξ, ω)û k i (ξ, x, ω; n) dΓ ξ = σk (x o , x, ω) : A : σ(x o , ω) -ρω 2 ûk (x o , x, ω)•u(x o , ω) (x ∈ S) (39)
The above formula and can be interpreted as an integral equation governing the elastodynamic state (V , T ) in the cavity-free, reference body Ω (with homogeneous boundary conditions on S) due to a point source at x o , generated by the superposition of a concentrated strain A : σ(x o , ω) and a concentrated force -ρω 2 u(x o , ω). The right-hand side of (39), and therefore its solution depends on the grid point x o , the angular frequency ω and the shape of the vanishing cavity.

Here it is useful to note that the small-cavity approximation (29) is inherently related to the lowfrequency approximation of elastic scattering in infinite media [START_REF] Dassios | Low frequency scattering[END_REF]. In particular, the latter theory, based on series expansions of field variables with respect to the elastic wave numbers, involves canonical elastostatic problems that are similar to (32).

Topological derivative by the direct approach

To evaluate the topological derivative T (x o , ω) as defined by [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF], one needs the O(ε 3 ) (leading) contribution of ũε on S obs ⊂ S. For the cavity identification problem using external boundary measurements, the topological derivative can be expressed directly from [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF] in terms of V as

T (x o , ω) = S obs Re ∂ϕ ∂u u(x), x •V (x) dΓ x . (40)

ADJOINT FIELD APPROACH

In what follows, an alternative formulation of the elastodynamic topological derivative for finite bodies, eq. ( 40), is developed by means of the adjoint field method. By maintaining the rigor of the earlier elastostatic developments [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF] while elevating their transparency, it will be shown that the application of the adjoint field technique to elastodynamics leads to an elegant formulation whose compactness is matched by its computational efficiency.

To arrive at an adjoint counterpart of (40), let ⋆ u(ξ) denote the elastodynamic displacement field in the finite, cavity-free reference body Ω due to prescribed surface traction ⋆ p(ξ) acting on its boundary. In mathematical terms this solution, herein termed the adjoint field, satisfies the field equations and mixed boundary conditions

∇•(C : ∇ ⋆ u) = -ρω 2 ⋆ u, ξ ∈ Ω, ⋆ t = ⋆ p, ξ ∈ S N , ⋆ u = 0, ξ ∈ S D . (41)
Here it should be noted that a general definition of the adjoint state permits an arbitrary distribution of body (i.e. interior) and surface sources, with the latter specialization adopted here to facilitate the formulation. With reference to ( 16) defining the scattered field ũε over the "punctured" domain Ω ε (bounded by S ∪ Γ ε ) due to primary excitation p, application of the Graffi's reciprocal theorem in elastodynamics [START_REF] Wheeler | Some theorems in classical elastodynamics[END_REF] in terms of ⋆ u and ũε over Ω ε yields

p(ξ) =    ∂ϕ ∂u (u(ξ), ξ), ξ ∈ Γ obs ⊂ S N 0, ξ ∈ S N \Γ obs ( 43
)
where ϕ is the featured misfit function, e.g. given by (5); ∂ϕ/∂u is defined through [START_REF] Éa | The shape and topological optimization connection[END_REF]; u is the free field in Ω due to primary excitation p (see [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF]); and, as examined earlier, the external boundary S of the actual body Ω true (and thus that of the cavity-free reference domain Ω) is assumed to contain the measurement surface so that Γ obs ⊂ S. In practical terms, the assumed distribution of ⋆ p over Γ obs can be interpreted as being proportional to a measure of the misfit between the measured displacement field (u obs ) and the (reference) free-field u. In particular, if ϕ is given by ( 5), ∂ϕ/∂u = W (uu obs ) so that for W = cI one has ⋆ p = c(u obs -u), where c is a suitable constant and u obsu is the observed scattered field in Ω true due to B true . On the basis of ( 13), ( 16) and (42), one finds that

J (Ω ε ; p) = J (Ω; p) -Re Γε ũε • ⋆ t dΓ - Γε ⋆ u • tε dΓ + o( ũε ) = J (Ω; p) -Re Γε ũε • ⋆ t dΓ + Γε ⋆ u •t dΓ + o( ũε ) as ε → 0. ( 44
)
For a given primary excitation p, the last integral over Γ ε in (44b) is known since the adjoint field ⋆ u is determined, through (43), solely in terms of u, u obs and the chosen cost function density ϕ. In fact, it can be shown by virtue of (10) that

Γε ⋆ u•t dΓ = - Bε ∇•( ⋆ u•σ) dΩ = Bε ρω 2 ⋆ u•u -∇ ⋆ u : σ dΩ = ε 3 |B| ρω 2 ⋆ u•u -∇ ⋆ u : σ (x o , ω) + o(ε 3 ) (45) = ε 3 |B| ρω 2 ⋆ u•u - 1 2µ ⋆ σ : σ - ν 1+ν tr ⋆ σ trσ (x o , ω) + o(ε 3 ) as ε → 0 where ⋆ σ = C : ∇ ⋆ u
is the stress tensor associated with the adjoint field, and the minus sign in front of the first integral over B ε appears because t is defined in terms of the unit normal n pointing to the interior of B ε .

To complete the expansion of J (Ω ε ; p) in (44b), one is left with evaluating the leading asymptotic behavior of the "bubble" integral Γε ũε • ⋆ t dΓ for vanishing ε, a task that can be effected in a straightforward manner with the aid of the results of Section 4. In particular, by employing (30) and (33), one obtains

Γε ũε • ⋆ t dΓ = ε 3 ⋆ σ ij (x o , ω) S U i ( ξ) n j ( ξ) dϑξ + o(ε 3 ) (46) = ε 3 ⋆ σ ij (x o , ω) σ kℓ (x o , ω) S U kℓ i ( ξ) n j ( ξ) dϑξ + o(ε 3 ) as ε → 0. ( 47
)
On substituting (47) and ( 45) into (44), the leading asymptotic behavior of J (Ω ε ; p) for vanishing ε can be further shown to take the form

J (Ω ε ; p) = J (Ω; p) + ε 3 |B| Re ⋆ σ : A : σ -ρω 2 ⋆ u•u (x o , ω) + o(ε 3 ) as ε → 0, (48)
where the constant fourth-order tensor A is again defined through (37), i.e. by

A ijkℓ = 1 2µ I ijkℓ - ν 1 + ν δ ij δ kℓ - 1 |B| S U kℓ i ( ξ)n j ( ξ) dϑξ
Finally, by virtue of ( 8) and ( 48), the adjoint-field formula for the topological derivative can be written in the compact form as

T (x o , ω) = Re ⋆ σ : A : σ -ρω 2 ⋆ u•u (x o , ω) (49)
where u and ⋆ u are the elastodynamic displacement fields solving the free-field problem [START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF] and the adjoint problem (41), respectively.

SHAPE SENSITIVITY APPROACH

The concept of topological sensitivity investigated herein is related to the asymptotic behavior of a cost function with respect to the introduction of a new, infinitesimal cavity. In contrast, the concept of shape sensitivity [START_REF] Sokolowski | Introduction to shape optimization[END_REF] is related to infinitesimal geometric perturbations of pre-existing finite cavities. Despite their obvious differences, however, these two subjects are related. In this section, a third method for computingT (x o , ω) is established by exploiting this relationship, providing additional insight into the concept of topological sensitivity. For that purpose, one notes that (8) implies

d dε J (Ω ε ; p) = 3ε 2 |B| T (x o , ω) + o(ε 2 ), as ε → 0. ( 50
)
and that the derivative dJ /dε can be expressed for any meaningful ε > 0 by means of a shape sensitivity formula [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF][START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF]:

d dε J (Ω ε ; p) = u(x o , ω).
Considerations similar to those presented in this Section have been developed in [START_REF] Éa | The shape and topological optimization connection[END_REF] and [START_REF] Novotny | Topological sensitivity analysis[END_REF] for scalar (potential) problems and spherical vanishing cavities. Equation ( 53) is a generalization of these results to vector (elastodynamic) problems and vanishing cavities of arbitrary shape. This formula can also be established by means of a direct calculation (i.e. without explicitly invoking a shape sensitivity argument), as shown in the Appendix. It should however be emphasized that the computational procedure based on exploiting the connection between shape and topological sensitivities, i.e. on definition (52) of T (x o , ω), is neither simpler nor more efficient than that developed in Section 5, contrary to what is suggested in [START_REF] Novotny | Topological sensitivity analysis[END_REF].

DISCUSSION

Computational issues

Solution methods. To establish expressions (40) and (49) for elastodynamic topological derivative, it was necessary to investigate the limiting behavior of the scattered field for vanishingly small cavities, a task for which the BIE formulation of the scattering problem is well suited. Similarly, the solution of the canonical exterior elastostatic problems (34) is best computed using a boundary element method.

However, the actual computation of the field T (x o , ω) can be achieved by means of any numerical solution technique available for elastodynamic problems, such as the finite element method (FEM), the finite difference method (FDM) or the boundary element method (BEM). This is especially true for the adjoint field approach of Section 5, where the adjoint problem is, like the forward problem, a conventional elastodynamic boundary-value problem. On the other hand, obtaining an accurate evaluation of T (x o , ω) by means of FEM or FDM using the direct approach of Section 4 may prove more difficult because the associated elastodynamic problems involve interior point sources (see eq. ( 39)).

The present work is a part of the ongoing study on defect identification from elastodynamic measurements, for which the BEM is a convenient numerical tool allowing easy re-meshing as the trial defect evolves. For that reason, the numerical experiments in Section 8 are performed by means of a boundary element method. In what follows, the BEM-based computational procedures for evaluating T (x o , ω) are described in some detail.

Elastodynamic BEM equations. On introducing a boundary element discretization of the surface S (together with the boundary displacements and tractions supported thereon) into integral equation [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundary integral equation[END_REF] and performing pertinent element quadrature [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Dominguez | Boundary elements in dynamics[END_REF], one arrives at a linear system of equations

[H(ω)]{X} = [G(ω)]{P } ( 54
)
where {X} and {P } are 'vectors' collecting respectively the N unknown degrees of freedom (DOFs) and prescribed nodal tractions on S N ; [H(ω)] and [G(ω)] are fully-populated matrices. The unknown DOFs are thus found by solving (54) for {X}, with the matrix [H(ω)], of size N × N , being invertible (based on the assumption that ω does not define an eigenfrequency for the reference body) and nonsymmetric. Likewise, the integral representation formulas for the displacement vector and the stress tensor at a grid point x o have respectively the forms

u(x o , ω) = [G u (x o , ω)]{P } -[H u (x o , ω)]{X} σ(x o , ω) = [G σ (x o , ω)]{P } -[H σ (x o , ω)]{X} (55)
where the influence matrices

[G u (x o , ω)], [G σ (x o , ω)], [H u (x o , ω)]
and [H σ (x o , ω)] result from the BEM discretization of the relevant boundary integrals. Here, it is assumed that the set of internal grid points x o , over which T (x o , ω) is to be computed, is chosen beforehand. One may note that solving the forward problem (54) and evaluating the representation formulas (55) for N grid grid points via traditional methods entail O(N 3 ) and O(N × N grid ) computation times, respectively. However, both stages can be considerably accelerated for large N and N grid by resorting to fast solution strategies such as the fast multipole method [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF].

BEM-based computational procedures for the topological derivative. Both the direct and the adjoint field approaches entail the following common operations:

1. Solve the BEM matrix equation ( 54) for the free-field DOFs {X}; 2. Compute integral representations (55) of the free-field for all grid points x o ; 3. Compute the 6 × 6 constant tensor A defined by (37), which in general entails solving the set of six canonical problems (34). For spherical vanishing cavities, A has been precomputed analytically and is given by (38).

Additional operations required by the direct approach include:

4D. On the basis of Steps 2 and 3, compute the right-hand side {F (x o , ω)} of (39) for all grid points; 5D. Solve problem (39) in discretized form, i.e. Multiple excitations. A generalization of the cost functional (4) and its topological derivative as given by ( 40), ( 49) or (53) to multiple excitations p q (q = 1, 2, . . . Q) is straightforward and involves external summation in the form of

J (Ω c ) ≡ Q q=1 J (Ω c ; p q ) and T (x o , ω) ≡ Q q=1 T (x o , ω; p q ).
(56)

Procedural similarities between topological and shape sensitivity formulations

In addition to the previously emphasized mathematical link, there are also strong similarities between the computational procedures devised for topological and shape sensitivity analyses. In both situations, the direct approach involves finding the point-wise shape (resp. topological) sensitivity of the primary field variable, which entails solving one sensitivity boundary value problem per geometrical parameter (resp. per grid point). In contrast, the adjoint field approach requires solving but one new boundary value problem and provides the most direct route for computing the sensitivity of the objective function at hand, although it does not furnish sensitivity information on the primary field variable. There is, however, a noteworthy difference between the direct approaches for topological and shape sensitivity analyses. The sensitivity boundary value problems associated with topological sensitivity analysis involve singular loadings (in the form of point sources), whereas the latter are typically nonsingular in shape sensitivity analysis.

Explicit expression of T (x o , ω) in terms of Green's tensors

As examined earlier, alternative expressions (40), (49) and (53) for T (x o , ω) are implicit in that they rely upon solutions of boundary-value problems on the cavity-free reference body Ω.

Explicit expressions for T (x o , ω) are obtained when the elastodynamic fundamental solution ûk i (ξ, x, ω), tk i (ξ, x, ω; n) is a Green's tensor for Ω, i.e. when it satisfies the boundary conditions

ûk i (ξ, x, ω) = 0 (ξ ∈ S D , ξ = x), tk i (ξ, x, ω; n) = 0 (ξ ∈ S N , ξ = x)
In this case, the free and adjoint displacement fields are given by

u k (x) = S N p i (ξ)û k i (ξ, x, ω) dΓ ξ , ⋆ u k (x) = S N ⋆ p i (ξ)û k i (ξ, x, ω) dΓ ξ
and the expression (49) becomes explicit as well. This constitutes a formal generalization of the results of [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], where the topological derivative for infinite and semi-infinite elastic media was considered. However, such Green's tensors are in general not known explicitly for finite domains.

Relationship to the Eshelby solution

The well-known Eshelby solution, i.e. the elastostatic response of a three-dimensional infinite body Ω = R 3 to a uniform initial strain distribution ε I prescribed over an ellipsoidal region B ⊂ Ω, is such that the resulting (total) strain ε E is uniform within B, i.e.

ε E (ξ) = S : ε I (ξ ∈ B), (57)
where the constant fourth-order Eshelby tensor S is known analytically [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]. As a result, the induced stress field, σ E , is also uniform within B, with σ E = C : (S -I sym 4 ) : ε I . Since the tractions must be continuous across the boundary Γ = ∂B of the "loaded" region, the Eshelby displacement solution in R 3 \ B coincides with the solution of the elastostatic exterior problem when the cavity surface Γ is subjected to tractions σ E • n, where n is the normal on Γ outward to R 3 \ B. It follows that, when the shape of the unit cavity B is ellipsoidal, the solution of the elastostatic exterior problem (32) coincides with that of the Eshelby inclusion problem with the initial strain ε I taken so that σ(x o , ω) = -C : (S -I sym 4 ) : ε I . On the basis of (35) and earlier discussion (immediately following eq. ( 39)), the free-field term σk (x o , x, ω):A : σ(x o , ω) appearing in the right-hand side of (39) is thus expected to coincide with the limiting form (as |B| = ε 3 |B| → 0) of the elastodynamic displacement field u E (ξ, ω)/(ε 3 |B|) generated in R 3 \ B by the Eshelby initial strain ε I =-(S -I sym 4 ) -1 : C -1 : σ(x o , ω). Owing to the fact that the latter displacement field is given by

u E k (ξ, ω) = B ε I (ξ) : σk (ξ, x, ω) dΩ ξ = -(S -I sym 4 ) -1 : C -1 : σ(x o , ω) : B σk (ξ, x, ω) dΩ ξ (58) = ε 3 |B| σk (x o , x, ω) : -(S -I sym 4 ) -1 : C -1 : σ(x o , ω) + o(ε 3 ) as ε → 0. the implied coincidence holds only if A = [(I sym 4 -S ) -1 : C -1 ].
For a spherical cavity, this equality is easily checked by direct inspection, using (38) and the known expression of S , i.e.

S = 1 15(1 -ν) (8 -10ν)I sym 4 -(1 -5ν)I 2 ⊗ I 2 . ( 59
)
According to the foregoing reasoning, the closed form of tensor A can be obtained in the same manner from the known value of S for any cavity of ellipsoidal shape. Moreover, on considering the limiting cases of infinitely thin ellipsoidal cavities, one should be able to establish topological sensitivity formulae for circular or elliptical infinitesimal cracks. To provide a basis for comparison, numerical values of the topological derivative T (r o , ω) are computed by means of a three-dimensional boundary element solution [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF] applied to both direct approach (40) and adjoint-field formula (49) with p = µ, ν = 0.3, R/a true = 3 and ωR/c s = 3 where c s = µ/ρ denotes the shear wave speed in the solid. In the numerical model, the inner and outer surfaces of the annular sphere were discretized using 486 and 150 eight-noded boundary elements, respectively. For the excitation frequency chosen, this mesh provides no less than 12 element lengths per shear wavelength. A comparison between the 3D numerical results and the exact solution given by ( 72) is shown in Fig. 2. As can be seen from the display, both numerical techniques for computing T (r o , ω) yield results that are in reasonable agreement with the closed-form solution. In the Figure, the RMS errors for the featured set of grid points are 2.6 10 -3 for the direct approach and 1.5 10 -3 for the adjoint-field method, respectively. It should be mentioned, however, that the accuracy of numerical estimates of T (r o , ω) was found to deteriorate for values of r o /R close to unity, primarily as a result of the near-singular nature of the integral representation for stresses at observation points close to the boundary.

Cavity embedded in a cube

To illustrate the utility of topological derivative as a preliminary tool for elastic-wave sounding of finite bodies, the next example deals with the delineation of a spherical cavity of radius a true = 0.5d, hidden inside a solid cube of size 6d × 6d × 6d. Similar to the previous example, the elastic solid is characterized by the shear modulus µ, mass density ρ, and Poisson's ratio ν = 0.3. With reference to a Cartesian coordinate system aligned with the box edges (see Fig. 3), the cube and the cavity are centered respectively at (0, 0, 0) and (d, 1.5d, d). The cubical body, with external surface S, is fixed (u(ξ, ω) = 0, ξ ∈ S D ) over the bottom patch S D = {ξ ∈ S| -2.4d<ξ 1 <2.4d, -2.4d<ξ 2 <2.4d, ξ 3 = -3d}. Prior to the application of time-harmonic excitation used to illuminate the cavity, the rest of the external surface, S N = S \ S D , is traction free. As indicated in Fig. 3, surfaces of the cube and the cavity are discretized using 600 and 64 eight-node (quadratic) boundary elements, respectively. In sequence, a virtual elastodynamic experiment is performed on each of the four vertical faces of the cube. With reference to the vertical testing face ξ 1 = -3d, the experiment parameters are chosen Table I. Cavity embedded in a cube: comparison of direct and adjoint method for the evaluation of T (x o , ω). so that i) Γ obs = {ξ ∈ S N , |ξ 1 = -3d, -3d<ξ 2 <3d, -3d<ξ 3 <3d}, and ii) the cavity is illuminated in sequence by five localized (pyramid-shaped) distributions p q (ξ) (ξ ∈ Γ obs , q = 1, 2, . . . 5) of surface tractions centered at (-3d, 0, 0) and (-3d, ±1.8d, ±1.8d), respectively, applied over square patches of four boundary elements. In the context of the employed boundary element mesh, these localized tractions, each of resultant P = µd 2 and acting in the positive ξ 1 -direction, are generated by setting the nodal traction at the respective central node of each square patch (e.g. at (-3d, 1.8d, 1.8d)) to a value of 25µ/9 while assigning half that value to the four closest (mid-side) boundary nodes. Virtual experiments on the remaining vertical faces, i.e. ξ 1 = 3d, ξ 2 = -3d and ξ 2 = 3d, are performed in an analogous fashion with the applied normal tractions (pressure) acting respectively in the negative ξ 1 -, positive ξ 2 -, and negative ξ 2 -direction, respectively. With reference to (3) and (56), the cost function J and its topological derivative T are calculated with W = I 2 and by summing the contribution from all Q = 4 × 5 = 20 experiments. The topological derivative T (x o , ω) is evaluated at the nodes of a regular cubic grid of 11×11×11 = 1331 interior points such that -2.5d ≤ x o i ≤ 2.5d (i = 1, 2, 3). As mentioned earlier, all elastodynamic calculations are performed using the boundary element analysis in [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF].

To check the computational developments in the absence of suitable analytical solutions, Table I shows a comparison between the values of topological derivative for the cavitated cube problem computed by the direct differentiation approach (40) and the adjoint field approach (49) at two representative frequencies, ω = 1 and ω = 2, where ω = ω d ρ/µ ≡ ωd/c s , with c s denoting the shear wave speed. As can be seen from the display, relative error between the two estimates does not exceed 1.5 percent despite significant differences in the respective computational schemes.

In view of the findings in structural shape optimization [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF][START_REF] Garreau | The Topological Asymptotic for PDE Systems: The Elasticity Case[END_REF], the elastodynamic topological derivative developed in this study is expected to facilitate the identification of internal defects (voids) in finite solids through an algorithm that would entail i) evaluation of T (x o , ω) over a suitable grid of sampling points spanning the reference (cavity-free) body Ω, and ii) delineation of possible cavities through regions in Ω where T takes the largest negative values. To investigate such possibility, Figs. 4 and 5 illustrate the distribution of topological derivative (computed via the adjoint field approach) corresponding to excitation frequencies ω = 1 and ω = 2, respectively. In both displays, the distribution is plotted in three coordinate (cutting) planes containing the center of the true cavity. To provide a reference, intersection of each plane with the true cavity surface is outlined in white. One may observe that each distribution indeed points toward the true cavity through negative values of T (x o , ω). Unfortunately, the "low-frequency" map suffers from blurring, while an increased resolution of the "higher-frequency" map is diminished by the appearance of spurious minima. One should mention, however, that both excitation frequencies fall into the so-called resonance region [START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF] where the wave lengths are larger than the size of the scatterer. In particular, (λ s /D)| ω=1 = 2π and (λ s /D)| ω=2 = π, where D is the cavity diameter and λ s denotes the shear wave length. As examined in [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], the particular effectiveness of long wavelengths for preliminary imaging by way of topological derivative is not surprising, since the assumption of an infinitesimal cavity, implicit to [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF], is better conformed with by finite cavities that are 'small' relative to the probing wavelength.

With diagrams such as those in Figs. 4 and5, an algorithm for identifying plausible cavity locations could be thus devised on the basis of the non-zero distribution of an auxiliary function

T (x o ) = T (x o , ω), T < C, 0, T ≥ C, (60)
where C < 0 denotes a suitable threshold value. With such definition, it is also possible to combine a) b) c) c) the individual advantages of different probing wavelengths by employing the product of (60) at several frequencies [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF]. As an illustration of the latter approach, Fig. 6

ξ 1 /d ξ 1 /d ξ 2 /d ξ 2 /d ξ 3 /d ξ 2 /d ξ 3 /d
ξ 1 /d ξ 1 /d ξ 2 /d

SUMMARY

In this study, the concept of topological derivative that has its origins in elastostatics and shape optimization is extended to 3D elastic-wave imaging of finite solids. On taking the limiting form of the boundary integral equation governing the scattered field caused by a cavity with diminishing size, the topological derivative, which quantifies the sensitivity of the featured cost functional due to the creation of an infinitesimal hole, is formulated (in terms of an elastodynamic field for the reference, void-free body) using three alternative methodologies, namely i) the direct differentiation approach, ii) adjoint field method, and iii) limiting form of the shape sensitivity analysis. The proposed derivation of the topological gradient has two noteworthy features. Firstly, it is organized so as to ensure adequate rigor while avoiding the use of non-essential mathematical concepts that would put it beyond the reach of most researchers and practitioners in engineering sciences. Secondly, it follows a pattern which is generic, i.e. transposable to i) other numerical techniques (e.g. BIE, BEM, or FDM), ii) variety of physical problems described by linear field equations such as acoustics, heat transfer, and electromagnetics, and iii) other types of infinitely small objects, such as inhomogeneities or cracks.

It is also shown that the cavity nucleation, central to topological derivative, is closely related to the Eshelby solution for the response of an infinite elastic solid due to a localized application of uniform initial strain.

For the example of a pressurized spherical shell that admits a closed-form solution, both the direct and the adjoint field estimates of topological derivative, obtained using a three-dimensional BIE method, are found to agree well with the exact result. To examine the potential of topological sensitivity for preliminary elastic-wave sounding of finite bodies, a numerical example for the solid cube hiding a spherical cavity is also presented. The results suggest the utility of the proposed imaging algorithm, wherein a plausible void region is delineated through negative values of the topological derivative at internal (sampling) points. It should be noted, however, that the proposed imaging approach is essentially of heuristic nature since the topological sensitivity is by definition related to infinitesimal cavities, whereas the obstacles being sought are of finite size. It is also found that the use of multiple excitation frequencies as a tool to illuminate the cavity may enrich the experimental data set and thus elevate the performance of topological derivative as a preparatory imaging tool. Next, with the help of the equality where u and σ are given by ( 65) and (66), (67), respectively. Therefore, the topological derivative T (r o , ω) at the grid point x o = r o er is given by: T (r o , ω) = -1 p u(R, ω) -u obs (ω) γ 2 4µ(3γ 2 -1) (σrr + 2σϕϕ) 2 (r o , ω)

+ γ 2 4µ(9γ 2 -1) 2 3 (σrr -σϕϕ) 2 (r o , ω) -4µγ 2 Q 2 u 2 (r o , ω) (72)
One may note that expression (72) evaluated for r o = 0 coincides with the result (71) obtained by the direct approach.

  [H]{V } = {F (x o , ω)}, for all grid points x o ; 6D. Compute T (x o , ω) for all grid points x o using (40). Additional operations required by the adjoint field approach include: 4A. Compute the adjoint nodal tractions { ⋆ P } using (43) together with experimental observations and free-field displacements computed along Γ obs ∈ S in Step 1; compute the adjoint solution { ⋆ X} by solving (54) with {P } replaced with { ⋆ P }; 5A. Compute the adjoint displacements and stresses, ⋆ u(x o , ω) and ⋆ σ(x o , ω), for all internal grid points x o from (55) with {X} and {P } replaced by { ⋆ X} and { ⋆ P }; 6A. Compute T (x o , ω) for all grid points x o using (49).

Figure 1 .Figure 2 .

 12 Figure 1. Pressurized spherical shell: geometry and notation.

Figure 3 .

 3 Figure 3. Cavity embedded in a solid cube: boundary element mesh.

ω

  x o /d T (x o , ω)/d T (x o , ω)/d Relative (5, -0.5, -2.5) +1.9531 e-3 +1.9526 e-3 2.56 e-4 (0.5, -0.5, -2.0) -3.0256 e-3 -3.0156 e-3 -3.31 e-3 (0.5, -0.5, -1.5) -6.2075 e-3 -6.1476 e-3 -9.70 e-3

Figure 4 .

 4 Figure 4. Distribution of T (x o , ω)/d in coordinate planes: a) ξ3 = d, b), ξ2 = 1.5d, and c) ξ1 = d containing the center of the true cavity (ω = 1).

Figure 5 .

 5 Figure 5. Distribution of T (x o , ω)/d in coordinate planes: a) ξ3 = d, b), ξ2 = 1.5d, and c) ξ1 = d containing the center of the true cavity (ω = 2).

  plots the distribution of the product T (x o , c s /d) × T (x o , 2c s /d) in three coordinate planes containing the center of the true cavity with C = 0.5 inf x o T (x o , ω). Notwithstanding the limited accuracy characterizing the individual maps for ω = 1 and ω = 2, the combined result stemming from (60) clearly points to a single cavity with its center and size closely approximating the true void configuration.

Figure 6 .

 6 Figure 6. Distribution of (1/d 2 ) T |ω=1 × T |ω=2 in coordinate planes: a) ξ3 = d, b), ξ2 = 1.5d, and c) ξ1 = d containing the center of the true cavity.

  and to MTS Systems Corporation for providing the opportunity for M. Bonnet to visit the University of Minnesota through the MTS Visiting Professorship of Geomechanics. First, due to the Neumann boundary condition in (32) and the local equilibrium equation (34a), one has σij(x o , ω) S ⋆ U inj dϑξ = -S ⋆ U iΣij nj dϑξ =ξ) = C : ∇ξU and O = R 3 \ B. In (62), both ξ 2 ∇ξU and ξ 2 ∇ξ ⋆ U are bounded as a consequence of the exterior problem (34), so that the integral over the normalized exterior domain O is convergent.

Copyright c 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2004; 61:2344-2373 Prepared using nmeauth.cls

Γε ( ⋆ u• tεũε • ⋆ t) dΓ = S N ũε • ⋆ p dΓ,(42)which makes use of the fact that u = ⋆ u ≡ 0 on S D , tε ≡ 0 on S N and ⋆ t = ⋆ p on S N . To utilize (42) as a tool for obtaining the formula for topological derivative, it is further convenient to specialize the surface traction ⋆ p generating the adjoint state ⋆ u according to ⋆

Γερω 2 u ε • ⋆ u εσ ε : ∇ ⋆ u ε (ξ)θ n (ξ) dΓ ξ(51)where ⋆ u ε denotes the adjoint displacement field for the cavitated body, which is governed by equations (2) with the load p on S N replaced by the adjoint load ⋆ p defined by (43). In addition, the normal transformation velocity θ n (ξ) must be chosen so that the mapping ξ → ξ + θ n (ξ)n(ξ)δε transforms Γ ε into Γ ε+δε (neglecting o(δε) contributions), which is achieved with the choice θ n = ξ•n. From (51) and (50), an alternative expression of the topological derivative is then inferred asT (x o , ω) = lim ε→0 1 3ε 2 |B| Γε ρω 2 u ε • ⋆ u εσ ε : ∇ ⋆ u ε (ξ)θ n (ξ) dΓ ξ(52)On using notations and results of Section 4, the leading contribution of the above integral for vanishing ε is readily established, leading to the more explicit formulaT (x o , ω) = 1 3 |B| S ρω 2 u(x o , ω)• ⋆ u(x o , ω) σ(x o , ω) + Σ( ξ) : ∇ ξ ⋆ u(x o , ω) + ∇ξ ⋆ U ( ξ) θ n ( ξ) dϑξ (53)where Σ( ξ) = C : ∇ξU ( ξ), U ( ξ) is given by (33), and ⋆ U ( ξ) is the normalized elastostatic displacement given by (33) with σ(x o , ω) replaced by ⋆ σ(x o , ω) = C : ∇ ξ ⋆

3Σij ⋆ U i,j = ⋆ Σij Ui + Σij ⋆ U i -⋆ Σij U i,k + Σij ⋆ U i,k ξk + Σ ik ⋆ U i,k ξj ,jwhich holds for any pair of elastostatic displacement fields U and ⋆ U without body forces, the divergence formula is applied to the domain integral of (62), leading toσij (x o , ω) S ⋆ U inj dϑξ = 1 3 σij (x o , ω) S ( ⋆ U i -⋆ U i,k ξk )nj dϑξ + 1 3 ⋆ σij(x o , ω) S (Ui -U i,k ξk )nj dϑξ -1 3 S (Σ ik ⋆ U i,k )( ξjnj ) dϑξ (63)where the boundary conditionsΣ•n = -σ(x o , ω)•n and ⋆ Σ•n = -⋆ σ(x o , ω)•n have also been used. Then, using the identity(Vi -V i,k ξk )nj = (Vi ξk ),j n k -(Vi ξk ) ,k nj -Vi,j ξk n k + 3Vinj with V = U or V = ⋆U and invoking the following variant of the Stokes formulaS (nj f ,k -n k f,j ) dϑξ = 0 j, k = 1, 2, 3which holds true for any piecewise closed surface S and any continuous and piecewise differentiable function f ( ξ), equation (63) becomes, after rearranging terms:⋆ σij(x o , ω) S Uinj dϑξ = 1 3 S ⋆ σij (x o , ω)Ui,j + σij (x o , ω) ⋆ Ui,j + Σij ⋆ U i,j ( ξk n k ) dϑξ(64)Finally, on substituting (64) and the identityσ : ∇ ξ ⋆ u -ρω 2 u• ⋆ u (x o , ω) = 1 3 |B| S ρω 2 u• ⋆ uσ : ∇ ξ ⋆ u (x o , ω) ξ•n dϑξ into (61), the alternative expression (53) for the topological derivative is established.II. EXACT SOLUTION FOR THE SPHERICAL SHELL PROBLEMIn this example, the reference body Ω is a spherical ball of radius R. Let Ωa = Ω \ Ba denote the spherical shell bounded by concentric spheres of radii a and R > a. The external surface r = R is subjected to a uniform pressure, i.e. p = -per for r = R (where (r, ϕ, θ) denote spherical coordinates and er is the unit radial vector), while the cavity surface, r = a, is traction-free. The resulting elastodynamic state (u a , σ a ) in the spherical shell is radially symmetric, so that u a (ξ, ω) = u a (r, ω)er and σ a (ξ, ω) = σ a (r, ω)er ⊗ er + τ a (r, ω)[eϕ ⊗ eϕ + e θ ⊗ e θ ]. Copyright c 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2004; 61:2344-2373 Prepared using nmeauth.cls

ACKNOWLEDGEMENTS Special thanks are extended to École Polytechnique for providing a Visiting Scientist appointment to B. Guzina, Copyright c 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2004; 61:2344-2373 Prepared using nmeauth.cls

APPENDIX I. PROOF OF EXPRESSION (53) FOR T (x o , ω)

From equations (34), (37) and (49), one has

The integral over S can be transformed via integration by parts, noting that both U and ⋆ U denote elastostatic states without body forces.

Free field. For the foregoing configuration, the nonzero components of the free field (a = 0) can be calculated as

(1 -γ 2 Q 2 y 2 ) sin(Qy) -Qy cos(Qy) sin(Qy) -Qy cos(Qy) (65)

where

Expansion of the displacement for a vanishingly small cavity. The displacement on the external surface for the case of a finite cavity of radius a, u a (R, ω), is given by

where x = a/R and

For a sufficiently small cavity radius, denoted by a = ε for consistency, a Taylor expansion of (68

where u(R, ω) = u 0 (R, ω), i.e. the free field on S, is obtained by setting x = 0 in (68) and

Topological derivative at the origin by direct differentiation approach. Let (u true , σ true ) correspond to the 'true' cavity B true defined by the radius a true . The observation surface S obs is taken as the whole external surface S, so that u obs (ω) ≡ u obs (ω)er = u true (R, ω)er, and the cost function J(a) is defined by:

where the convention Ω0 = Ω is used. The topological derivative of J at the origin, T (0, ω), is given by