# Shape identification using acoustic measurements: a numerical investigation using BIE and shape differentiation

Marc Bonnet

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

## To cite this version:

Marc Bonnet. Shape identification using acoustic measurements: a numerical investigation using BIE and shape differentiation. IUTAM Symposium: Inverse Problems in Engineering Mechanics, 1993, Tokyo, Japan. pp.191-200, 10.1007/978-3-642-52439-4_19. hal-00122006

HAL Id: hal-00122006

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

Submitted on 20 Oct 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Shape Identification Using Acoustic Measurements:
A Numerical Investigation Using BIE and Shape
Differentiation

Marc BONNET

Laboratoire de Mécanique des Solides (CNRS, Ecole Polytechnique, Mines, Ponts et Chaussées),
Ecole Polytechnique, 91128 Palaiseau Cedex, FRANCE.

## 1 Preliminaries.

The purpose of this paper is to study numerically the problem of identifying a rigid 3D finite body Ω imbedded in an infinite acoustical medium Ωₑ = ℝ³ − Ω (wave velocity c). The cavity is subjected to a known harmonic incident pressure pᵢ(x) exp(-iωt), which satisfies Helmholtz' equation (Δ + k²)pᵢ = 0 inside Ωₑ (with k = ω/c. The scattered pressure field pΓ(x) induced by the presence of the obstacle (i.e. the solution of the direct problem) satisfies the so-called 'state equation':

$$ egin{cases} (\Delta + k^2)p = 0 & \text{in } \Omega_e \\ p_n + p_n' = 0 & \text{on } \Gamma = \partial\Omega_e = \partial\Omega (\cdot, n = \frac{\partial}{\partial n}(\cdot)) \\ \text{(radiation condition)} \end{cases} \quad (1) $$

where the unit normal n is directed outside Ωₑ, i.e. is interior to Γ.

The shape of Ωₑ, ie the surface Γ, is unknown, and p defined by (1) depends on Γ: p = pΓ. Extra data is necessary if one is to solve the inverse problem, e.g.

$$ p(x) = \hat{p}(x) \quad \text{on } C \quad (2) $$

C being a surface, or curve, exterior to Γ, on which the pressure field is measured. Thus Γ is classically sought as

$$ \min_{\Gamma} J(\Gamma), \quad J(\Gamma) \equiv J(p_{\Gamma}, \Gamma) \quad (3) $$

where:

$$ J(p, \Gamma) = \frac{1}{2} \int_C |p(y) - \hat{p}(y)|^2 dC_y \quad (4) $$

The nonlinear minimization problem (3) is best solved (in terms of both computational efficiency and accuracy) using gradient methods, such as BFGS of conjugate gradient [11]. These algorithms need repeated computations of the derivative of J(Γ) with respect to Γ (or, in practice, the design parameters which define the current location of Γ). These derivatives may be computed using finite-difference methods. However, this is computationally expensive (because the evaluation of each partial derivative needs a complete solution of (1) on a perturbed geometry Γ + δG) and may be poor in terms of accuracy.

Hence an analytical differentiation of (4) is investigated instead. Its derivation relies on the *shape differentiation approach*, and two approaches are presented here: the *direct differentiation approach (DDA)* (section 3) and the *adjoint problem approach (section 4).*

The basic nature of the inverse problem under consideration (to search an unknown surface) suggests the use of boundary integral equations (BIE) and boundary elements for the numerical
modelling of the direct problem (1). In the present problem, this indication is made even stronger by the infinite character of the domain $\Omega_e$.

Any solution $p(\mathbf{x})$ of Helmholtz' equation satisfies a well-known BIE derived from third Green's formula. In the present paper, we use its regularized version [1], [2], which reads:

$$ \gamma p(\mathbf{x}) + \int_{\Gamma} \{p(\mathbf{y})[G_{,n}(\mathbf{x}, \mathbf{y}) - G_{,n}^{0}(\mathbf{x}, \mathbf{y})] + [p(\mathbf{y}) - p(\mathbf{x})]G_{,n}^{0}(\mathbf{x}, \mathbf{y}) - q(\mathbf{y})G(\mathbf{x}, \mathbf{y})\} dS_y = 0 \quad (5) $$

In (5), $G(\mathbf{x}, \mathbf{y}) = e^{ikr}/(4\pi r)$ and $G^0(\mathbf{x}, \mathbf{y}) = 1/(4\pi r)$ are the dynamic and static Green functions, while $r = \|\mathbf{x} - \mathbf{y}\|, (\cdot),_n \equiv n_s(\mathbf{y})(\cdot),_s \equiv n_s(\mathbf{y})\partial(\cdot)/\partial y_s, q \equiv \partial p/\partial n$. The coefficient $\gamma$ in (5) depends only on the boundedness of the domain under consideration: $\gamma = 1$ (infinite medium), 1/2 (half-space), 0 (bounded domain). For a given domain, eqn. (5) holds for interior and boundary points $\mathbf{x}$ using the same $\gamma$, which is therefore not to be mistaken with the conventional free-term coefficient. In the sequel, as infinite domains are considered, eqn. (5) will be used with $\gamma = 1$. The integrand in eqn. (5) is weakly singular, provided $p(\mathbf{x}) \in C^{0,\alpha}(\Gamma)$, thanks to the regularizing effect of the term $[p(\mathbf{y}) - p(\mathbf{x})]$, and contains no Cauchy principal value (CPV) integral. This point will prove crucial for the derivation of the rate BIE in section 3 below, as well as for all numerical computations. For $p(\mathbf{x})$ solution of (1), one puts $q(\mathbf{y}) = -p_{,n}^I(\mathbf{y})$ in (5).

