# Topological derivative applied to cavity identification from elastodynamic data

Marc Bonnet, Bojan B Guzina

<figure><img src="image_1.png" /></figure>

To cite this version:

Marc Bonnet, Bojan B Guzina. Topological derivative applied to cavity identification from elastodynamic data. ECCOMAS, 2004, Paris, France. pp.355-362, 10.3166/reef.13.425-436. hal-00116389

HAL Id: hal-00116389

https://hal.science/hal-00116389v1

Submitted on 22 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Topological derivative applied to cavity identification from elastodynamic surface measurements

Marc Bonnet* — Bojan B. Guzina**

* Laboratoire de Mécanique des Solides (UMR CNRS 7649)
Ecole Polytechnique, F-91128 Palaiseau cedex
bonnet@lms.polytechnique.fr

** Department of Civil Engineering, University of Minnesota
500 Pillsbury Drive, Minneapolis, MN 55455-0116, USA
guzina@wave.ce.umn.edu

**ABSTRACT.** This article is concerned with the use of topological derivative as a tool for preliminary elastic-wave probing of bounded or unbounded solids for buried objects. A formulation for computing the topological derivative field, based on an adjoint solution, is presented. A set of numerical results is included to illustrate the utility of topological derivative for outlining the cavity location and size prior to doing an actual inversion of measurements. The results presented here were obtained from a BIE solution, but the proposed methodology is applicable to other computational platforms such as the finite element method.

**RÉSUMÉ.** Cet article porte sur l'application de la notion de gradient topologique à l'identification d'objets enfouis à l'aide de données élastodynamiques. On présente une formulation, reposant sur un état adjoint, pour le calcul du champ de dérivée topologique. Des exemples numériques sont présentés dans le but d'illustrer l'utilité de la dérivée topologique pour l'estimation de l'emplacement et de la taille d'une cavité, notamment dans le but d'obtenir une estimation initiale préliminaire à une procédure d'inversion. Les résultats ont été obtenus au moyen d'un programme fondé sur les éléments de frontière, mais le calcul de la dérivée topologique peut être effectué via toute autre technique, en particulier les éléments finis.

**KEYWORDS:** cavity identification, elastodynamics, inverse scattering, adjoint field method, shape sensitivity, topological derivative.

**MOTS-CLÉS :** identification de cavité, élastodynamique, diffraction inverse, méthode de l'état adjoint, dérivée géométrique, dérivée topologique.
# 1. Introduction

Three-dimensional imaging of cavities embedded in a solid using elastic waves is a topic of intrinsic interest in a number of applications ranging from nondestructive material testing to oil prospecting and underground object detection. In situations when detailed mapping of objects (defense facilities, buried waste) is required and only a few measurements can be made, the use of surface discretization-based boundary integral equation (BIE) techniques provides the most direct link between the surface measurements and the buried geometrical objects. Such an approach is well established for acoustic problems [COL92]. More recently, a BIE-based analytical and computational framework for the identification of cavities in a semi-infinite solid from surface elastodynamic measurements has been developed and reported in [GUZ03, NIN03]. Computing time considerations are such that global search methods are currently impractical for three-dimensional elastodynamic inverse scattering problems, and hence led to using a conventional gradient-based ‘blind’ optimization scheme. As often with such schemes, the highly non-convex character of the cost function makes the identification result dependent on the initial guess.

This has prompted the authors to investigate the use of topological derivative as a tool for preliminary probing. The formulation presented in this article is applicable to bounded or unbounded elastic media. A set of numerical results is included to illustrate the utility of topological derivative for outlining the cavity location, shape and size prior to doing an actual inversion of measurements. Despite the fact that the results presented here were obtained from a BIE solution, the proposed methodology is readily applicable to other computational platforms such as the finite element method.

# 2. Preliminaries

Let $\Omega^{\text{true}}$ denote an elastic body bounded by the external surface $S$, divided into complementary subsets $S_N$ and $S_D$ supporting prescribed tractions and displacements, respectively. An unknown cavity (or a set thereof) $B^{\text{true}}$ bounded by the closed surface(s) $\Gamma^{\text{true}}$ is embedded in $\Omega^{\text{true}}$, so that $\Omega^{\text{true}} = \Omega \setminus \bar{B}^{\text{true}}$ where $\Omega$ is the reference, i.e. cavity-free, counterpart of $\Omega^{\text{true}}$. On applying a steady-state traction $p$ on $S_N$ with angular frequency $\omega$, an elastodynamic state $u^{\text{true}}$ arises which solves the problem

$$ \mathcal{L}(\omega)u^{\text{true}} = 0 \quad (\text{in } \Omega^{\text{true}}), \\ t^{\text{true}} = p \quad (\text{on } S_N), \quad t^{\text{true}} = 0 \quad (\text{on } \Gamma^{\text{true}}), \quad u^{\text{true}} = 0 \quad (\text{on } S_D). \qquad [1] $$

