A variational approach to the inverse imaging of composite elastic materials

TL;DR

This paper presents a variational framework for inverse imaging of composite elastic materials using surface acoustic wave data, combining mathematical modeling, numerical methods, and gradient flow analysis to recover interior composition.

Contribution

It introduces a novel variational approach with gradient-based optimization for reconstructing composite material interiors from boundary acoustic data.

Findings

Gradient flow effectively recovers interior composition

Numerical methods successfully solve the inverse problem

Framework demonstrates accurate reconstruction in simulations

Abstract

We introduce a framework for performing the inverse imaging of composite elastic materials. Our technique uses surface acoustic wave (SAW) boundary observations within a minimization problem to express the interior composition of the composite elastic materials. We have approached our target problem by developing mathematical and computational methods for investigating the numerical solution of the corresponding inverse problem. We also discuss a mathematical model for expressing the propagation of elastic waves through composite elastic bodies, and develop approximation schemes for investigating its numerical solutions. Using these methods, we define a cost functional for measuring the difference between simulated and given SAW data. Then, using a Lagrangian approach, we are able to determine the gradient of the cost functional and analyze the inverse imaging problem's solution as a gradient flow. The cost functional's gradient is composed of solutions to a state equation, as well as of solutions to related adjoint problems. We thus developed numerical methods for solving these problems and investigated the gradient flow of the cost functional. Our results show that the gradient flow is able to recover the interior composition of the composite, and we illustrate this fact using the numerical realization of our proposed framework.

Figures

Tables

$𝒇(𝒙,t)={(f_1,f_2)}^T$	: a time dependent outer force in the ${(x_1,x_2)}^T$ direction
$(l_x,l_y)=(0.5,0.98)$	: central focus of the outer force
$l_N=4$	: number of oscillations within the time interval $(0,T)$
$A={10}^{10}$	: a parameter controlling the magnitude of the outer force.
$l_{width}=1300$	: a parameter controlling the extent of the outer force.

Equations

$$\Omega ={\Omega }_1\cup {\Omega }_2.$$