## 2 The shape differentiation approach.

The shape differentiation approach deals with derivatives of functionals with respect to variable domains or boundaries (e.g. $J(\Gamma)$) involving fields which themselves depend on the geometry, notably the solutions of boundary-value problems like $p_{\Gamma}(\mathbf{y})$.

Let $\Gamma$ denote the current location of the unknown boundary during the minimization process, and consider a further (small) evolution of the surface $\Gamma$ defined by means of a normal 'velocity' field $\theta(\mathbf{y})$:

$$ \mathbf{y} \in \Gamma \rightarrow \mathbf{y} + \theta(\mathbf{y})\mathbf{n}(\mathbf{y})\tau \quad \text{i.e. } \Gamma(\tau) = \Gamma + \theta\mathbf{n}\tau \quad (\tau \ge 0, \tau \text{'small'}) \quad (6) $$

while the measurement area $C$ is kept fixed ($\theta(\mathbf{y}) = 0 \quad \forall \mathbf{y} \in C$). Definition (6) is considered only for small values of $\tau$, consistently with the fact that we are only interested here by derivatives for $\tau = 0$, i.e. for the current $\Gamma$ (in the sequel, derivatives with respect to $\tau$ are always taken for $\tau = 0$).

The SDA defines several kinds of derivatives with respect to $\tau$ for fields $u(\mathbf{y}, \tau)$. In a BIE approach of shape differentiation, it seems natural to use 'material' derivatives, i.e. to 'follow' the field $u$ while the field point $\mathbf{y}$ moves according to (6), in order to use only the information available on the boundary as it 'moves'. When, as here, the geometrical transformation is described by means of a normal velocity field $\theta$, the 'transformation derivative' $\dot{u}^*$ [9] of $u = u(\mathbf{y}, t)$ is introduced:

$$ \dot{u}^*(\mathbf{y}, 0) = u_{,\tau}(\mathbf{y} + \theta n \tau, \tau) \Big|_{\tau=0} = u_{,\tau}(\mathbf{y}, 0) + \theta(\mathbf{y}) u_{,n}(\mathbf{y}, 0) \quad (7) $$

Various formulas are given in the literature (see e.g. [9], [7]) for the derivative of integrals with respect to variable volumes $\Omega$ or surfaces $\Gamma$, among which:

$$ \frac{d}{d\tau} \int_{\Omega} a(\mathbf{y}, \tau) dV_{\mathbf{y}} = \int_{\Omega} a_{,\tau}(\mathbf{y}, \tau) dV_{\mathbf{y}} + \int_{\partial\Omega} a(\mathbf{y}, \tau) \theta(\mathbf{y}) dS_{\mathbf{y}} \quad (8) $$

$$ \frac{d}{d\tau} \int_{\Gamma} a(\mathbf{y}, \tau) dS_{\mathbf{y}} = \int_{\Gamma} \{a_{,\tau}(\mathbf{y}, \tau) + (a_{,n}(\mathbf{y}, \tau) - 2K(\mathbf{y})a(\mathbf{y}, \tau))\theta(\mathbf{y})\} dS_{\mathbf{y}} \quad (9) $$

$$ = \int_{\Gamma} \{\dot{a}(\mathbf{y}, \tau) - 2K(\mathbf{y})a(\mathbf{y}, \tau)\theta(\mathbf{y})\} dS_{\mathbf{y}} \quad (10) $$

where $K(\mathbf{y})$ denotes the mean curvature at $\mathbf{y} \in \Gamma$. Equations (9), (10) hold only for a closed smooth surface, while in equation (8) $\theta$ refers to the unit normal $\mathbf{n}$ directed towards the exterior of $\Omega$. Generalization of above formulas to piecewise smooth surfaces is available [9] but will not be used here.
### 3 Rate BIE: the direct differentiation approach.

In this approach, (4) is differentiated with respect to $\tau$ directly, using eqn. (10). material derivative of integrals [9]. This yields:

$$ \frac{d}{d\tau} J(\Gamma) = \int_C (p\Gamma(y) - \hat{p}(y)) \dot{p}_{\Gamma}(y) dC_y \quad (11) $$

in which the fact that $C$ is fixed has been taken into account.

The use of (11) implies, in turn, the calculation of $\dot{p}_{\Gamma}$. This is achievable by deriving from the state BIE (5) a rate BIE, in which the fields $(\dot{p}_{\Gamma}, \dot{q}_{\Gamma})$ appear. The rate BIE is indeed obtained upon application of formula (10) to (5). The validity of the differentiation of (5) using (10) relies upon the weakly singular character of the acoustic BIE (4) (as opposed to conventional CPV BIEs, for which extreme care must be taken when considering their differentiation with respect to a parameter).

As a result of this differentiation process, the rate BIE reads, in the present context:

$$ \begin{aligned} & \dot{p}(\mathbf{x}) + \int_{\Gamma} \{ \dot{p}(\mathbf{y})[G_{,n}(\mathbf{x}, \mathbf{y}) - G_{,n}^{0}(\mathbf{x}, \mathbf{y})] + [\dot{p}(\mathbf{y}) - \dot{p}(\mathbf{x})]G_{,n}^{0}(\mathbf{x}, \mathbf{y}) - \dot{q}(\mathbf{y})G(\mathbf{x}, \mathbf{y}) \} dS_y = \\ & = \int_{\Gamma} q(\mathbf{y}) \{[\theta(\mathbf{y})n_s(\mathbf{y}) - \theta(\mathbf{x})n_s(\mathbf{x})]G_{,s}(\mathbf{x}, \mathbf{y}) - 2K(\mathbf{y})\theta(\mathbf{y})G(\mathbf{x}, \mathbf{y})\} dS_y \\ & \quad - \int_{\Gamma} [\theta(\mathbf{y})n_r(\mathbf{y}) - \theta(\mathbf{x})n_r(\mathbf{x})] \{D_{rs}p(\mathbf{y})G_{,s}(\mathbf{x}, \mathbf{y}) - k^2p(\mathbf{y})G(\mathbf{x}, \mathbf{y})\} dS_y \end{aligned} \quad (12) $$

where $D_{rs}$ is the tangential differential operator given by eqn. (30) of Appendix A. The derivation of (12) uses the integration by parts formula (35) and:

$$ \frac{d}{d\tau} G(\mathbf{x}, \mathbf{y}) = [\theta(\mathbf{y})n_s(\mathbf{y}) - \theta(\mathbf{x})n_s(\mathbf{x})]G_{,s}(\mathbf{x}, \mathbf{y}) $$

**Remark 1** The rate BIE as given by (12) is valid for a smooth surface, without edges or corners (but the present approach can be extended to piecewise smooth surfaces).

**Remark 2** All integrands in eqn. (12) are weakly singular, thanks to the regularizing effect of $[\dot{p}(\mathbf{y}) - \dot{p}^*(\mathbf{x})]$, $[\theta(\mathbf{y})n_s(\mathbf{y}) - \theta(\mathbf{x})n_s(\mathbf{x})]$ and $[p(\mathbf{y}) - p(\mathbf{x})]$.

**Remark 3** The validity of the derivation of (12) relies upon the validity of the differentiation of (5) using (10). In this respect, the weakly singular character of the acoustic BIE (5) is crucial. At any stage of the derivation, all integrals are at most weakly singular. In contrast, a very close attention should be exercised with respect to the handling of exclusion neighbourhoods and subsequent limit processes, were the same approach to be applied to conventional strongly singular BIEs.

As a consequence of (1), the expression of $\dot{q}_{\Gamma} = -(p_{,n}^{I})^*$ on $\Gamma$ is given (see appendix B) by:

$$ (p_{,n}^{I})^* = (2K p_{,n}^{I} - k^2 p^{I}) \theta - D_s (\theta D_s p^{I}) \quad (13) $$

and the right-hand side $\mathcal{F}(p, q)\theta$ in (12) becomes

$$ \begin{aligned} \mathcal{F}(p, q)\theta = & - \int_{\Gamma} p_{,n}^{I}(\mathbf{y})(\theta(\mathbf{y})n_s(\mathbf{y}) - \theta(\mathbf{x})n_s(\mathbf{x}))G_{,s}(\mathbf{x}, \mathbf{y}) dS_y \\ & + \int_{\Gamma} \{D_s(\theta D_s p^{I})(\mathbf{y}) - (\theta(\mathbf{y})n_r(\mathbf{y}) - \theta(\mathbf{x})n_r(\mathbf{x}))D_{rs}p(\mathbf{y})\} G_{,s}(\mathbf{x}, \mathbf{y}) dS_y \end{aligned} \quad (14) $$

The representation formula for interior values $p(\mathbf{x})$ may also be differentiated with respect to the variable domain. The result is:

$$ \begin{aligned} \dot{p}(\mathbf{x}) = & - \int_{\Gamma} \{[\theta(\mathbf{y})p_{,n}^{I}(\mathbf{y}) + \dot{p}_{\Gamma}(\mathbf{y})]G_{,n}(\mathbf{x}, \mathbf{y}) \\ & + \theta(\mathbf{y}) (k^2(p_{\Gamma} + p^{I})(\mathbf{y})G(\mathbf{x}, \mathbf{y}) - D_s(p_{\Gamma} + p^{I})D_s G(\mathbf{x}, \mathbf{y}))\} dS_y \end{aligned} \quad (15) $$
Then, in order to evaluate the derivative (11), one has to

1. Solve the state BIE (5) for the unknown $p_{\Gamma}$.

2. Solve the rate BIE (14) for the unknown $p_{\Gamma}^*$. This must be done successively with each $\theta$ in turn (or rather, in practice, with a finite number of $\theta$s associated to design variables).

3. Compute the field $p_{\Gamma}^*$ on the measurement surface $C$ using representation formula (15).

4. Insert, for each $\theta$, $p_{\Gamma}$ and $p_{\Gamma}^*$ in $\frac{dJ(\Gamma)}{dt}$ (11).