Here $\rho$ is the mass density, $C$ is the fourth-order isotropic elasticity tensor characterized by the shear modulus $\mu$ and Poisson's ratio $\nu$, $\mathcal{L}(\omega)$ is the Navier linear partial differential operator, i.e. $\mathcal{L}(\omega)w = \text{div}(C: \nabla w) + \rho\omega^2 w$, and $t^{\text{true}} = \sigma^{\text{true}} \cdot n = (C: \nabla u^{\text{true}}) \cdot n$ denotes the traction vector associated with the displacement $u^{\text{true}}$ through Hooke's law. For simplicity, it is assumed that $\omega$ is not an eigenfrequency of any of the boundary-value problems appearing in the ensuing developments.

For the inverse problem of interest, where the location, topology and geometry of $B^{\text{true}}$ (or equivalently $\Gamma^{\text{true}}$) is being sought, the trace of $u^{\text{true}}$ on $S$ (denoted hereafter
as $\mathbf{u}^{\text{obs}}$) is assumed to be available over the measurement region $S^{\text{obs}} \subset S_{\text{N}}$. Let $\mathbf{u}^c$ denote the solution to the forward problem for a given excitation $\mathbf{p}$ and a trial cavity $B_c$ bounded by $\Gamma$: $\mathbf{u}^c$ is then defined over $\Omega_c = \Omega \setminus \bar{B}_c$, and governed by the equations

$$ \mathcal{L}(\omega)\mathbf{u}^c = \mathbf{0} \quad (\text{in } \Omega_c), \quad \mathbf{t}^c = \mathbf{p} \quad (\text{on } S_{\text{N}}), \quad \mathbf{t}^c = \mathbf{0} \quad (\text{on } \Gamma), \quad \mathbf{u}^c = \mathbf{0} \quad (\text{on } S_{\text{D}}). \quad [2] $$

where $\mathbf{t}^c$ is the traction vector associated with $\mathbf{u}^c$. To solve the inverse problem, a misfit cost function is set up in order to minimize the difference between $\mathbf{u}^{\text{obs}}$ and $\mathbf{u}^c$. Generic cost functions having the format

$$ \mathcal{J}(\Omega_c) = \int_{S^{\text{obs}}} \varphi(\mathbf{u}^c(\xi), \xi) \, d\Gamma_{\xi} \quad [3] $$

are considered. The output least-squares cost function, commonly used for such purpose, corresponds to the particular case where $2\varphi(\mathbf{w}, \xi) = (\overline{\mathbf{w}(\xi)} - \mathbf{u}^{\text{obs}}(\xi)) \cdot \mathbf{W}(\xi) \cdot (\mathbf{w}(\xi) - \mathbf{u}^{\text{obs}}(\xi))$, where $\mathbf{W}(\xi)$ is a $3 \times 3$ matrix-valued weighting function, assumed to be symmetric and positive definite (the simplest choice being $\mathbf{W}(\xi) = \mathbf{I}_2$), while overbar denotes complex conjugation.

### 3. Topological derivative

To aid the gradient-based minimization of $\mathcal{J}(\Omega)$, a tool often used for identifying $B^{\text{true}}$ on the basis of $\mathbf{u}^{\text{obs}}$, the development of topological derivative for the cost functionals of form [3], which would facilitate a rational selection of the necessary initial “guess” in terms of the location, topology and geometry of $B^{\text{true}}$, is investigated. To this end, let $B_{\varepsilon}(\mathbf{x}^o) = \mathbf{x}^o + \varepsilon\mathcal{B}$, where $\mathcal{B} \subset \mathbb{R}^3$ is a fixed bounded open set with boundary $\mathcal{S}$ and volume $|\mathcal{B}|$ containing the origin, define the region of space occupied by a cavity of (small) size $\varepsilon > 0$ containing a fixed sampling point $\mathbf{x}^o$. Following [GAR01, SOK99], one is in particular interested in the asymptotic behavior of $\mathcal{J}(\Omega_{\varepsilon})$ for infinitesimal $\varepsilon > 0$, where $\Omega_{\varepsilon} = \Omega \setminus \overline{B_{\varepsilon}(\mathbf{x}^o)}$, and $\overline{B_{\varepsilon}(\mathbf{x}^o)}$ is the closure of $B_{\varepsilon}(\mathbf{x}^o)$. The topological derivative $\mathcal{T}(\mathbf{x}^o)$ of the cost functional $\mathcal{J}(\Omega)$ at $\mathbf{x}^o$ for a cavity-free body is hence defined through an expansion of the form:

$$ \mathcal{J}(\Omega_{\varepsilon}) = \mathcal{J}(\Omega) + f(\varepsilon) |\mathcal{B}| \mathcal{T}(\mathbf{x}^o) + o(f(\varepsilon)) \quad (\varepsilon \ll \text{Diam}(\Omega), \ B_{\varepsilon}(\mathbf{x}^o) \subset \Omega) \quad [4] $$