$$\{\begin{array}{ll}\rho (\mathbf{x}){\mathbf{u}}_{tt}=\mathrm{div}(c(\mathbf{x})\otimes \nabla \mathbf{u})+\mathbf{f}(\mathbf{x},t)& \text{in}(0,T)\times \Omega \\ u_i=g_i& \text{on}(0,T)\times {\Gamma }_{g_i}(\text{Dirichlet B.C.})\\ (c(\mathbf{x})\otimes \nabla \mathbf{u}{)}_{ij}{n}_j=h_i& \text{on}(0,T)\times {\Gamma }_{h_i}(\text{Neumann B.C.})\\ [u_i(\mathbf{x},t)]=0& \text{on}(0,T)\times \partial {\Omega }_1\cap \partial {\Omega }_2\\ [(c(\mathbf{x})\otimes \nabla \mathbf{u}{)}_{ij}]=0& \text{on}(0,T)\times \partial {\Omega }_1\cap \partial {\Omega }_2\\ u_i(0,\mathbf{x})=u_{0i}(\mathbf{x})& \mathbf{x}\in \Omega \\ u_{i,t}(0,\mathbf{x})={\dot{u}}_{0i}(\mathbf{x})& \mathbf{x}\in \Omega .\end{array}$$

$$\partial \Omega \backslash \{\partial {\Omega }_1\cap \partial {\Omega }_2\}=\left\{\bigcup _i{\Gamma }_{h_i}\right\}\cup \left\{\bigcup _i{\Gamma }_{g_i}\right\}$$

$$\nabla \mathbf{u}=\left(\begin{array}{ccc}\frac{\partial u_1}{\partial x}& \frac{\partial u_1}{\partial y}& \frac{\partial u_1}{\partial z}\\ \frac{\partial u_2}{\partial x}& \frac{\partial u_2}{\partial y}& \frac{\partial u_2}{\partial z}\\ \frac{\partial u_3}{\partial x}& \frac{\partial u_3}{\partial y}& \frac{\partial u_3}{\partial z}\end{array}\right),$$

$$c(\mathbf{x})\otimes \nabla \mathbf{u}=\sum _{i,j=1}^{3}c(\mathbf{x}{)}_{ijkl}\nabla {\mathbf{u}}_{ij}.$$

$$\theta (\mathbf{x})=\{\begin{array}{ll}-{\text{dist}}_{\partial {\Omega }_1}(\mathbf{x})& \mathbf{x}\in {\Omega }_1\\ {\text{dist}}_{\partial {\Omega }_1}(\mathbf{x})& \mathbf{x}\in {\Omega }_2,\end{array}$$

$${\text{dist}}_E(\mathbf{x})=\underset{y\in E}{\inf }||\mathbf{x}-\mathbf{y}||\mathbf{x}\in \Omega .$$

$${\Omega }_1:=\{\mathbf{x}\in \Omega |\theta (\mathbf{x})<0\}$$

$${\Omega }_2:=\{\mathbf{x}\in \Omega |\theta (\mathbf{x})\ge 0\}.$$

$$\varphi (\theta )=\frac{1}{2}\left(\tanh \left(\frac{\theta }{\epsilon }\right)+1\right),$$

$$\rho (\theta )$$

$$c(\theta )$$

$$\frac{d\varphi (\theta )}{d\theta }=\frac{1}{2\epsilon }{\text{sech}}^2\left(\frac{\theta }{\epsilon }\right).$$

$$\frac{{\partial }^2\mathbf{u}}{\partial t^2}=\frac{\mathbf{u}-2{\mathbf{u}}_n+{\mathbf{u}}_{n-1}}{h^2}+O(h^2),\text{(as h→0)}$$

$$\frac{{\partial }^2\mathbf{u}}{\partial t^2}\approx \frac{\mathbf{u}-2{\mathbf{u}}_n+{\mathbf{u}}_{n-1}}{h^2},$$

$${\mathbf{u}}_{-1}={\mathbf{u}}_0-h{\mathbf{v}}_0.$$

$$\rho (\theta )\mathbf{u}=h^2\nabla \cdot (c(\theta )\otimes \nabla \mathbf{u})+h^2{\mathbf{f}}_n+\rho (\theta )(2{\mathbf{u}}_n-{\mathbf{u}}_{n-1}).$$

$${\int }_{\Omega }\rho (\theta )\mathbf{u}\cdot \mathbf{v}d\mathbf{x}$$

$$=$$

$$h^2{\int }_{\Omega }\nabla \cdot (c(\theta )\otimes \nabla \mathbf{u})\cdot \mathbf{v}d\mathbf{x}=-h^2{\int }_{\Omega }(c(\theta )\otimes \nabla \mathbf{u}):\nabla \mathbf{v}d\mathbf{x}.$$

$${\int }_{\Omega }\rho (\theta )\mathbf{u}\cdot \mathbf{v}d\mathbf{x}=-h^2{\int }_{\Omega }(c(\theta )\otimes \nabla \mathbf{u}):\nabla \mathbf{v}+h^{2}f\cdot \mathbf{v}+\rho (\theta )(2{\mathbf{u}}_n-{\mathbf{u}}_{n-1})\cdot \mathbf{v}d\mathbf{x},$$

$$\partial {\Omega }_1=\{\mathbf{x}\in {\mathbf{R}}^2:||\mathbf{x}-(0.5,0.5)||=0.1\}.$$

$${\chi }_{{\Omega }_2}(\mathbf{x})=\{\begin{array}{l}1\mathbf{x}\in {\Omega }_2\\ 0\text{otherwise.}\end{array}$$

$${\nu }_1=0.26,E_1=180\times 10^9$$

$${\nu }_2=0.25,E_2=70\times 10^9.$$

$${\rho }_1=4\times 10^3$$

$${\rho }_2=8\times 10^3.$$

$$r(\mathbf{x})=\sqrt{(x-l_x{)}^2+(y-l_y{)}^2)}\text{ and }\psi (\mathbf{x})=\arctan \left(\frac{y-l_y}{x-l_x}\right).$$

$$\mathbf{f}(\mathbf{x},t)=A\cos \left(\frac{2\pi t}{T}{l}_N\right)\exp (-l_{width}r)\left(\begin{array}{c}\cos \psi \\ \sin \psi \end{array}\right).$$

$$\text{dilation}\mathbf{u}=\frac{\mathrm{div}\mathbf{u}}{|\nabla \mathbf{u}|}.$$

$$\theta (\mathbf{x})=\{\begin{array}{ll}-{\text{dist}}_{{\Omega }_1}(\mathbf{x})& \mathbf{x}\in {\Omega }_1\\ {\text{dist}}_{{\Omega }_2}(\mathbf{x})& \mathbf{x}\in {\Omega }_2,\end{array}$$

$$\{\begin{array}{ll}\rho (\mathbf{x}){\mathbf{u}}_{tt}=\mathrm{div}\left(\sigma (\mathbf{x})\right)+\mathbf{f}(\mathbf{x},t)& \text{ in }(0,T)\times \Omega \\ u_i=g_i& \text{ on }(0,T)\times {\Gamma }_{g_i}\text{ (Dirichlet)}\\ {\sigma }_{ij}{n}_j=h_i& \text{ on }(0,T)\times {\Gamma }_{h_i}\text{ (Neumann)}\\ u_i(0,\mathbf{x})=u_{0i}(\mathbf{x})& \mathbf{x}\in \Omega \\ u_{i,t}(0,\mathbf{x})={\dot{u}}_{0i}(\mathbf{x})& \mathbf{x}\in \Omega \\ [u_i]=0& \text{on}\partial {\Omega }_1\cap \partial {\Omega }_{2}t\in (0,T)\\ [{\sigma }_{ij}{n}_j]=0& \text{on}\partial {\Omega }_1\cap \partial {\Omega }_{2}t\in (0,T)\end{array}$$

$$\sigma (\mathbf{x})=c(\mathbf{x}):ϵ(\mathbf{u}),$$

Full text

We also discuss a mathematical model for expressing the propagation of elastic waves through composite elastic bodies, and develop approximation schemes for investigating its numerical solutions. Using these methods, we define a cost functional for measuring the difference between simulated and given SAW data. Then, using a Lagrangian approach, we are able to determine the gradient of the cost functional and analyze the inverse imaging problem’s solution as a gradient flow. The cost functional’s gradient is composed of solutions to a state equation, as well as of solutions to related adjoint problems. We thus developed numerical methods for solving these problems and investigated the gradient flow of the cost functional. Our results show that the gradient flow is able to recover the interior composition of the composite, and we illustrate this fact using the numerical realization of our proposed framework.

1 Introduction

Designing accurate simulations of natural and social phenomena undoubtably requires the effective integration of quantitative and qualitative information. Such information is often gained through both theoretic and experimental methods, and techniques for its incorporation represent an important theme in data assimilation. It is also well known that relating theoretical assumptions and experimental observations to simulations often leads to the mathematical formulation of an inverse problem (see e.g., [6]). From CT (computed tomography) scanning techniques in the medical sciences, to sonic cloaking and sonar technology, inverse problems abound in applied mathematics and engineering.

In [7], the authors designed an experiment for obtaining surface acoustic wave (SAW) data of a phononic crystal and used k-space imaging techniques to analyze their experimental results. In addition to estimating band gaps of the phononic crystal, the authors showed that the crystal’s approximate surface geometry can be imaged using the acoustic field data. The current research investigates the possibility of using SAW data to image not only the exterior of composite media, but the interior as well. In particular, we design and analyze an inverse problem for imaging the interior of composite elastic materials, using SAW data.

This paper is outlined as follows. A mathematical model for describing elastic wave propagation through composite media is introduced in section 2. This section also discusses a level-set approach for describing the location of the different materials. Section 3 briefly describes the numerical approximation of solutions to the model equation, and these techniques are used in section 4 to illustrate the target phenomena’s simulation. The main target of this research, which is an inverse problem for describing the interior of composite elastic materials, is introduced in section 5. The inverse problem is formulated as a minimization between boundary (experimental) data and simulated data. We determine the gradient of the corresponding cost functional in section 6, and a level set method for approximating solutions of the inverse problem is formulated in section 7. The numerical results are shown in section 8, where we show that the gradient flow of the cost function indeed contains information on the interior of the composite media. We give our conclusions, including areas for future work, in section 9.

2 Modeling of wave propagation through composite elastic materials

This section describes the mathematical model equation used in simulating the propagation of elastic waves through composite elastic media. We also define a level set framework for expressing the location of material parameters. We will discuss our method for constructing numerical solutions to the model equations and show the numerical behavior of the model equation under an induced outer force.

We represent the composite media as a bounded, smooth, open subset $\Omega$ of ${\bf{R}^{3}}$ or ${\bf{R}^{2}}$ which is formed by a union:

[TABLE]

Here $\Omega_{1}$ and $\Omega_{2}$ are themselves smooth, open subsets of the ambient space which represent the two media (see figure 1 for a schematic representation in the two dimensional setting). In what follows, for the sake of clarity, we will assume the three dimensional setting, writing $\bm{x}=(x,y,z)$.

The elastic wave propagation is then assumed to be described by the displacement vector field $\bm{u}:[0,T)\times\Omega\rightarrow{\bf{R}^{3}}$, which evolves as a solution to the composite elastic wave equation (see e.g., [2] or [4]):

[TABLE]

When $\Omega$ is a subset of ${\bf{R}^{3}}$ the indices take the values $i=1,2,3$ and $j=1,2,3$. Similarly, whenever $\Omega$ is a subset of ${\bf{R}^{2}}$, the indices take the values $i=1,2$ and $j=1,2$. We also formally assume that the exterior of $\Omega$:

[TABLE]

is composed of regions $\Gamma_{g_{i}}$ and $\Gamma_{h_{i}}$, where Dirichlet and Neumann boundary conditions are prescribed, repectively (see figure, 2). In addition to the exterior boundary conditions on $\Omega$, interfacial boundary conditions holding on boundaries between $\Omega_{1}$ and $\Omega_{2}$ are prescribed. These are expressed as jump conditions, where $[\cdot]$ denotes the jump of its quantity across the normal to the interface (see e.g., [4])．

In the above, $\nabla\bm{u}$ refers to the matrix:

[TABLE]

$c(\bm{x})$ is a spatially dependent stiffness tensor, and $\rho(\bm{x})$ is a spatially dependent mass density function. The notation $c(\bm{x})\otimes\nabla\bm{u}$ denotes:

[TABLE]

The spatial dependence of $\rho(\bm{x})$ and $c(\bm{x})$ are described using a signed distance level set function:

[TABLE]

where $\text{dist}_{E}(\bm{x})$ is defined to be the distance from the point $\bm{x}$ to the set $E$:

[TABLE]

Here $||\cdot||$ denotes the euclidean norm. We remark that the level set function $\theta(\bm{x})$ allows us to determine the location of $\Omega_{1}$ and $\Omega_{2}$ as follows:

[TABLE]

The level set function is given as the argument of a sigmoid function (see figure 3):

[TABLE]

where $0<\varepsilon\ll 1$ is an interpolation parameter.

This allows us to control the mass density and elastic tensors using the interpolation function:

[TABLE]

where $\rho_{1},\rho_{2}$ and $c_{1},c_{2}$ correspond to the parameters within $\Omega_{1}$ and $\Omega_{2}$, respectively. For later use in our numerical computations, we remark that the derivative of the sigmoid function is an approximation of a delta function (see figure 3):

[TABLE]

The term $\bm{f}(\bm{x},t)$ in equation (2) denotes a time dependent outer force, whereas $u_{0i}(\bm{x})$ and $\dot{u}_{0i}(\bm{x})$ designate the initial displacement and velocity fields.

3 Approximation of solutions to the model equation

This section briefly discusses the methods used in constructing numerical solutions of the composite elastic wave equation (2).

Time is discretized using a time step $h>0$ and we have:

[TABLE]

where $\bm{u}_{n}$ and $\bm{u}_{n-1}$ are given vector valued functions approximating $\bm{u}(\bm{x},nh)$ and $\bm{u}(\bm{x},(n-1)h)$, respectively. That is, for suitably small $h>0,$ we have the approximation:

[TABLE]

where $n$ is a nonnegative integer. When $n=0$, $\bm{u}_{0}(\bm{x})$ denotes the initial displacement field and $\bm{u}_{-1}$ is constructed using the initial velocity field $\bm{v}_{0}(\bm{x})$ within a backward differencing:

[TABLE]

Inserting approximation (14) into the model equation (2) and rearranging yields the following elliptic partial differential equation for $\bm{u}(\bm{x})$:

[TABLE]

Here $\bm{f}_{n}(x)$ denotes $\bm{f}(\bm{x},nh)$. For prescribed boundary conditions, we then define $\bm{u}_{n}$ to be the solution of equation (16) and are thus able to approximate the solution of the model equation (2) as the sequence of functions $\{\bm{u}_{n}\}_{n=0}^{N},$ where $N$ denotes the final time step.

The numerical realization of each $\bm{u}_{n}$ is obtained as a weak solution of (16). In particular, assuming Dirichlet zero boundary conditions for simplicity, we multiply equation $(\ref{eq:sabun})$ by an arbitrary test function $\bm{v}\in C^{\infty}_{0}(\Omega;\mathbb{R}^{d})$ and integrate over $\Omega$:

[TABLE]

Applying Gauss-Green’s theorem (see e.g., [10]) and using the boundary conditions expresses

[TABLE]

The weak form of (16) is then expressed:

[TABLE]

for all test functions $\bm{v}$ in $C^{\infty}_{0}(\Omega;{\bf{R}}^{d})$.

Using the Lagrange $P^{1}$ finite element method (see e.g., [3]), we are thus able to obtain each $\{\bm{u}_{n}\}_{n=0}^{N}$ as the solution of a system of linear equations stemming from (19).

4 Numerical behavior of the model equation

Using the approximation method described in the previous section, we will now briefly illustrate the numerical behavior of solutions to the model equation (2) in the two dimensional setting. A Delaunay triangulation is employed on 5266 nodes to partition the domain $\Omega=(0,1)\times(0,1)$ into a finite number of elements, each with a target maximum edge length of 0.015 (figure 5 illustrates the mesh). We prescribe a time step approximately equal to $h=3\times 10^{-}6$.

The mass density function and stiffness tensors are prescribed using equation (2), where $\varepsilon=0.05$ and

[TABLE]

The corresponding level set function is given by equation (6) and figure 5 shows the truncation of the level set function:

[TABLE]

The purpose of the following computations is to show the influence of the interior structure on the wave propagation. For our first computation, the domain is assumed to be composed of two isotropic materials. For the sake of illustration, we prescribe $\Omega_{2}$ to represent the region of a glass and $\Omega_{1}$ to represent a steel region. Omitting units, the stiffness tensors are then expressed in terms of their Poisson’s ratio and Young’s modulus:

[TABLE]

Correspondingly, the density functions are set with the values:

[TABLE]

Elastic waves are generated by imparting an outer force $\bm{f}(\bm{x},t)$ near the top boundary of the domain, described as follows:

We also set

[TABLE]

The outer force used in our simulation is then defined (the vector field in figure 6 represents an instant in time of one such outer force):

[TABLE]

Using the approximation method described in the previous section, we solve the model equation (2) for $\bm{u}$ under Neumann boundary conditions, until approximately the time $T=4\times 10^{-4}$. The numerical results are shown in figure 7, where we display the dilation of the displacement field:

[TABLE]

Close inspection of figure 7 shows the spatial dependence imparted by the inclusion $\Omega_{1}$.

For comparison, we perform the same simulation where $\Omega$ is composed of entirely $\Omega_{2}$. The result is shown in figure 8, where we again display the dilation $\hat{\bm{u}}$. The difference between the dilations are shown in figure 9, where we clearly observe differences at the location of the inclusion.

5 The inverse problem

Having stated the model equation for elastic wave propagation through composite media and developed numerical methods for investigating its numerical simulation, we can now state the target inverse problem of this research. Moreover, for sake of agreement with the numerical investigations of the inverse problem, we will assume the case of isotropic materials. We remark that, in this setting, the model equation (2) becomes the composite linear elastic wave equation (see e.g., [4]).

As before, we assume that the domain $\Omega$ is composed of $\Omega_{1}$ and $\Omega_{2}$ and that their corresponding level set function has been constructed:

[TABLE]

where $\text{dist}_{E}(\bm{x})$ refers to the distance function (7) to the set $E$.

For a given layout of $\Omega_{1}$ and $\Omega_{2}$, we define a target data $\delta$ (defined below) obtained from the solution of the model equation:

[TABLE]

for given initial and boundary conditions, and where $\sigma(\bm{x})$ denotes the spatially dependent stress tensor:

[TABLE]

where $\epsilon(\bm{u})$ denotes the strain tensor:

[TABLE]

Here we remark that the density and stress tensor functions can be controlled through the level set function $\theta(\bm{x})$. The interfacial boundary condition $[u_{i}]=0$ specifies continuity in the displacement field, whereas $[\sigma_{ij}n_{j}]=0$ designates continuity of the stresses normal to the interface (see e.g., [4]). Here $n_{j}$ denotes the $j^{th}$ component of the normal to the interface, $\bm{n}$, and we assume that the outer force $\bm{f}(\bm{x},t)$ is a given time dependent vector field, such as equation (26).

Given a solution to the model equation (29), the boundary data $\delta$ is defined as a single component of the displacement field along an external boundary $\Gamma$ of $\Omega$:

[TABLE]

For example, when the closure of $\Omega$ is the closed unit cube $[0,1]\times[0,1]$, we can define $\Gamma:=\{\bm{x}\in\Omega\hskip 2.0pt|\hskip 2.0pty=1\}$ and $\delta(x,t)=u_{2}(x,1,t).$ Figure 10 shows the evolution of $\delta$ on the time interval $[0,9.0\times 10^{-4}$, where $\Omega_{1}$, $\Omega_{2}$, as well as their boundary conditions and material parameters are prescribed as in the simulations of the previous section (see equation (20)). Here we have prescribed $\bm{f}(\bm{x},t)$ using equation (26). We now state the target inverse problem of this study.

Given a boundary data $\delta:\Gamma\times(0,T)\rightarrow{\bf{R}}$ and a penalty parameter $\tilde{\epsilon}>0$, find $\theta(\bm{x})$ minimizing the cost functional

[TABLE]

where $\bm{u}(\bm{x},t)$ solves the model equation:

[TABLE]

Here the level set function $\theta({\bm{x}})$ controls the composition of $\Omega$, including the density function and elastic coefficients:

[TABLE]

6 The method of Lagrange multipliers for the inverse problem

In this section we derive the gradient of the cost functional $\mathcal{E}$ (see equation (33)) for investigating solutions of the target inverse problem (34). In particular, upon determining this gradient, we will investigate the solution of the inverse problem as the gradient flow of the cost functional. In order to do this, we will use a level set function $\theta(\bm{x})$ to describe the shape and location of the composite materials. This function will serve as the control parameter for minimizing $\mathcal{E}$. We will now employ the method of Lagrange multipliers, suitably adapted to the functional setting (see, e.g., [1], [5]).

6.1 Gradient of the cost functional $\mathcal{E}$

In order to use the method of Lagrange multipliers, we define the following Lagrangian:

[TABLE]

where $\bm{u}$ and $\bm{v}$ are arbitrary vector valued functions and $\theta(\bm{x})$ denotes the level set function (6). Here, $\mathcal{E}$ denotes the cost functional (33) for a given penalty parameter $\tilde{\epsilon}>0:$

[TABLE]

and $\mathcal{L}_{S}$ is expressed by the following equation:

[TABLE]

In the above, $T$ is a positive value corresponding to the final time in (2) and the vector valued function $\bm{f}(\bm{x},t)$ represents an outer force (for example, given by (26)). For $A,B\in\mathbb{R}^{n\times n}$, the symbol $:$ in $(c(\theta)\otimes\nabla\bm{u}):\nabla\bm{v}$ refers to the quantity:

[TABLE]

The Frechet derivative of the Lagrangian $\mathcal{L}$ is expressed:

[TABLE]

where $\vartheta,\tilde{\bm{u}},$ and $\tilde{\bm{v}}$ are arbitrary variations of $\theta,\bm{u}$, and $\bm{v}$, respectively.

The partial Frechet derivative of $\mathcal{L}$ with respect to $\bm{v}$ is expressed:

[TABLE]

A formal computation then yields:

[TABLE]

When $\bm{u}$ is a solution of equation (2), we remark that $\mathcal{L}_{S}(\theta,\bm{u},\bm{v})$ takes the value 0 for any arbitrary $\bm{v}$.

Similarly, the partial Frechet derivative of $\mathcal{L}$ with respect to $\bm{u}$ is computed:

[TABLE]

Here we remark that $\mathcal{L}_{\bm{u}}(\theta,\bm{u},\bm{v})[\tilde{\bm{u}}]$ takes the value 0 whenever $\bm{v}$ is a solution of the following adjoint problem:

[TABLE]

In this way, having set $\bm{u}$, and $\bm{v}$, and since $\vartheta$ is an arbitrary variation of $\theta$, we can obtain the gradient of the cost functional as the remaining term in equation (39):

[TABLE]

Here $\rho^{\prime}(\theta),c^{\prime}(\theta)$ are regularizations of delta functions supported on interfaces between the regions $\Omega_{1}$ and $\Omega_{2}$:

[TABLE]

In the above, $\varepsilon$ is a positive parameter which prescribes the cost functional’s gradient to act as a regularized delta function near the interface. We remark that the particular value of $\varepsilon$ is important when performing numerical calculations.

Having determined the gradient of the cost functional, we now turn to designing an approximation method for solving the inverse problem.

7 An approximation method for the inverse problem

In this section, using the cost functional’s gradient information obtain by the method of Lagrange multipliers from the previous section, we will develop an approximation scheme for investigating the numerical behavior of solutions to our target inverse problem. A brief explanation of our method is as follows. We first construct the boundary data $\delta$ for a given target composition. Second, we set the composition of the region to be different from the target and solve the corresponding model equation. This gives us the displacement field arising in the adjoint problem (41). After solving the adjoint problem we are then able to construct the gradient of the cost functional. This gradient information is then used within the level set method to evolve the interface and decrease the value of the cost functional. The details are explained in the following.

7.1 Algorithm for approximating solutions of the inverse problem

The numerical solution of the inverse problem is expressed as a gradient flow of the cost functional (33), which is performed as follows.

Step 1: Construct the boundary data $\delta$.

Set the target composition of $\Omega$ and construct its level set function:

[TABLE]

where $\Omega=\tilde{\Omega}_{1}\cup\tilde{\Omega}_{2}$ denotes the target domain composition. 2. 2.

Solve the model equation (47) corresponding to $\theta_{target}(\bm{x})$:

[TABLE]

for a given outer force $\bm{f}(\bm{x},t)$, and initial and boundary conditions. 3. 3.

Extract the boundary data from the solution:

[TABLE]

where $\Gamma$ denotes an exterior boundary of $\Omega$.

Step 2: Set the initial condition for the level set method.

[TABLE]

Here, $\Omega_{1}$ and $\Omega_{2}$ are generally assumed to be different from the target compositions $\tilde{\Omega}_{1}$ and $\tilde{\Omega}_{2}$.

Step 3: Construct the cost functional gradient of $\mathcal{E}$ with respect to $\theta(\bm{x})$.

Solve the model equation (50) corresponding to $\theta_{init}(\bm{x})$:

[TABLE] 2. 2.

Solve the adjoint problem $(\ref{eq:optim3})$, using $u_{2}(\bm{x},t)$ and $\delta$:

[TABLE] 3. 3.

Construct the cost functional gradient of $\mathcal{E}$ as the function:

[TABLE]

where $\rho(\theta_{init})$ and $c(\theta_{init})$ are given by equation (2).

Step 4: Evolve the level set function using the gradient $g(\bm{x})$ within the level set equation.

Solve the level set equation for a time $\tau^{*}>0$:

[TABLE] 2. 2.

Evolve the level set function:

[TABLE]

Step 5: Repeat steps 3 and 4 until $\mathcal{E}$ converges.

7.2 Properties of the cost functional gradient

As we have explained, the gradient of the cost functional depends on various parameters through solutions to partial differential equations. For this reason, we have performed numerical experiments to confirm its basic properties. In particular, in this section, we will show that the cost functional gradient indeed holds information on the location of inclusions. In particular, using the domain triangulation and parameters set forth in section 4, we set set the initial condition for the level set method as in equation (20) and investigate the behavior of $g(\bm{x})$ for seven different target compositions.

We employ an immersed boundary technique in the adjoint equation for $\bm{v}$. In particular, by setting Neumann zero boundary conditions on the boundary of $\Omega$, our approach models the Dirichlet boundary condition on $\Gamma$ by using a regularized delta function as an outer force:

[TABLE]

where we set $\tilde{\epsilon}=1/1000$, and $\epsilon=1/25$.

Figure 11 shows the cost functional gradient $g$ responding to various target compositions. The first row of the figure shows the initial condition (the light blue region) overlaid with each target (shown in dark blue). The second row shows the corresponding gradients.

Six cases were considered. In the first (figure 11 a), the $\theta_{target}$ is placed above $\theta_{init}$ and near the upper boundary. Inspection of the gradient $g$ confirms that the level set equation evolves the interface with a normal velocity in the direction of the target. As will be shown in the next section, evolving the level set function by the level set equation decreases both the distance between the target and the value of the cost functional. Figure 12 displays the value of the cost functional $\mathcal{E}(\theta_{init}-\tau g)$ for $\tau$ within the interval $[-0.01,0.01]$ where we are able to verify that the gradient descent of $\mathcal{E}$ is in the direction $-g(\bm{x})$.

8 Numerical solution of the inverse problem

This section examines the numerical solution of the inverse problem by observing the gradient flow of the cost functional. We begin by employing our approximation method on each of the five configurations shown in figure 11 a-e.

The level set function $\theta_{init}(\bm{x})$, the outer force $\bm{f}(\bm{x},t)$, and the remaining parameters are set as in the previous sections. We set $\theta_{target}(\bm{x})$ (a disk with the same radius as $\theta_{init}(\bm{x})$, shown as a light shadow in the figures) directly above the initial condition and near the boundary of the computational domain (see figure 11 a). Repeatedly solving the level set equation (53) yields the evolution shown in figure 13. We have displayed twelve instants of the evolution, where we observe the level set equation’s solutions to converge to the prescribed target. Figure 14 shows a plot of the cost functional’s value as a function of the $n^{th}$ iteration of the approximation method. We confirm the value of the cost functional to decrease along the flow.

Our next computation sets $\theta_{target}(\bm{x})$ directly below the initial condition, near the bottom boundary (see figure 11 e). We apply our approximation method and observe the motion of the level set function as it converges to its target. Figure 15 shows 12 instants in time of the evolution, and figure 16 shows the value of the cost functional at each iteration step. We again confirm that the gradient flow decreases the value of the cost functional.

A similar result can be obtained when the target is located to the right of the initial condition (see figure 11 c). Our method yields an approximation of the gradient flow, shown in figure 17. We again observe the values of the cost functional to decrease with each iteration of the optimization (figure 18).

We also examine the gradient flow of the cost functional when the target is placed in the upper right corner of the computational domain (as in figure 11 b). The evolution is shown in figure 19. The level set function converges to the approximate location of $\theta_{target}$ while decreasing the value of the cost functional (see figure 20). Although the shape of the limit function is not as precise as those obtained in our the previous computations, we remark that its overall shape and location is correct. We also recall that these results are obtained through the use of the boundary data $\delta$, and nothing more. Boundary conditions could be effecting the computational result, but we also expect that mesh and parameter dependence may be playing a role as well. Similar results are obtained when the target is place in the bottom right corner of the computational domain (see figures 21 and 22).

We also considered the case where $\theta_{init}$ overlaps $\theta_{target}$ with a larger radius. The numerical results are shown at six instances in figure 23. Here, we observe the level set function to shrink inward while converging to its target. The evolution of the cost functional’s values under the gradient flow are displayed in figure 24 which shows a linear decrease.

Our final numerical experiment sets the target as two disconnected disks with radii equal to that of $\theta_{init}$ (see figure 25). The numerical results show the level set function spread towards the target and finally splits into two disconnected regions. Although the cost functional’s value decreases with each iteration of our approximation method, we note that the limit function is visibly different from the target. Nevertheless, the limit function visually does represent an approximation of the target and this fact suggests the possible existence of local minimizers of the cost functional near the target.

The computations in this section clearly show the ability of our method to approximate solutions of the inverse problem stated in section 5. Overall, these results show that the cost functional’s gradient indeed contains information about the interior composition of the domain.

9 Conclusion and future work

A mathematical technique for determining the internal structure of composite elastic bodies was developed. We showed how a single component of the displacement field along a portion of the boundary of the elastic body can be used to gain information about the interior composition. This was done by using the method of Lagrange multipliers to derive the gradient of a cost functional whose minimization corresponds to our target inverse problem.

Numerical solvers were also developed for computing solutions to model and adjoint equations. These were used to construct numerical realizations of the cost functional’s gradient. In addition, we designed a level set method which incorporates the cost functional’s gradient to evolve the level set function as the gradient descent of the cost functional. This allowed us to examine the numerical solution of the inverse problem and our results showed that our framework is able to recover the interior composition of composite elastic materials, using only a single component of the displacement field’s boundary data.

We would like to continue improving our framework for use on problems using real experimental data, such as the surface acoustic wave data mentioned in the introduction (see also [7]). Here we are particularly interested in whether or not our method can be used to image the interior of composite elastic bodies which contain sub-micron scale inclusions. Relatedly, we aim to develop numerical methods for treating the three dimensional phenomena as well. Since the numerical methods developed in this study tend to be computationally heavy, we would like to improve their efficiency. Finally, we are also planning to investigate inverse problems involving a larger number of composite materials.

10 Acknowledgements

The work of E. Ginder was supported by JST Presto “Multiphase shape optimization in phononic crystal design.” The author would also like to acknowledge support from JSPS KAKENHI Grant Number 17K14229 and JSPS KAKENHI Grant Number 15KT0099.