Indeed, as soon as the state BIE (5) is solved, the right-hand side $\mathcal{F}(p, q)\theta$ of (12) is completely known. The successive solution of (14) for various $\theta$s may seem at first sight to be a formidable computational task. However, as the same integral operator governs the fields $p$ in (5) and $\dot{p}$ in (12), the actual solution of (12) needs only, for each $\theta$, the building of its right-hand side $\mathcal{F}(p, q)\theta$ followed by a backsubstitution using the already factored matrix. In other words, each computation of the derivative needs the building and factorization of only one matrix.

The field $\dot{p}_{\Gamma}^*$ appears clearly to be a linear form over $\theta$. Thus, with sufficiently regular $\Gamma$, $\theta$ and $p^!$, Riesz representation theorem may be used to state the existence of a kernel $\frac{\partial}{\partial \Gamma} p(x, y)$ (the derivative of $p$ with respect to $\Gamma$, such that:

$$ \dot{p}(x) = \int_{\Gamma} \frac{\partial p}{\partial \Gamma}(x, y)\theta(y) dS_y \quad (16) $$

hence the terminology of 'shape differentiation'. As a result, $dJ/d\tau$ itself is, from (11) and the previous remarks, a linear form over $\theta$, and the shape derivative kernel $\frac{\partial}{\partial \Gamma} J(y)$ may be defined as well.

The DDA has also been applied to Galerkin type BIEs in [8], for crack identification problems governed by Laplace's equation. As Galerkin-type BIEs are weakly singular (they are regularized by integrating by parts twice), the DDA approach does not raise difficulties.

## 4 The adjoint problem approach.

This is an alternative approach for deriving an analytical expression, which is known e.g. in the field of structural shape optimization [7] (see also [6] for thermal inverse problems). The problem (3) may be viewed as a constrained optimization problem, where $J(p, \Gamma)$ is to be minimized under the constraint $p = p_{\Gamma}$ (1). The latter can be expressed in weak form as:

$$ A(p, w; \Gamma) = \int_{\Omega_e} (\nabla p \cdot \nabla w - k^2 p w) dV_y + \int_{\Gamma} w p_{,n}^I dS_y = 0 \quad \forall w \in V \quad (17) $$

where $V = \{w \in H^1(\Omega), w(y) = 0 (y \in S)\}$. A Lagrangian functional $\mathcal{L}$ is introduced:

$$ \mathcal{L}(p, w, \Gamma) = J(p, \Gamma) + A(p, w; \Gamma) \quad (18) $$

where $w$ is the Lagrange multiplier. Upon application of formulas (8), (9), taking into account identity (13) and using identity (34), the stationarity of $\mathcal{L}$ is expressed as

$$ \frac{d}{d\tau} \mathcal{L} = \frac{\partial \mathcal{L}}{\partial p} p_{,\tau} + \frac{\partial \mathcal{L}}{\partial \Gamma} \theta \quad (19) $$

$$ \frac{\partial \mathcal{L}}{\partial p} p_{,\tau} = \int_C p_{,\tau}(p - \hat{p}) dC_y + \int_{\Omega_e} (\nabla p_{,\tau} \cdot \nabla w - k^2 p_{,\tau} w) dV_y \quad (20) $$

$$ \frac{\partial \mathcal{L}}{\partial \Gamma} \theta = \int_{\Gamma} \theta \left[ \nabla_S w \cdot \nabla_S (p + p^!) - k^2 w (p + p^!) \right] dS_y \quad (21) $$
Now the choice of the Lagrange multiplier $w$ is restricted in such a way that $\frac{d}{d\tau}\mathcal{L} = 0$ for $\theta \equiv 0$, that is, we put:

$$ \frac{\partial \mathcal{L}}{\partial p} p_{,\tau} = 0, \quad \forall p_{,\tau} \in \mathcal{V} \qquad (22) $$

This leads, in view of eqn. (20), to introduce the adjoint field $w_\Gamma$ as the solution of the following adjoint problem:

$$ \left\{ \begin{array}{ll} (\Delta + k^2)w = -(p - \hat{p})\delta_C & \text{in } \Omega_e \\ w_{,n} = 0 & \text{on } \Gamma \\ \text{(radiation condition)} & \end{array} \right. \qquad (23) $$

The adjoint problem is equivalently formulated in terms of a BIE:

$$ \begin{aligned} & w(\mathbf{x}) + \int_{\Gamma} \{ w(\mathbf{y})n_s(\mathbf{y})[G_{,s}(\mathbf{x}, \mathbf{y}) - G_{,s}^{0}(\mathbf{x}, \mathbf{y})] + [w(\mathbf{y}) - w(\mathbf{x})]n_s(\mathbf{y})G_{,s}^{0}(\mathbf{x}, \mathbf{y}) \} \, dS_y \\ & = \int_{C} [p(\mathbf{y}) - \hat{p}(\mathbf{y})]G(\mathbf{x}, \mathbf{y}) \, dC_y \end{aligned} \qquad (24) $$

Finally the shape derivative of $J$ is given by:

$$ \frac{dJ}{d\tau} = \frac{d}{d\tau} \mathcal{L}(p_\Gamma, w_\Gamma, \Gamma) = \int_\Gamma \theta \left[ \nabla_S w_\Gamma \cdot \nabla_S (p_\Gamma + p^I) - k^2 w_\Gamma (p_\Gamma + p^I) \right] dS \qquad (25) $$

Then, in order to evaluate the derivative $dJ/d\tau$, one has to

1. Solve the state BIE (5) for the unknown $p_\Gamma$.

2. Solve the adjoint BIE (24) associated to the adjoint problem (23) for the unknown $w_\Gamma$.

3. Insert $p_\Gamma$ and $w_\Gamma$ in $dJ/d\tau$ (25).

**Remark 1** The adjoint problem approach is not specifically BIE-oriented. Indeed, its establishment uses weak forms of boundary-value problems and is therefore equally well suited for FEM numerical treatments.

**Remark 2** The adjoint field $w$ does not depend on $\theta$. Therefore, the adjoint problem approach needs the solution of two distinct boundary-value problems.

**Remark 3** As for the DDA, the state and adjoint problems use the same integral operator. Therefore, each computation of the derivative needs the building and factorization of only one matrix.

# 5 Numerical implementation and examples.

Numerical tests for the solution of the inverse problem using shape differentiation have been conducted, for 3D situations. The basic numerical tool is the regularized collocation BIE (5), which is implemented in our BEM research code ASTRID. Isoparametric 8-noded curved surface elements were used throughout the present study, and $\Gamma$ is made of 24 elements, which amounts to 74 nodes. The incident wave $p^I(\mathbf{y})$ is here taken as a plane wave propagating along $Ox_3$ in the positive direction.

Here the unknown surface $\Gamma$ is searched as an ellipsoid. This choice reduces the inverse problem to the search of 9 design parameters $d_1, \dots, d_9$ (coordinates $x_G, y_G, z_G$ of the center, Euler angles $\phi, \theta, \psi$ and principal axes $a, b, c$, in this order), it has been made in order to test the method on situations with moderate number of design parameters. In the next step (left for future work), mesh nodes coordinates will be taken as the design parameters. The following one-to-one mapping between $\Gamma$ and the unit sphere $S$ is used:

$$ Y = (Y_1, Y_2, Y_3) \in S \rightarrow y = (y_1, y_2, y_3) \qquad \left\{ \begin{array}{l} y_1 = x_G + r_{11}aY_1 + r_{12}bY_2 + r_{13}cY_3 \\ y_2 = y_G + r_{21}aY_1 + r_{22}bY_2 + r_{23}cY_3 \\ y_3 = y_G + r_{31}aY_1 + r_{32}bY_2 + r_{33}cY_3 \end{array} \right. \qquad (26) $$
where $r_{ij} = r_{ij}(\phi, \theta, \psi)$ are the components of the rotation which transforms the coordinate axes onto the principal axes of the ellipsoid (several $(\phi, \theta, \psi)$ triplets may define the same rotation). Hence, taking $\tau = d_i$, $\theta_i(y)$ is defined by considering an increment $\delta d_i$ added to $d_i$:

$$ \delta y = \frac{\partial y}{\partial d_i} \delta d_i = \theta_i(y) \delta d_i \qquad (27) $$

To be more precise, the sphere $S$ is meshed using 24 elements (3 per octant) and the mesh of $\Gamma$ is defined by applying the mapping (26) between the nodes of $S$ and of $\Gamma$. Hence eqn. (27) defines the nodal values of $\theta$, the latter being is interpolated using the same shape functions as the geometry.

In all computations, the measurement surface $C$ was a sphere of radius 10 metres, centered at the origin, and the wave velocity was equal to unity. Two meshes of measurement points were used (respectively 24 elements, 74 points and 96 elements, 290 points). The elements (and the shape functions associated) were used for the computation of integrals like (4) or (25) over $C$. In some cases, measurement noise has been artificially introduced by multiplication of the data values by $1+r$, $r$ being random numbers uniformly distributed in $[-\epsilon, \epsilon]$, with $\epsilon = 10^{-3}, 10^{-2}, 10^{-1}$).

Both conjugate gradient (CG) and BFGS variable metric algorithms have been applied to the minimization problem (3), using programs from [10]. All numerical computations have been done in double precision complex arithmetic, on HP-Apollo 400 type workstations.

Numerical results are presented below for four situations. Each of them is defined by the data of 'true' and initial values of $(d_1, \dots, d_9)$, according to the table 1 below. In the last example (case 4), the 'true' obstacle is not an ellipsoid but a rectangular box of sides $2a, 2b, 2c$, in order to test the performance of the method on a case where the searched-for shape differs from the true one.

As the same ellipsoid can result from many combinations of Euler angles and permutations of principal axes, it is difficult to measure the accuracy of the identification of $\Gamma$ by means of a mere comparison of the identified parameters $d_i$ with those defining the 'true' $\Gamma$, used to compute the simulated data. Instead, the relative errors $e_v, e_A, e_I$ for the volume, boundary area and geometrical inertia tensor (with respect to the fixed origin and axes have been computed (the indicator $e_I$ being very sensitive to the orientation of $\Gamma$ in space), together with the ratio $J_n/J_0$, where $J_n = J(\Gamma_n), \Gamma_n$ being the current $\Gamma$ after the n-th iteration of the minimization process. Our numerical values of $J_{final}/J_0, e_v, e_A, e_I$ obtained for the cases defined in table 1 are displayed in tables 2, 3, 4, 5.

* remark 1 Cases 1,2,3, where the 'true' cavity is also an ellipsoid, exhibit very good convergence and accuracy, especially for non-perturbed data.

* remark 2 The convergence and accuracy remain good for case 4, where the 'true' cavity is a rectangular box and exact convergence is hence impossible. The 'final' ellipsoid found by the algorithm has very similar volume and inertia tensor than the box (see table 5) and slightly different area. It even seems to be less sensitive to data noise.

* remark 3 At least in the range $\epsilon = 10^{-3}$ to $10^{-1}$], the error indicators $e_v, e_A, e_I$ vary linearly with $\epsilon$ in the results presented here, while $J_{final}/J_0$ vary quadratically. The numerical solution of the inverse problem hence behaves well with respect to measurement noise. This is probably a consequence of the strong assumption made on the unknown geometry, which is described using only 9 parameters.

* remark 4 Both CG and BFGS minimization methods have been tried, using library routines and without any optimization attempt. Neither of them seems to be significantly more accurate than the other: the final convergence and accuracy are similar. The CG method seems to exhibit slower convergence in the final stage.
# 6 Concluding comments.

A very strong assumption has been made here on the shape of $\Gamma$, and this is probably one of the reasons of the good convergence of the nonlinear optimization here. and of the fairly good behaviour with respect to data noise.

Hence, a crucial step in the numerical tests will be to relax this a priori information and to allow the search for more general surfaces, notably by using the node coordinates as the design variables. In this event, difficulties are expected:

* The general inverse problem is expected to be ill-posed (high sensitivity to data noise), and undesirable oscillations in the recovered shape of $\Gamma$ are to be expected. Hence, in the formulations outlined in sections 4 and 3 above, the functional $J(\Gamma)$ may be replaced by $J(\Gamma) + \alpha\Omega(\Gamma)$, where $\alpha > 0$ and $\Omega(\Gamma)$ is a (positive) stabilizing functional (Tikhonov regularization [12]), in order to cater for the ill-posedness by adding qualitative a priori information. For instance, one may consider $\Omega(\Gamma) = \int_{\Gamma} (D_s n_s)^2 dS$, which is expected to damp oscillations. This achieves a trade-off between accuracy and stability with respect to data noise.

* Element distortion may occur as the minimization proceeds, and must be controlled.

The assumption made on $\Gamma$ in the results presented here may indeed be viewed as a regularization, as some selection has been done in the space of parameters (i.e. of all possible $\Gamma$).

The work presented here Finally, let us mention that rate BIEs or the adjoint problem approach can be derived for elastodynamics as well; work is under progress.

# References

[1] Bonnet M. - Méthode des équations intégrales régularisées en élastodynamique tridimensionnelle. PhD thesis (Ecole Nationale des Ponts et Chaussées, Paris, France), Bulletin EDF/DER série C, n° 1/2, 1987.

[2] Bonnet M. - Regularized Boundary Integral Equations for Three-dimensional Bounded or Unbounded Elastic Bodies Containing Curved Cracks of Arbitrary Shape Under Dynamic Loading. In "Boundary Element Techniques: Applications in Engineering (ed. C.A. Brebbia & N.G. Zamani), Computational Mechanics Publications (Southampton), 1989.

[3] Bonnet M., Bui H.D. - On Some Inverse Problems for Determining Volumic Defects by Electric Current Using B.I.E. Approaches : an Overview. 6th National Japonese Conf. on Boundary Elements Methods., Tokyo, December 1989.

[4] Bonnet M. - Shape differentiation of regularized BIE: application to 3-D crack analysis by the virtual crack extension approach. In "Boundary Elements in Mechanical and Electrical Engineering" (C.A. Brebbia and A. Chaudouet, eds.), Computational Mechanics Publications, Springer-Verlag, 1990.

[5] Bonnet M. and Bui H.D. - Regularization of the Displacement and Traction BIE for 3D Elastodynamics using indirect methods, To appear in: *Advances in Boundary Element Techniques*, J.H. Kane, G. Maier, N. Tosaka and S.N. Atluri, eds., Springer-Verlag, 1992.

[6] Bonnet M. and Bui H.D. - Identification of Heat Conduction Coefficient: Application to Nondestructive Testing. *IUTAM Symposium on Inverse Problems in Engineering Mechanics* (Tokyo, 11-15 may 1992), H.D. Bui & M. Tanaka, eds., Springer-Verlag, 1992.

[7] Haug, Choi, Komkov - Design Sensitivity Analysis of Structural Systems, Academic Press, 1986.

[8] Nishimura N, Kobayashi S. - A boundary integral equation method for an inverse problem related to crack detection. Int. J. Num. Meth. Eng. 32, 1371-1387, 1991.
[9] Petryk H., Mroz Z. - Time derivatives of integrals and functionals defined on varying volume and surface domains. Arch. Mech. vol. 38, n°5-6, pp 697-724, 1986.

[10] Press W.H., Flannery B.P., Teukolsky S.A., Vetterling W.T. - *Numerical recipes: the art of scientific computing*. Cambridge press, 1986.

[11] Stoer J., Bulirsch R. - *Introduction to numerical analysis*. Springer-Verlag, 1980.

[12] Tikhonov A.N., Arsenin V.Y. - Solutions to ill-posed problems. Winston-Wiley, New York, 1977.

<figure><img src="image_3.png" /><figcaption>Figure 1: Convergence of $J_{final}/J_0$, $\epsilon_V$, $\epsilon_A$, $\epsilon_I$, for Case 2, with $\epsilon = 0$, $10^{-3}$, $10^{-2}$, $10^{-1}$.</figcaption></figure>
<table><caption>Table 1: Wavenumber values and 'true' and initial ($d_1, ..., d_9$) used in the numerical calculations.</caption><thead><tr><th></th><th></th><th>x<sub>G</sub></th><th>y<sub>G</sub></th><th>z<sub>G</sub></th><th>φ</th><th>θ</th><th>ψ</th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><td rowspan="2">Case 1<br>(ω/c = 1)</td><td>True</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>1.0</td><td>2.0</td><td>1.0</td></tr><tr><td>Initial</td><td>2.0</td><td>-1.0</td><td>1.0</td><td>1.0</td><td>2.0</td><td>3.0</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td rowspan="2">Case 2<br>(ω/c = .25)</td><td>True</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.4</td><td>0.9</td><td>0.6</td><td>1.0</td><td>2.0</td><td>1.0</td></tr><tr><td>Initial</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>5.0</td><td>5.0</td><td>5.0</td></tr><tr><td rowspan="2">Case 3<br>(ω/c = 1)</td><td>True</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.4</td><td>0.9</td><td>0.6</td><td>1.0</td><td>2.0</td><td>1.0</td></tr><tr><td>Initial</td><td>1.0</td><td>1.0</td><td>0.5</td><td>0.4</td><td>0.9</td><td>0.6</td><td>1.2</td><td>2.4</td><td>1.2</td></tr><tr><td rowspan="2">Case 4 (box)<br>(ω/c = .5)</td><td>True</td><td>1.0</td><td>0.0</td><td>-2.0</td><td>0.4</td><td>0.9</td><td>0.6</td><td>1.0</td><td>3.0</td><td>1.0</td></tr><tr><td>Initial</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>1.5</td><td>1.5</td><td>1.5</td></tr></tbody></table>

<table><caption>Table 2: Results for case 1 (convergence after about 40 CG iterations).</caption><thead><tr><th>Case 1</th><th>ε = 0.</th><th>ε = 10<sup>-3</sup></th><th>ε = 10<sup>-2</sup></th><th>ε = 10<sup>-1</sup></th></tr></thead><tbody><tr><td>J<sub>final</sub>/J<sub>0</sub></td><td>3.04 10<sup>-5</sup></td><td>3.43 10<sup>-6</sup></td><td>1.35 10<sup>-5</sup></td><td>7.88 10<sup>-4</sup></td></tr><tr><td>e<sub>V</sub></td><td>6.98 10<sup>-4</sup></td><td>1.28 10<sup>-3</sup></td><td>6.47 10<sup>-3</sup></td><td>5.80 10<sup>-2</sup></td></tr><tr><td>e<sub>A</sub></td><td>5.44 10<sup>-4</sup></td><td>9.42 10<sup>-4</sup></td><td>4.59 10<sup>-3</sup></td><td>4.00 10<sup>-2</sup></td></tr><tr><td>e<sub>I</sub></td><td>6.58 10<sup>-3</sup></td><td>6.7810<sup>-3</sup></td><td>1.33 10<sup>-2</sup></td><td>1.08 10<sup>-1</sup></td></tr></tbody></table>

<table><caption>Table 3: Results for case 2 (convergence after about 90 BFGS iterations).</caption><thead><tr><th>Case 2</th><th>ε = 0.</th><th>ε = 10<sup>-3</sup></th><th>ε = 10<sup>-2</sup></th><th>ε = 10<sup>-1</sup></th></tr></thead><tbody><tr><td>J<sub>final</sub>/J<sub>0</sub></td><td>1.57 10<sup>-12</sup></td><td>7.14 10<sup>-11</sup></td><td>6.5710<sup>-9</sup></td><td>6.7010<sup>-7</sup></td></tr><tr><td>e<sub>V</sub></td><td>1.48 10<sup>-5</sup></td><td>5.65 10<sup>-4</sup></td><td>5.3510<sup>-3</sup></td><td>5.3110<sup>-2</sup></td></tr><tr><td>e<sub>A</sub></td><td>1.34 10<sup>-4</sup></td><td>6.16 10<sup>-4</sup></td><td>3.8710<sup>-3</sup></td><td>1.8410<sup>-2</sup></td></tr><tr><td>e<sub>I</sub></td><td>3.20 10<sup>-4</sup></td><td>1.29 10<sup>-3</sup></td><td>6.5410<sup>-4</sup></td><td>7.7310<sup>-2</sup></td></tr></tbody></table>

<table><caption>Table 4: Results for case 3 (convergence after about 25 CG iterations).</caption><thead><tr><th>Case 3</th><th>ε = 0.</th><th>ε = 10<sup>-3</sup></th><th>ε = 10<sup>-2</sup></th><th>ε = 10<sup>-1</sup></th></tr></thead><tbody><tr><td>J<sub>final</sub>/J<sub>0</sub></td><td>4.78 10<sup>-9</sup></td><td>3.44 10<sup>-7</sup></td><td>3.45 10<sup>-5</sup></td><td>3.50 10<sup>-3</sup></td></tr><tr><td>e<sub>V</sub></td><td>3.48 10<sup>-6</sup></td><td>5.47 10<sup>-4</sup></td><td>5.42 10<sup>-3</sup></td><td>5.43 10<sup>-2</sup></td></tr><tr><td>e<sub>A</sub></td><td>7.61 10<sup>-6</sup></td><td>3.48 10<sup>-4</sup></td><td>3.41 10<sup>-3</sup></td><td>3.40 10<sup>-2</sup></td></tr><tr><td>e<sub>I</sub></td><td>1.75 10<sup>-5</sup></td><td>7.82 10<sup>-4</sup></td><td>7.73 10<sup>-3</sup></td><td>7.82 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>Table 5: Results for case 4 (convergence after 32 BFGS iterations).</caption><thead><tr><th>Case 4</th><th>ε = 0.</th><th>ε = 10<sup>-1</sup></th></tr></thead><tbody><tr><td>J<sub>final</sub>/J<sub>0</sub></td><td>2.79 10<sup>-4</sup></td><td>8.85 10<sup>-4</sup></td></tr><tr><td>e<sub>V</sub></td><td>1.68 10<sup>-2</sup></td><td>3.66 10<sup>-2</sup></td></tr><tr><td>e<sub>A</sub></td><td>1.51 10<sup>-1</sup></td><td>1.25 10<sup>-1</sup></td></tr><tr><td>e<sub>I</sub></td><td>5.71 10<sup>-2</sup></td><td>2.37 10<sup>-2</sup></td></tr></tbody></table>
A Tangential differential operators and integration by parts.

Let $S$ be a twice continuously differentiable closed ($C^2$) surface, of unit normal $n$ (open surfaces can be considered as well, see e.g. [5]). Consider a scalar field $u(y)$, $y \in S$, which may be undefined outside $S$ (e.g. $u = n_i(y)$, $u = \theta(y)$). In this case, the cartesian derivatives $u_i$ are generally meaningless, and one has to introduce tangential differential operators. The domain of definition of $u$ is extended in a neighbourhood $V$ of $S$ by introducing a continuation $\hat{u}$ of $u$ outside $S$ defined as: $\forall (y \in V)$, $\hat{u}(y) = u(P(y))$, where $P(y)$ is the orthogonal projection of $y$ onto $S$. Clearly the restriction of $\hat{u}$ to $S$ is equal to $u$. Moreover the normal derivative of $\hat{u}$ is equal to zero, i.e. the vector $\nabla \hat{u}$ is tangent to $S$; therefore it may be used to define the tangential gradient $\nabla_S u$ of the function $u$;

$$ \nabla_S u = \nabla_S \hat{u} = \nabla \hat{u} \qquad (28) $$

If $u$ is an arbitrary scalar function defined in $V$, one has, consistently with (28):

$$ \nabla_S u = \nabla u - n u_n = e_r \hat{D}_r u = e_r (u_r - n_r u_n) \qquad (29) $$

which defines the tangential partial derivatives $\hat{D}_r u$ (using the notation $\cdot, n = \partial/\partial n(\cdot)$). In the following, the symbol ($\cdot$) will be omitted, keeping in mind when necessary the extension. The operator

$$ D_{rs} f = (n_r f_s - n_s f_r) \qquad (30) $$

is also introduced. From (29), $D_{rs} f = n_r D_s f - n_s D_r f$: $D_{rs} f$ is a tangential differential operator.

An interesting consequence of (29) is the following identity for the Laplace operator:

$$ \Delta u = u_{nn} - 2K u_n + D_s D_s u \qquad (31) $$

$$ 2K = -D_s n_s \qquad (32) $$

The classical Stokes' identity for a vector field $U$ defined over $V$ reads:

$$ \int_S n \cdot \operatorname{rot}(U) dS = 0 \qquad (33) $$

Application of identity (33) to $U = (n \wedge e_j) f$ and $U = n \wedge (e_j \wedge n) f$ yields integration by parts identities associated to the tangential differential operators (29), (30):

$$ \int_S (-n_r K f + D_r f) dS_y = 0 \qquad (34) $$

$$ \int_S D_{rs} u \, dS_y = 0 \quad \text{for any fixed pair } r, s, \quad r, s = 1, 2, 3 \qquad (35) $$

Identity (35) is very interesting for BEM formulations: it allows integration by parts on surfaces using ordinary partial derivatives (i.e. without separation of tangential and normal derivatives), thanks to eqn. (30). It is a key tool for regularization techniques [6].

B Some auxiliary formulas for shape differentiation.

The following formulas hold (they can be found e.g. in [9]):

$$ (u_s)^* = \dot{u}_s - u_n \theta_s - \theta D_r u D_s n_r \qquad (36) $$

$$ \dot{n} = -\nabla_S \theta = -D_r \theta e_r \qquad (37) $$

A combined application of (36) and (37) leads to:

$$ (u_n)^* = \dot{u}_n - u_n \theta_n - D_s u D_s \theta \qquad (38) $$

Moreover, if $u(y, \tau) = u(y)$, then:

$$ \dot{u}(y) = \theta(y) u_n(y) \qquad (39) $$

$$ (u_n)^* = \theta u_{nn} - D_s u D_s \theta \qquad (40) $$

In the particular case $u = p^I$, eqn. (40) holds and identity (31), combined with Helmholtz' equation, gives eqn. (13).