where the function $f$ defines the leading asymptotic behavior of $\mathcal{J}(\Omega_{\varepsilon})$ and is such that $f(\varepsilon) \to 0$ as $\varepsilon \to 0$. The definition [4] is not restricted to spherical infinitesimal cavities (for which $\mathcal{B}$ is the unit ball, $\mathcal{S}$ the unit sphere and $|\mathcal{B}| = 4\pi/3$). In general, the value $\mathcal{T}(\mathbf{x}^o)$ is expected to depend on the shape of $\mathcal{B}$.

With reference to [4], the evaluation of $\mathcal{J}(\Omega_{\varepsilon})$ requires the knowledge of the elastodynamic solution $\mathbf{u}^{\varepsilon}$ to the forward problem [2] with $B_c$ replaced by $B_{\varepsilon} \equiv B_{\varepsilon}(\mathbf{x}^o)$. It is thus convenient to decompose the total displacement field $\mathbf{u}^{\varepsilon}$ as $\mathbf{u}^{\varepsilon} = \mathbf{u} + \tilde{\mathbf{u}}^{\varepsilon}$, where $\tilde{\mathbf{u}}^{\varepsilon}$ denotes the *scattered field* and $\mathbf{u}$ is the *free field* defined as the response of the void-free (reference) solid $\Omega$ due to given excitation (i.e. boundary traction) $\mathbf{p}$, so that

$$ \mathcal{L}(\omega)\mathbf{u} = \mathbf{0} \quad (\text{in } \Omega), \quad \mathbf{t} = \mathbf{p} \quad (\text{on } S_{\text{N}}), \quad \mathbf{u} = \mathbf{0} \quad (\text{on } S_{\text{D}}), \quad [5] $$
where $t$ is the traction vector associated with $u$, and

$$ \mathcal{L}(\omega)\tilde{u}^{\epsilon} = \mathbf{0} \quad (\text{in } \Omega_{\epsilon}), \\ \tilde{t}^{\epsilon} = \mathbf{0} \quad (\text{on } S_{\text{N}}), \quad \tilde{u}^{\epsilon} = \mathbf{0} \quad (\text{on } S_{\text{D}}), \quad \tilde{t}^{\epsilon} = -(\mathbf{C}:\nabla u) \cdot \mathbf{n} \quad (\text{on } \Gamma_{\epsilon}) \quad [6] $$

where $\sigma = \mathbf{C}:\nabla u$ is the stress tensor associated with the free field [5], $\Gamma_{\epsilon}$ is the boundary of $B_{\epsilon}$, and $\mathbf{n}$ is the normal on $S \cup \Gamma_{\epsilon}$ outward to $\Omega_{\epsilon}$. For infinitesimal $\epsilon$ the scattered field is expected to vanish, i.e. $\lim_{\epsilon \to 0} |\tilde{u}^{\epsilon}(x)| = 0$ ($x \in \Omega_{\epsilon}$), whereas the free-field, by its definition [5], does not depend on $\epsilon$. One may expand $\mathcal{J}(\Omega_{\epsilon})$ with respect to $\tilde{u}^{\epsilon}$ as

$$ \begin{aligned} \mathcal{J}(\Omega_{\epsilon}) &= \int_{S^{\text{obs}}} \varphi(\mathbf{u}^{\epsilon}(\xi), \xi) \, d\Gamma_{\xi} \\ &= \int_{S^{\text{obs}}} \left[ \varphi(\mathbf{u}(\xi), \xi) + \text{Re}\left(\frac{\partial\varphi}{\partial\mathbf{u}}(\mathbf{u}(\xi), \xi) \cdot \tilde{u}^{\epsilon}(\xi)\right) + o(|\tilde{u}^{\epsilon}(\xi)|) \right] \, d\Gamma_{\xi} \\ &= \mathcal{J}(\Omega) + \int_{S^{\text{obs}}} \text{Re}\left(\frac{\partial\varphi}{\partial\mathbf{u}}(\mathbf{u}(\xi), \xi) \cdot \tilde{u}^{\epsilon}(\xi)\right) \, d\Gamma_{\xi} + o(\|\tilde{u}^{\epsilon}\|) \end{aligned} \quad [7] $$

where

$$ \frac{\partial \varphi}{\partial \mathbf{w}} \equiv \frac{\partial \varphi}{\partial \mathbf{w}_R} - i \frac{\partial \varphi}{\partial \mathbf{w}_I} \quad (\mathbf{w}_R = \text{Re}(\mathbf{w}), \quad \mathbf{w}_I = \text{Im}(\mathbf{w})) \quad [8] $$

By means of [4] and [7], the topological derivative of $\mathcal{J}(\Omega)$ can be recast as:

$$ \mathcal{T}(\mathbf{x}^{\circ}) = \lim_{\epsilon \to 0} \frac{1}{f(\epsilon) |\mathcal{B}|} \int_{S^{\text{obs}}} \text{Re} \left( \frac{\partial \varphi}{\partial \mathbf{u}}(\mathbf{u}(\xi), \xi) \cdot \tilde{u}^{\epsilon}(\xi) \right) \, d\Gamma_{\epsilon}. \quad [9] $$

One is then left with the task of evaluating the leading contribution of the integral in [9] for $\epsilon \ll 1$.

### 4. Adjoint field approach

Define the adjoint field $\tilde{u}^{\epsilon}(\xi)$ as the solution to the elastodynamic problem

$$ \begin{aligned} \mathcal{L}(\omega)\tilde{u}^{\epsilon} &= \mathbf{0} \quad (\text{in } \Omega), \\ \tilde{t}^{\epsilon} &= \frac{\partial \varphi}{\partial \mathbf{u}}(\mathbf{u}, \cdot) \quad (\text{on } S^{\text{obs}}), \quad \tilde{t}^{\epsilon} = \mathbf{0} \quad (\text{on } S_{\text{N}} \setminus S^{\text{obs}}), \quad \tilde{u}^{\epsilon} = \mathbf{0} \quad (\text{on } S_{\text{D}}). \end{aligned} \quad [10] $$

where $\mathbf{u}$ is the free-field defined by [5], the prescribed traction is defined in terms of the cost function density $\varphi$, and the convention [8] is employed. Application of the elastodynamic reciprocity theorem to $\tilde{u}^{\epsilon}$ and $\tilde{u}^{\epsilon}$ over $\Omega_{\epsilon}$ yields

$$ \int_{\Gamma_{\epsilon}} (\tilde{u}^{\epsilon} \cdot \tilde{t}^{\epsilon} - \tilde{u}^{\epsilon} \cdot \tilde{t}^{\epsilon}) \, d\Gamma_{\epsilon} = \int_{S_{\text{obs}}} \tilde{u}^{\epsilon} \cdot \frac{\partial \varphi}{\partial \mathbf{u}}(\mathbf{u}, \cdot) \, d\Gamma_{\epsilon}, \quad [11] $$

which makes use of the boundary conditions in [6] and [10]. In practical terms, the assumed distribution of $\tilde{p}^{\epsilon}$ over $S^{\text{obs}}$ can be interpreted as being proportional to a
measure of the misfit between the measured displacement field $\mathbf{u}^{\text{obs}}$ and the (reference) free-field $\mathbf{u}$. In particular, in the least-squares case, $\partial\varphi/\partial\mathbf{u} = \mathbf{W}^{\text{T}}(\mathbf{u} - \mathbf{u}^{\text{obs}})$ so that $\dot{\mathbf{p}}$ depends linearly on the observed scattered field $\mathbf{u} - \mathbf{u}^{\text{obs}}$. On the basis of [6], [7] and [11], one finds that

$$ 
\begin{aligned}
\mathcal{J}(\Omega_{\varepsilon}) &= \mathcal{J}(\Omega) + \operatorname{Re} \left[ \int_{\Gamma_{\varepsilon}} \dot{\mathbf{u}}^{*} \cdot \tilde{\mathbf{t}}^{\varepsilon} \, \mathrm{d}\Gamma - \int_{\Gamma_{\varepsilon}} \tilde{\mathbf{u}}^{\varepsilon} \cdot \dot{\tilde{\mathbf{t}}} \, \mathrm{d}\Gamma \right] + o(\|\tilde{\mathbf{u}}^{\varepsilon}\|) \\
&= \mathcal{J}(\Omega) - \operatorname{Re} \left[ \int_{\Gamma_{\varepsilon}} \dot{\mathbf{u}}^{*} \cdot \mathbf{t} \, \mathrm{d}\Gamma + \int_{\Gamma_{\varepsilon}} \tilde{\mathbf{u}}^{\varepsilon} \cdot \dot{\tilde{\mathbf{t}}} \, \mathrm{d}\Gamma \right] + o(\|\tilde{\mathbf{u}}^{\varepsilon}\|) \quad \text{as } \varepsilon \to 0. \quad [12]
\end{aligned}
$$

For a given primary excitation $\mathbf{p}$, the first integral over $\Gamma_{\varepsilon}$ in [12] is known since the adjoint field $\dot{\mathbf{u}}^*$ is determined solely in terms of $\mathbf{u}$, $\mathbf{u}^{\text{obs}}$ and the cost function density $\varphi$ (see [10]). In fact, by virtue of [5], one has

$$ 
\begin{aligned}
\int_{\Gamma_{\varepsilon}} \dot{\mathbf{u}}^{*} \cdot \mathbf{t} \, \mathrm{d}\Gamma &= - \int_{B_{\varepsilon}} \nabla \cdot (\dot{\mathbf{u}}^{*} \cdot \mathbf{C} : \nabla \mathbf{u}) \, \mathrm{d}\Omega = \int_{B_{\varepsilon}} (\rho \omega^2 \dot{\mathbf{u}}^{*} \cdot \mathbf{u} - \nabla \dot{\mathbf{u}}^{*} : \mathbf{C} : \nabla \mathbf{u}) \, \mathrm{d}\Omega \\
&= \varepsilon^3 |\mathcal{B}| (\rho \omega^2 \dot{\mathbf{u}}^{*} \cdot \mathbf{u} - \nabla \dot{\mathbf{u}}^{*} : \mathbf{C} : \nabla \mathbf{u})(\mathbf{x}^{\circ}) + o(\varepsilon^3) \quad \text{as } \varepsilon \to 0 \quad [13]
\end{aligned}
$$

where the minus sign in front of the first integral over $B_{\varepsilon}$ appears because $\mathbf{t}$ is defined in terms of the unit normal $\mathbf{n}$ pointing to the interior of $B_{\varepsilon}$.

To complete the expansion of $\mathcal{J}(\Omega_{\varepsilon})$ in [12], one is left with evaluating the leading asymptotic behavior of the last integral for vanishing $\varepsilon$. To accomplish this task, the asymptotic behavior of $\tilde{\mathbf{u}}^{\varepsilon}$ on $\Gamma_{\varepsilon}$ must be investigated. This study [GUZ04] is based on an asymptotic analysis of the governing BIE for $\tilde{\mathbf{u}}^{\varepsilon}$ and yields the following result:

$$ \tilde{\mathbf{u}}^{\varepsilon}(\boldsymbol{\xi}) = \varepsilon \sigma_{k\ell}(\mathbf{x}^{\circ}) \mathcal{U}^{k\ell}(\bar{\boldsymbol{\xi}}) + o(\varepsilon) \quad (\boldsymbol{\xi} \in \Gamma_{\varepsilon}) \quad [14] $$

$$ \text{with} \quad \left\{ \begin{array}{l} \nabla_{\bar{\xi}} \cdot (\mathbf{C} : \nabla_{\bar{\xi}} \mathcal{U}^{k\ell}) = \mathbf{0} \quad (\bar{\xi} \in \mathbb{R}^3 \setminus \bar{\mathcal{B}}), \\ (\mathbf{C} : \nabla_{\bar{\xi}} \mathcal{U}^{k\ell}) \cdot \mathbf{n} = -(n_k e_\ell + n_\ell e_k)/2 \quad (\bar{\xi} \in \mathcal{S}), \end{array} \right. $$

where $\bar{\boldsymbol{\xi}} = \boldsymbol{\xi}/\varepsilon$ is a normalized position vector. The six canonical displacements $\mathcal{U}^{k\ell}(\bar{\boldsymbol{\xi}})$ are seen to solve normalized exterior elastostatic problems that do not depend on $\mathbf{x}^{\circ}$ and $\varepsilon$. Then, by employing [14], one obtains

$$ \int_{\Gamma_{\varepsilon}} \tilde{\mathbf{u}}^{\varepsilon} \dot{\mathbf{t}}^{*} \, \mathrm{d}\Gamma = \varepsilon^3 \dot{\sigma}_{ij}(\mathbf{x}^{\circ}) \sigma_{k\ell}(\mathbf{x}^{\circ}) \int_{\mathcal{S}} \mathcal{U}_{i}^{k\ell}(\bar{\boldsymbol{\xi}}) n_j(\bar{\boldsymbol{\xi}}) \, \mathrm{d}\vartheta_{\bar{\xi}} + o(\varepsilon^3) \quad \text{as } \varepsilon \to 0. \quad [15] $$

where $\mathrm{d}\vartheta_{\bar{\xi}}$ denotes the area differential element on the unit sphere. On substituting [15] and [13] into [12], the leading asymptotic behavior of $\mathcal{J}(\Omega_{\varepsilon})$ for vanishing $\varepsilon$ is found to take the form

$$ \mathcal{J}(\Omega_{\varepsilon}) = \mathcal{J}(\Omega) + \varepsilon^3 |\mathcal{B}| \operatorname{Re} \left[ (\dot{\boldsymbol{\sigma}} : \mathbf{A} : \boldsymbol{\sigma} - \rho \omega^2 \dot{\mathbf{u}}^{*} \cdot \mathbf{u})(\mathbf{x}^{\circ}) \right]. \quad [16] $$

Note in particular, with reference to [4], that the analysis shows that $f(\varepsilon) = \varepsilon^3$, which is also known to hold for traction-free cavities in the 3-D elastostatic case [GAR01, SOK99]. The constant fourth-order tensor $\mathbf{A}$ is defined by

$$ \mathcal{A}_{ijk\ell} = \frac{1}{2\mu} \left\{ I_{ijk\ell} - \frac{\nu}{1+\nu} \delta_{ij} \delta_{k\ell} \right\} - \frac{1}{|\mathcal{B}|} \int_{\mathcal{S}} \mathcal{U}_{i}^{k\ell}(\bar{\boldsymbol{\xi}}) n_j(\bar{\boldsymbol{\xi}}) \, \mathrm{d}\vartheta_{\bar{\xi}} \quad [17] $$
Finally, by virtue of [4] and [16], the sought adjoint-field formula for the topological derivative is

$$ \mathcal{T}(\mathbf{x}^o) = \mathrm{Re} \left[ (\hat{\boldsymbol{\sigma}}^* : \mathcal{A} : \boldsymbol{\sigma} - \rho \omega^2 \hat{\mathbf{u}}^* \cdot \mathbf{u})(\mathbf{x}^o) \right] \quad [18] $$

For an arbitrarily shaped infinitesimal cavity, the six canonical problems in [14] 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 $\mathbf{x}^o$. For the particular case of a spherical infinitesimal cavity, the problems in [14] have an analytical solution, from which the following closed-form expression of $\mathcal{A}$ is obtained:

$$ \mathcal{A}_{ijk\ell} = \frac{3(1-\nu)}{4\mu(7-5\nu)} \left[ 5(\delta_{ik}\delta_{j\ell} + \delta_{jk}\delta_{i\ell}) - \frac{1+5\nu}{1+\nu}\delta_{ij}\delta_{k\ell} \right] \quad [19] $$

Expression [18] of the topological derivative is implicit in that it relies upon solutions of boundary-value problems on the cavity-free reference body $\Omega$. An explicit expression for $\mathcal{T}(\mathbf{x}^o)$ is obtained in terms of the elastodynamic Green's tensor $\hat{\mathbf{u}}^k(\xi, \mathbf{x}, \omega)$, i.e. the elastodynamic displacement field generated by a unit point force applied at $\mathbf{x} \in \Omega$ along the $k$-th direction and satisfying the boundary conditions

$$ \hat{u}_i^k(\xi, \mathbf{x}, \omega) = 0 \quad (\xi \in S_D, \xi \neq \mathbf{x}), \qquad \hat{t}_i^k(\xi, \mathbf{x}, \omega; \mathbf{n}) = 0 \quad (\xi \in S_N, \xi \neq \mathbf{x}) $$

In this case, the free and adjoint displacement fields are given by

$$ u_k(\mathbf{x}) = \int_{S_N} p_i(\xi) \hat{u}_i^k(\xi, \mathbf{x}, \omega) \mathrm{d}\Gamma_{\xi} , \qquad \hat{u}_k(\mathbf{x}) = \int_{S_N} \hat{p}_i(\xi) \hat{u}_i^k(\xi, \mathbf{x}, \omega) \mathrm{d}\Gamma_{\xi} \quad [20] $$

and the expression [18] becomes explicit as well. Such Green's tensors are explicitly known for simple unbounded media, but not in general for finite domains.

# 5. Numerical examples

## 5.1. Pressurized annular sphere

To validate the foregoing developments, the elastodynamic Neumann problem for a spherical shell (outer radius $R$, inner radius $a^{\mathrm{true}} < R$), shown in Fig. 1, is considered. The shell, $\Omega^{\mathrm{true}}$, is subjected to a uniform time-harmonic pressure $p$ acting over its external surface $S$. For this problem, the closed-form expression for topological derivative, $\mathcal{T}(r^o)$, where $0 < r^o < R$ denotes the radial coordinate of a sampling point inside the void-free, i.e. reference body $\Omega$ (see Fig. 1), can be obtained for the least squares cost function $\mathcal{J}$ with $S^{\mathrm{obs}} = S$ and $\mathbf{W} = \mathbf{I}_2$. To provide a basis for comparison, numerical values of the topological derivative $\mathcal{T}(r^o)$ are computed by means of a three-dimensional boundary element solution [PAK99] applied to the adjoint-field formula [18] with $p = \mu$, $\nu = 0.3$, $R/a^{\mathrm{true}} = 3$ and $\omega R\sqrt{\rho/\mu} = 3$. In the numerical model, the inner and outer surfaces of the annular sphere were discretized using 486 and 150 eight-noded boundary elements, respectively, which provides at least 12 element lengths per shear wavelength for the excitation frequency chosen. The adjoint
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Spherical shell: geometry and notation (top); topological derivative $T(r^o)$ along a radial line (bottom)</figcaption></figure>

field technique for computing $\mathcal{T}(r^o)$ yields results that are in good agreement with the closed-form solution (see Fig. 1), the RMS error for the featured set of grid points being $1.5 \times 10^{-3}$. It should be mentioned, however, that the accuracy of numerical estimates of $\mathcal{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.

## 5.2. Cavity embedded in a half-space

The configuration is as depicted in Fig. 2, with the 'true' spherical cavity, of diameter $D=0.4d$, centered at $(d, 0, 3d)$. In succession, the cavity is illuminated by 16 axial point sources acting on the surface of a semi-infinite solid. For each source location, the surface motion $u^{\text{obs}}$ is monitored at 25 sensors distributed over the square testing grid; here, this data is simulated using the BEM formulation of the forward problem. Four excitation frequencies $\bar{\omega} \equiv \omega d\sqrt{\rho/\mu} = 1, 2, 4, 8$ have been considered.

For this testing configuration, the values of $\mathcal{T}(x^o)$ are computed over the horizontal surface $S_h = \{\xi \in \Omega | -5d < \xi_1 < 5d, -3d < \xi_2 < 3d, \xi_3 = 3d\}$ passing through the centroid of the 'true' cavity and plotted in Fig. 3 for the above-defined set of frequencies. The computational grid is chosen so that the sampling points $x^o$ are spaced by $0.25d$ in both $\xi_1$ and $\xi_2$ directions. In the display, the red tones indicate negative values of $\mathcal{T}$ and thus possible cavity location; for comparison, the true cavity is outlined in white in each of the diagrams. The results clearly demonstrate the usefulness of the topological derivative as a computationally efficient tool for exposing the approximate cavity location, with "higher" frequencies ($\bar{\omega} = 2, 4$) providing in general better resolution. From the diagram for $\bar{\omega} = 8$ where $\lambda_s/D \approx 1$, however, it is also
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Cavity embedded in a half-space: definition and notation</figcaption></figure>

evident that the infinitesimal-cavity assumption embedded in [4] performs best when used in conjunction with wave lengths exceeding the cavity diameter.

For completeness, the variation of $\mathcal{T}(\mathbf{x}^o)$ across the vertical planar region $S_v = \{\xi \in \Omega | -5d < \xi_1 < 5d, \xi_2 = 0, 0.25d < \xi_3 < 6d\}$ is given in Fig. 4. Similar to the earlier diagram, the sampling points $\mathbf{x}^o$ are spaced by $0.25d$ in the $\xi_1$ and $\xi_3$ directions. The fact that both source and receiver points are limited to a single planar surface clearly leads to a diminished resolution compared to the previous result. The contour plots for $\bar{\omega} = 2$ and $4$ exhibit greater accuracy than that for $\bar{\omega} = 1$, but are also plagued with local minima that are absent in the former diagram. The non-informative distribution

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Cavity embedded in a half-space: distribution of $(\mu d)^{-1} \mathcal{T}(\mathbf{x}^o)$ in the $\xi_3 = 3d$ (horizontal) plane</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4. Cavity embedded in a half-space: distribution of $(\mu d)^{-1} \mathcal{T}(x^o)$ in the $\xi_2 = 0$ (vertical) plane</figcaption></figure>

of $\mathcal{T}$ for $\bar{\omega} = 8$ indicates that the use of topological derivative in elastic-wave imaging is most effective at ‘low’ excitation frequencies, i.e. those inside the resonance region.

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

$$ \hat{T}(x^o) = \mathcal{T}(x^o)H(-C - \mathcal{T}(x^o)) \quad [21] $$

where $C < 0$ denotes a suitable threshold value and $H(\cdot)$ is the Heaviside step function. With such definition, it is also possible to combine the individual advantages of different probing wavelengths by employing the product of [21] at several frequencies. As an illustration of the latter approach, Fig. 5 plots the distribution of the product of $\hat{T}|_{\bar{\omega}=1}$ and $\hat{T}|_{\bar{\omega}=2}$ in the vertical plane, with $C$ set to approximately 40% of the global minima of the respective distributions in Fig. 4. Despite the limited accuracy and multiple minima characterizing respectively the individual solutions for $\bar{\omega} = 1$ and $\bar{\omega} = 2$, the combined result stemming from [21] points to a single cavity with its centre and size closely approximating the true void configuration.

### 5.3. Cavity embedded in a cube

To illustrate the utility of topological derivative as a preliminary tool for elastic-wave sounding of finite bodies, this example deals with the delineation of a spherical
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Distribution of $(\mu d)^{-2} \hat{T}|_{\bar{\omega}=1} \times \hat{T}|_{\bar{\omega}=2}$ in the $\xi_2 = 0$ plane</figcaption></figure>

cavity of radius $a^{\text{true}} = 0.5d$, hidden inside a solid cube of size $6d \times 6d \times 6d$. Similar to the previous example, the elastic solid is characterized by the shear modulus $\mu$, mass density $\rho$, and Poisson's ratio $\nu = 0.3$. With reference to a Cartesian coordinate system aligned with the box edges (see Fig. 6), 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(\xi)=0, \xi \in S_D)$ over the bottom patch $S_D = \{\xi \in S|-5.4d<\xi_1<5.4d, -5.4d<\xi_2<5.4d, \xi_3 = -3d\}$. Prior to the application of time-harmonic excitation used to illuminate the cavity, the rest of the external surface, $S_N = S \setminus S_D$, is traction free. As indicated in Fig. 6, 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 $\xi_1 = -3d$, the experiment parameters are chosen so that i) $\Gamma^{\text{obs}} = \{\xi \in S_N, |\xi_1 = -3d, -3d < \xi_2 < 3d, -3d < \xi_3 < 3d\}$, and ii) the cavity is illuminated in sequence by five localized (pyramid-shaped) distributions $p^q(\xi)$ ($\xi \in \Gamma^{\text{obs}}, q = 1, 2, \dots, 5$) of surface tractions, each of resultant $P = \mu d^2$, centered at $(-3d, 0, 0)$ and $(-3d, \pm 1.8d, \pm 1.8d)$, respec-

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Cavity embedded in a cube: BEM mesh</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Distribution of $T(x^o)/d$ in coordinate planes: a) $\xi_3 = d$, b), $\xi_2 = 1.5d$, and c) $\xi_1 = d$ containing the center of the true cavity (top: $\bar{\omega} = 1$, bottom: $\bar{\omega} = 2$)</figcaption></figure>

tively, applied over square patches of four boundary elements. Virtual experiments on the remaining vertical faces, i.e. $\xi_1 = 3d$, $\xi_2 = -3d$ and $\xi_2 = 3d$, are performed in an analogous fashion with the applied normal tractions (pressure) acting respectively in the negative $\xi_1$-, positive $\xi_2$-, and negative $\xi_2$-direction, respectively. The cost function $\mathcal{T}$ is of least-squares format with $W=I_2$, and its topological derivative $\mathcal{T}$ is calculated by summing the contribution from all $4 \times 5 = 20$ experiments. Again, all elastodynamic calculations are performed using the boundary element analysis in [PAK99].

Figure 7 illustrates the distribution of topological derivative corresponding to excitation frequencies $\bar{\omega} = 1$ and $\bar{\omega} = 2$. 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* [COL92] where the wave lengths are larger than the size of the scatterer: $(\lambda_s/d)|_{\bar{\omega}=1}=2\pi$ and $(\lambda_s/d)|_{\bar{\omega}=2}=\pi$, where $d$ is the cavity diameter and $\lambda_s$ denotes the shear wave length. 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 [9], is better conformed with by finite cavities that are ‘small’ relative to the probing wavelength.

In addition, Fig. 8 plots the distribution of the product $\hat{T}(x^o)|_{\bar{\omega}=1} \times \hat{T}(x^o)|_{\bar{\omega}=2}$ in three coordinate planes containing the center of the true cavity, where $\hat{T}(x^o)$ is defined by [21] with $C = 0.5 \inf_{x^o} T(x^o)$. Again, this distribution clearly points to a single cavity with its center and size closely approximating the true void configuration.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distribution of $(1/d^2)\hat{T}|_{\bar{\omega}=1} \times \hat{T}|_{\bar{\omega}=2}$ in coordinate planes: a) $\xi_3 = d$, b), $\xi_2 = 1.5d$, and c) $\xi_1 = d$ containing the center of the true cavity</figcaption></figure>

**6. Summary**

In this study, the concept of topological derivative that has its origins in elasto-
statics and shape optimization is extended to 3D elastic-wave imaging of solids. The
proposed derivation 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 de-
scribed by linear field equations such as acoustics, heat transfer, and electromagnetics,
and iii) other types of infinitely small objects, such as inhomogeneities or cracks. Nu-
merical results suggest that the topological derivative is a useful tool for outlining
hidden objects from surface measurements at a modest computational cost.

Acknowledgements

This research is currently supported by a CNRS-NSF cooperation grant. In ad-
dition, special thanks are extended to École Polytechnique for providing a Visiting
Scientist appointment to B. Guzina, and to MTS Systems Corporation for funding
visits by M. Bonnet at the University of Minnesota as a MTS Visiting Professor.

**7. Bibliography**

[COL92] COLTON D., KRESS R., *Inverse acoustic and electromagnetic scattering theory*, Springer-Verlag, 1992.

[GAR01] GARREAU S., GUILLAUME P., MASMOUDI M., "The topological asymptotic for PDE systems: the elasticity case.", *SIAM J. Contr. Opt.*, vol. 39, 2001, p. 1756-1778.

[GUZ03] GUZINA B. B., NINTCHEU FATA S., BONNET M., "On the stress-wave imaging of cavities in a semi-infinite solid", *Int. J. Solids Struct.*, vol. 40, 2003, p. 1505-1523.

[GUZ04] GUZINA B. B., BONNET M., "Topological derivative for the inverse scattering of elastic waves", *Quart. J. Mech. Appl. Math.*, vol. 57, 2004, p. 161-179.

[NIN03] NINTCHEU FATA S., GUZINA B. B., BONNET M., "A Computational Basis for Elastodynamic Cavity Identification in a Semi-Infinite Solid", *Comp. Mech.*, vol. 32, 2003, p. 370-380.

[PAK99] PAK R. Y. S., GUZINA B. B., "Seismic soil-structure interaction analysis by direct boundary element methods.", *Int. J. Solids Struct.*, vol. 36, 1999, p. 4743-4766.

[SOK99] SOKOLOWSKI J., ZOCHOWSKI A., "On the topological derivative in shape optimization.", *SIAM J. Control Optim.*, vol. 37, 1999, p. 1251-1272.