Inverse diffraction problems for buried objects in the layered elastic media: the anti-plane problem

Abstract

The authors study direct and inverse scattering problems for buried objects in the stratified layered elastic media. It is assumed that a source load is applied on the upper boundary surface, to generate a certain wave field inside the medium. Then we build the solutions to the direct problems through specially constructed Green’s functions for every considered case. Using the known measured amplitude of the upper boundary oscillations over a certain finite-length interval, the inverse problems in a one-mode frequency range are solved. Some numerical examples demonstrate the stability of the proposed algorithm for the formulated inverse problems.

KEYWORDS:

• Buried object
• layered elastic medium
• inverse diffraction problem
• boundary integral equations
• object identification

AMS SUBJECT CLASSIFICATIONS:

• 45B05
• 45Q05
• 74B05

1. Introduction

Recently, there was a need to study layered solids which contain buried objects, to identify the geometries of the objects and to test the methods optimal for detection of such objects. An example of important problem of this kind is given by a black-box detection of damages in the aircrafts using the scattering of the electromagnetic waves. Another example is related to location and size recognition for defects in the elastic layers, for instance, the detection of cavities or voids in structures using ultrasonic waves. There are also other applications, such as mines detection and medical diagnostics. Mathematically, the problems in these fields are reduced to some boundary value problems for elliptic partial differential equations. The Laplace or the Helmholtz equations give the simplest examples of the governing equation for such problems, describing wave processes in elastic and acoustic media.

It is well-known that the boundary value problem for any linear isotropic elastic volume of arbitrary configuration can be reduced to a system of boundary integral equations (BIE), using Green’s function or Green’s tensor; this allows one to evaluate the displacement wave field in the considered domains. The works [1–7] give the governing equations and some efficient techniques for the direct and the inverse diffraction problems in the acoustic and elastic cases. Those methods are applied to simplest geometries of the objects to be identified as well as to simplest geometries of the domain under detection. The application of mathematical and ultrasonic methods to object identification in the elastic media is demonstrated in [8–11]. Very often the input data for the reconstruction in such inverse problems is the measured amplitude of oscillations over a certain part of the outer boundary [5]. The identification method for several cracks simultaneously by the ultrasonic (US) scanning is proposed in [12]. Some applications to the engineering practice of the Ultrasonic (US) impulse propagation in the micro-damaged materials can be found in [13].

The choice of the appropriate mathematical and numerical techniques to identify the geometrical objects, such as cracks, cavities or inclusions, is the subject of many investigations, see [2,14] where the reader can also find further helpful references. The most common approach to solve an inverse problem is its reduction to a system of non-linear operator equations, as a rule, of the first kind. Therefore, generally, the inverse problems are ill-posed, in the sense that small perturbations of the input data may cause large change of the solution. The numerical solution for these operator equations require a certain regularization method [15]. If the geometry of the object under identification is of a particular form defined by a finite number of parameters, like in the cases of linear cracks, spherical or other particular geometries, then the inversion of respective first-kind operators is studied on a finite-dimensional subset. It is known that in this case the inversion is well-posed and numerically stable. The main goal of the present work is to develop a stable numerical algorithm, rather than to prove existence and uniqueness of the formulated problems.

The aim of the present work is to study the inverse problems on reconstruction of buried objects in an isotropic elastic layered medium, based on the known fields measured at a certain part of the free boundary of the medium at a fixed frequency. Using the BIE method, the direct problems are reduced to singular integral equations over the boundary contour of the object. The solution of the BIE system is obtained using the numerical collocation method. Moreover, the problem is reduced to a system of linear algebraic equations with respect to displacement function at the nodes on the object boundary contour. This allows the authors to find the whole field in the medium, in particular, on the upper free surface of the medium. Therefore, the solution of the direct problem, perturbed by a random noise, may be used as the input data for the inverse problem.

The solution to the inverse problems, about identification of a buried object in the stratified elastic media using the input data on the displacement field measured over a certain part of the free surface of the medium is constructed in two steps: in the first step, we construct the dynamic displacement field in the single layer or in the case of two layers, by the known source load on the upper surface of the medium. While in the second step and based on this latter information, the inverse problem for identification of the geometry of the buried object and its location is solved. The consideration is restricted by elliptic voids, in the one-mode frequency range. This predetermines a finite-dimensional inversion, since any elliptic object is defined by no more that five geometric parameters. More specifically, the inverse problem for domain’s geometry described by a finite number of parameters can directly be reduced to a linearization of the well-posed initial system of non-linear equations, which automatically leads to a stable identification algorithm.

2. The governing equations for a buried object in the case of single- and double-layered elastic medium

2.1. Problem1

Here, we formulate the anti-plane (SH) problem of harmonic oscillations of an isotropic elastic layer containing a void. Let us consider the homogeneous elastic layer $|y_1|<\infty ,0<y_2<h$ of thickness h containing a buried void with a smooth contour l, see Figure 1. Let us choose the Cartesian coordinates $(y_1,y_2)$ such that the horizontal axis coincides with the lower fixed boundary $y_2=0$, and the upper free boundary $y_2=h$ of the layer is loaded by a point force applied at $y_1=x_0$ and harmonically oscillating with the circular frequency $\mathit{\omega }$.

Figure 1. A buried object in the single-layered elastic medium.

In the case of anti-plane (SH) deformation, the only nontrivial component of the displacement vector field is $u_3=w(y_1,y_2)e^{-i\mathit{\omega }t}$. The only nontrivial components of the stress tensor are: ${\mathit{\sigma }}_{13}=\mathit{\mu }\mathit{\partial }w/\mathit{\partial }{y}_1$, ${\mathit{\sigma }}_{23}=\mathit{\mu }\mathit{\partial }w/\mathit{\partial }{y}_2$, where $\mathit{\mu }$ is the shear elastic modulus. Then the problem is reduced to the following boundary value Problem 1, see [1]:(2.1) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_{y}w(y_1,y_2)+k^{2}w(y_1,y_2)=0,\left(k=\frac{\mathit{\omega }}{c},c^2=\frac{\mathit{\mu }}{\mathit{\rho }}\right),\end{array}$$(2.1) (2.2) $$\begin{array}{c}\hfill {\left.\frac{\mathit{\partial }w}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=\frac{p_0}{\mathit{\mu }}\mathit{\delta }(y_1-x_0),{\left.w\right|}_{y_2=0}=0,{\left.\frac{\mathit{\partial }w}{\mathit{\partial }{n}_y}\right|}_l=0,\end{array}$$(2.2)

where ${\mathrm{\Delta }}_y={\mathit{\partial }}^2/\mathit{\partial }{y}_1^2+{\mathit{\partial }}^2/\mathit{\partial }{y}_2^2$ is the Laplace operator, k is the wave number, c is the transverse wave speed, $\mathit{\rho }$ is the mass density, $p_0$ is the amplitude of the outer applied oscillating point force. To complete the formulation, it is necessary to add a radiation condition at infinity.

If the geometry of the contour l is known then the problem is a standard direct diffraction problem. If it is solved by any method, the basic wave field $w(y_1,y_2)$ can be determined at any observation point, including the upper free boundary surface $y_2=h$. If the void’s boundary contour l is unknown then one can formulate respective inverse problem to identify the geometry of the contour by knowing, say the amplitude of the oscillations over a finite interval of the free surface $y_2=h$. The latter can be taken as an input data from the results of the ‘measurements’ obtained from some natural or numerical experiments. In our theoretical treatment such an input data may be obtained from the numerical solution of respective direct problem subjected to a certain perturbation by a random noise.

Further, let us apply a standard technique of the potential theory. If one constructs Green’s function $G(y,\mathit{\xi })$,   $\mathit{\xi }=({\mathit{\xi }}_1,{\mathit{\xi }}_2)$, inside the layer(2.3) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_y{G}^{(1)}+k^2{G}^{(1)}=-\mathit{\delta }(y_1-{\mathit{\xi }}_1)\mathit{\delta }(y_2-{\mathit{\xi }}_2),\end{array}$$(2.3)

which satisfies, similarly to the two first boundary conditions (2.2(2.2) $$\begin{array}{c}\hfill {\left.\frac{\mathit{\partial }w}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=\frac{p_0}{\mathit{\mu }}\mathit{\delta }(y_1-x_0),{\left.w\right|}_{y_2=0}=0,{\left.\frac{\mathit{\partial }w}{\mathit{\partial }{n}_y}\right|}_l=0,\end{array}$$(2.2) ), the Neumann homogeneous condition over the upper face of the layer $y_2=h$ and the Dirichlet homogeneous condition over the lower face of the layer $y_2=0$:(2.4) $$\begin{array}{c}\hfill {\left.\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=0,{\left.G^{(1)}\right|}_{y_2=0}=0,\end{array}$$(2.4)

then a classical result of the potential theory claims that [16](2.5) $$\begin{array}{c}\hfill w(\mathit{\xi })={\int }_{\mathit{\partial }D}\left[w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}-G^{(1)}(y,\mathit{\xi })\frac{\mathit{\partial }w(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\right]\mathrm{d}{l}_y,\end{array}$$(2.5)

for any point $\mathit{\xi }$ inside the layer located outside contour l . The superscript (1) here above Green’s function $G^{(1)}(y,\mathit{\xi })$ indicates that the latter is related to Problem 1, $n_y$ is the outer unit normal to the contour $\mathit{\partial }D$, which itself is a union of contour l and the two faces of the layer: $\mathit{\partial }D=\{y_2=0\}\bigcup \{y_2=h\}\bigcup l$. Obviously, the integral over the lower layer’s face is trivial, since both the quantities w and $G^{(1)}$ are trivial over this horizontal axis, due to their boundary conditions. Analogously, the integral in (2.5(2.5) $$\begin{array}{c}\hfill w(\mathit{\xi })={\int }_{\mathit{\partial }D}\left[w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}-G^{(1)}(y,\mathit{\xi })\frac{\mathit{\partial }w(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\right]\mathrm{d}{l}_y,\end{array}$$(2.5) ) taken over the upper layer’s face $y_2=h$ can also be calculated explicitly, using respective boundary conditions for quantities $\mathit{\partial }w/\mathit{\partial }{n}_y$ and $\mathit{\partial }{G}^{(1)}/\mathit{\partial }{n}_y$. As a result, Equation (2.5(2.5) $$\begin{array}{c}\hfill w(\mathit{\xi })={\int }_{\mathit{\partial }D}\left[w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}-G^{(1)}(y,\mathit{\xi })\frac{\mathit{\partial }w(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\right]\mathrm{d}{l}_y,\end{array}$$(2.5) ) can be rewritten in the following form:(2.6) $$\begin{array}{cc}\hfill {\displaystyle w(\mathit{\xi })}& =f(\mathit{\xi })+{\int }_{l}w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y,\hfill \\ \hfill {\displaystyle f(\mathit{\xi })}& =\frac{p_0}{\mathit{\mu }}{\left.G^{(1)}(y,\mathit{\xi })\right|}_{y_1=x_0,y_2=h}=\frac{p_0}{\mathit{\mu }}{G}^{(1)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),\hfill \end{array}$$(2.6)

where f is the wave field in the homogeneous layer, free of any buried object. Note that ${\overline{n}}_y$ in Equation (2.6(2.6) $$\begin{array}{cc}\hfill {\displaystyle w(\mathit{\xi })}& =f(\mathit{\xi })+{\int }_{l}w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y,\hfill \\ \hfill {\displaystyle f(\mathit{\xi })}& =\frac{p_0}{\mathit{\mu }}{\left.G^{(1)}(y,\mathit{\xi })\right|}_{y_1=x_0,y_2=h}=\frac{p_0}{\mathit{\mu }}{G}^{(1)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),\hfill \end{array}$$(2.6) ) is the unit normal vector chosen outer to contour l. Now, one can easily deduce from Equation (2.6(2.6) $$\begin{array}{cc}\hfill {\displaystyle w(\mathit{\xi })}& =f(\mathit{\xi })+{\int }_{l}w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y,\hfill \\ \hfill {\displaystyle f(\mathit{\xi })}& =\frac{p_0}{\mathit{\mu }}{\left.G^{(1)}(y,\mathit{\xi })\right|}_{y_1=x_0,y_2=h}=\frac{p_0}{\mathit{\mu }}{G}^{(1)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),\hfill \end{array}$$(2.6) ) the basic boundary integral equation (BIE) of Problem 1, if one moves the point $\mathit{\xi }$ to the boundary l : $\mathit{\xi }\to x=(x_1,x_2)\in l$, from outside. Then, according to the properties of the potential of double layer [2], one deduces from Equation (2.6(2.6) $$\begin{array}{cc}\hfill {\displaystyle w(\mathit{\xi })}& =f(\mathit{\xi })+{\int }_{l}w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y,\hfill \\ \hfill {\displaystyle f(\mathit{\xi })}& =\frac{p_0}{\mathit{\mu }}{\left.G^{(1)}(y,\mathit{\xi })\right|}_{y_1=x_0,y_2=h}=\frac{p_0}{\mathit{\mu }}{G}^{(1)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),\hfill \end{array}$$(2.6) ) the following BIE:(2.7) $$\begin{array}{cc}\hfill {\displaystyle \frac{w(x)}{2}-{\int }_{l}K(y,x)w(y)\mathrm{d}{l}_y}& =f(x),x\in l,\hfill \\ \hfill {\displaystyle K(y,x)}& =\frac{\mathit{\partial }{G}^{(1)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}^{(1)}·n_y\}=n_y^1\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}.\hfill \end{array}$$(2.7)

It should be noted that generally the kernels in some formulations of BIEs may possess certain singularities [2]. However, for the anti-plane problem, it is known that in the case of smooth boundary line the Neumann boundary value problem is reduced to the second-kind integral equation whose kernel contains no singularity as $\mathit{\xi }\to x$. It can be proved that in this case the limiting value of the kernel is expressed in terms of curvature of the boundary curve which is always finite if the contour is smooth.

2.2. Problem 2

Let us consider the Problem 2 about harmonic oscillations of two homogeneous elastic layers of thickness $h_1,h_2$, respectively, see Figure 2. The cavity with a smooth contour l is located in the upper layer (I), while the lower layer (II) consists of a different elastic material. Assume that the oscillations in the media is caused again, like in the Problem 1, by a source load applied at the upper surface $y_2=h=h_1+h_2$ of the first layer, i.e. at a certain point $(x_0,h)$. It should be noted that the case when the defect is located in the lower layer (II) can be studied by analogy.

Figure 2. A buried object in the upper layer of the double-layered elastic medium.

The Problem 2 at hand is mathematically governed by the following equations. For the upper layer (I), it can be written as follows:(2.8) $$\begin{array}{c}\hfill {\displaystyle {\mathrm{\Delta }}_y{w}_1(y_1,y_2)+k_1^2{w}_1(y_1,y_2)=0,k_1=\frac{\mathit{\omega }}{c_1},}\\ \hfill {\displaystyle {\left.{\mathit{\mu }}_1\frac{\mathit{\partial }{w}_1}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=p_0\mathit{\delta }(y_1-x_0),{\left.\frac{\mathit{\partial }{w}_1}{\mathit{\partial }{n}_y}\right|}_l=0.}\end{array}$$(2.8)

For the lower layer (II) one writes:(2.9) $$\begin{array}{cc}& {\displaystyle {\mathrm{\Delta }}_y{w}_2(y_1,y_2)+k_2^2{w}_2(y_1,y_2)=0,k_2=\frac{\mathit{\omega }}{c_2},}\hfill \\ \hfill & {\displaystyle {\left.w_2\right|}_{y_2=0}=0,}\hfill \end{array}$$(2.9)

Besides, the boundary conditions on the interface $y_2=h_2$ should be added since the displacement and the stress fields are continuous when crossing this interface line:(2.10) $$\begin{array}{c}\hfill {\displaystyle {\left.w_1\right|}_{y_2=h_2}={\left.w_2\right|}_{y_2=h_2},{\left.{\mathit{\mu }}_1\frac{\mathit{\partial }{w}_1}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}={\left.{\mathit{\mu }}_2\frac{\mathit{\partial }{w}_2}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}.}\end{array}$$(2.10)

In the same way as in Problem 1, a specific Green’s function can be introduced for this boundary value Problem 2 ($h=h_1+h_2$):(2.11) $$\begin{array}{c}\hfill {\displaystyle G^{(2)}(y_1,y_2)=\left\{\begin{array}{c}{G}_1^{(2)}(y_1,y_2),h_2<y<h,\\ G_2^{(2)}(y_1,y_2),0<y<h_2.\end{array}\right.}\end{array}$$(2.11)

Green’s function in the upper layer (I) is defined as follows: (2.12a) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_y{G}_1^{(2)}(y,\mathit{\xi })+k_1^2{G}_1^{(2)}(y,\mathit{\xi })=-{\mathit{\delta }(y}_1-{\mathit{\xi }}_1){\mathit{\delta }(y}_2-{\mathit{\xi }}_2),{\left.\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=0,\end{array}$$(2.12a)

while for the lower layer:(2.12b) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_y{G}_2^{(2)}(y,\mathit{\xi })+k_2^2{G}_2^{(2)}(y,\mathit{\xi })=0,{\left.G_2^{(2)}\right|}_{y_2=0}=0,\end{array}$$(2.12b) with the same boundary conditions of the continuity on the interface line $y_2=h_2$:(2.13) $$\begin{array}{c}\hfill {\left.G_1^{(2)}\right|}_{y_2=h_2}={\left.G_2^{(2)}\right|}_{y_2=h_2},{\mathit{\mu }}_1{\left.\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}={\mathit{\mu }}_2{\left.\frac{\mathit{\partial }{G}_2^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}.\end{array}$$(2.13)

Using the constructed Green’s function, it can be proved that the following representation of the potential theory, like in (2.6(2.6) $$\begin{array}{cc}\hfill {\displaystyle w(\mathit{\xi })}& =f(\mathit{\xi })+{\int }_{l}w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y,\hfill \\ \hfill {\displaystyle f(\mathit{\xi })}& =\frac{p_0}{\mathit{\mu }}{\left.G^{(1)}(y,\mathit{\xi })\right|}_{y_1=x_0,y_2=h}=\frac{p_0}{\mathit{\mu }}{G}^{(1)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),\hfill \end{array}$$(2.6) ):(2.14) $$\begin{array}{c}\hfill {\displaystyle w_1(\mathit{\xi })=f(\mathit{\xi })+{\int }_l{w}_1(y)\frac{\mathit{\partial }{G}_1^{(2)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}l,f(\mathit{\xi })=\frac{p_0}{{\mathit{\mu }}_1}{G}_1^{(2)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),}\end{array}$$(2.14)

is valid at arbitrary point $\mathit{\xi }=({\mathit{\xi }}_1,{\mathit{\xi }}_2)$ in the upper layer, outside contour l. The standard approach leads again to the basic BIE in this Problem 2, analogous to Equation (2.7(2.7) $$\begin{array}{cc}\hfill {\displaystyle \frac{w(x)}{2}-{\int }_{l}K(y,x)w(y)\mathrm{d}{l}_y}& =f(x),x\in l,\hfill \\ \hfill {\displaystyle K(y,x)}& =\frac{\mathit{\partial }{G}^{(1)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}^{(1)}·n_y\}=n_y^1\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}.\hfill \end{array}$$(2.7) ) in Problem 1:(2.15) $$\begin{array}{cc}& {\displaystyle \frac{w_1(x)}{2}-{\int }_{l}K(y,x)w_1(y)\mathrm{d}{l}_y=f(x),x\in l,}\hfill \\ \hfill & {\displaystyle K(y,x)=\frac{\mathit{\partial }{G}_1^{(2)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}_1^{(2)}.n_y\}=n_y^1\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}.}\hfill \end{array}$$(2.15)

2.3. Problem 3

Further, let us consider the anti-plane Problem 3 about harmonic oscillations of two homogeneous elastic layers of thickness $h_1$ and $h_2$, respectively. On the interface boundary surface $y_2=h_2$ there is located a void whose smooth boundary contour $l=l_1\bigcup l_2$ is free of any load, see Figure 3. It is assumed that the upper part of the contour, $l_1$ is fully located in the first layer (I), and the lower part $l_2$ is fully located in the second layer (II). Therefore, the boundary value problem 3 can be written for each layer as follows:

Figure 3. A buried object on the interface of the double-layered elastic medium.

For the upper layer (I): (2.16a) $$\begin{array}{cc}& {\displaystyle {\mathrm{\Delta }}_y{w}_1(y_1,y_2)+k_1^2{w}_1(y_1,y_2)=0,}\hfill \\ \hfill & {\displaystyle {\left.{\mathit{\mu }}_1\frac{\mathit{\partial }{w}_1}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=p_0\mathit{\delta }(y_1-x_0),{\left.\frac{\mathit{\partial }{w}_1}{\mathit{\partial }{n}_y}\right|}_{l_1}=0.}\hfill \end{array}$$(2.16a)

For the lower layer (II):(2.16b) $$\begin{array}{c}\hfill {\displaystyle {\mathrm{\Delta }}_y{w}_2(y_1,y_2)+k_2^2{w}_2(y_1,y_2)=0,}\\ \hfill {\displaystyle {\left.w_2\right|}_{y_2=0}=0,{\left.\frac{\mathit{\partial }{w}_2}{\mathit{\partial }{n}_y}\right|}_{l_2}=0.}\end{array}$$(2.16b)

To complete the formulation of the Problem 3, one should add the boundary conditions (2.10(2.10) $$\begin{array}{c}\hfill {\displaystyle {\left.w_1\right|}_{y_2=h_2}={\left.w_2\right|}_{y_2=h_2},{\left.{\mathit{\mu }}_1\frac{\mathit{\partial }{w}_1}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}={\left.{\mathit{\mu }}_2\frac{\mathit{\partial }{w}_2}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}.}\end{array}$$(2.10) ) on the internal interface $y_2=h_2$, outside the buried object.

To solve the formulated boundary value Problem 3, we should build Green’s tensor in the two-layered medium. The components of Green’s tensor can be presented as follows:(2.17) $$\begin{array}{c}\hfill G^{(3)}(y,\mathit{\xi })=\left\{\begin{array}{c}{G}_{11}(y,\mathit{\xi }),\mathit{\xi }\in D_1,y\in D_1\\ G_{12}(y,\mathit{\xi }),\mathit{\xi }\in D_1,y\in D_2\\ G_{21}(y,\mathit{\xi }),\mathit{\xi }\in D_2,y\in D_1\\ G_{22}(y,\mathit{\xi }),\mathit{\xi }\in D_2,y\in D_2\end{array}\right.,\sim G^{(3)}=\left(\begin{array}{c}{G}_{11}{G}_{12}\\ G_{21}{G}_{22}\end{array}\right),\end{array}$$(2.17)

where domain $D_1$ contains all points of the first layer outside contour $l_1$, and domain $D_2$ is analogous in the second layer.

Thus, the pair of the boundary value problems can be formulated here: (2.18a) $$\begin{array}{c}\hfill {\displaystyle {\displaystyle {\mathrm{\Delta }}_y\left(\begin{array}{c}{G}_{11}\\ G_{12}\end{array}\right)+k_{1,2}^2\left(\begin{array}{c}{G}_{11}\\ G_{12}\end{array}\right)=\left\{\begin{array}{c}-{\mathit{\delta }(y}_1-{\mathit{\xi }}_1){\mathit{\delta }(y}_2-{\mathit{\xi }}_2),\mathit{\xi },y\in D_1,\\ {\displaystyle 0,\mathit{\xi }\in D_1,y\in D_2,}\end{array}\right.}}\\ \hfill {\displaystyle {\left.G_{12}\right|}_{y_2=0}=0,{\left.\frac{\mathit{\partial }{G}_{11}}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=0,(h=h_1+h_2),}\end{array}$$(2.18a)

and(2.18b) $$\begin{array}{cc}& {\mathrm{\Delta }}_y\left(\begin{array}{c}{G}_{21}\\ G_{22}\end{array}\right)+k_{1,2}^2\left(\begin{array}{c}{G}_{21}\\ G_{22}\end{array}\right)=\left\{\begin{array}{c}0,\mathit{\xi }\in D_2,y\in D_1,\\ {\displaystyle -{\mathit{\delta }(y}_1-{\mathit{\xi }}_1){\mathit{\delta }(y}_2-{\mathit{\xi }}_2),\mathit{\xi },y\in D_2,}\end{array}\right.\hfill \\ \hfill & {\displaystyle {\left.G_{22}\right|}_{y_2=0}=0,{\left.\frac{\mathit{\partial }{G}_{21}}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=0.}\hfill \end{array}$$(2.18b)

The boundary conditions on the internal interface$y_2=h_2$, are the conditions of continuity for the component of Green’s displacements tensor and their corresponding components of shearing stress that:(2.19) $$\begin{array}{c}\hfill G_{i1}=G_{i2},{\mathit{\mu }}_1\frac{\mathit{\partial }{G}_{i1}}{\mathit{\partial }{y}_2}={\mathit{\mu }}_2\frac{\mathit{\partial }{G}_{i2}}{\mathit{\partial }{y}_2},(\mathit{\xi }\in D_i,i=1,2).\end{array}$$(2.19)

By analogy to Problems 1 and 2, if Green’s tensor is constructed for the present Problem 3, the standard approach of the potential theory reduces this problem to the system of BIEs: (2.20a) & (2.20b) $$\begin{array}{cc}& {\displaystyle x\in l_1:\frac{{\mathit{\mu }}_1}{2}{w}_1(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{11}(x),f_{11}(\mathit{\xi })=p_0{G}_{11}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),}\hfill \\ \hfill & {\displaystyle x\in l_2:\frac{{\mathit{\mu }}_2}{2}{w}_2(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{22}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{21}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{21}(x),f_{21}(\mathit{\xi })=p_0{G}_{21}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2).}\hfill \end{array}$$(2.20a) & (2.20b)

In order to solve the direct Problems 1–3, one should build Green’s functions for the corresponding boundary value problems. The next section is devoted to construct Green’s functions and Green’s tensors for elastic layered medium in an explicit form.

3. Green’s functions for Problems 1–3

3.1. Problem 1

Let us start with Problem 1 for the single layer; this means to solve the boundary value problem (2.3(2.3) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_y{G}^{(1)}+k^2{G}^{(1)}=-\mathit{\delta }(y_1-{\mathit{\xi }}_1)\mathit{\delta }(y_2-{\mathit{\xi }}_2),\end{array}$$(2.3) )–(2.4(2.4) $$\begin{array}{c}\hfill {\left.\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=0,{\left.G^{(1)}\right|}_{y_2=0}=0,\end{array}$$(2.4) ).

The solution to Equation (2.3(2.3) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_y{G}^{(1)}+k^2{G}^{(1)}=-\mathit{\delta }(y_1-{\mathit{\xi }}_1)\mathit{\delta }(y_2-{\mathit{\xi }}_2),\end{array}$$(2.3) ), for any fixed point $\mathit{\xi }=({\mathit{\xi }}_1,{\mathit{\xi }}_2)$ inside the layer, may be obtained using the integral Fourier transform along the horizontal axis $y_1$:(3.1) $$\begin{array}{cc}\hfill {\displaystyle {\tilde{G}}^{(1)}(s,y_2)}& ={\int }_{-\infty }^{\infty }{G}^{(1)}(y_1,y_2)e^{is{y}_1}d{y}_1,\hfill \\ \hfill {\displaystyle G^{(1)}(y_1,y_2)}& =\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{\tilde{G}}^{(1)}(s,y_2)e^{-is{y}_1}\mathrm{d}s.\hfill \end{array}$$(3.1)

Then, using the classical property of Dirac’s delta function:$e^{is{\mathit{\xi }}_1}={\int }_{-\infty }^{\infty }\mathit{\delta }(y_1-{\mathit{\xi }}_1)e^{is{y}_1}d{y}_1$, the general solution to homogeneous equation and a particular solution of the full equation, one deduces:(3.2) $$\begin{array}{c}\hfill {\tilde{G}}^{(1)}(s,y_2)=A(s)e^{\mathit{\gamma }{y}_2}+B(s)e^{-\mathit{\gamma }{y}_2}+\frac{e^{is{\mathit{\xi }}_1}{e}^{-\mathit{\gamma }\left|{\mathit{\xi }}_2-y_2\right|}}{2\mathit{\gamma }},\mathit{\gamma }(s)=\sqrt{s^2-k^2},\end{array}$$(3.2)

the quantities $A(s),B(s)$should be determined from the boundary conditions (2.4(2.4) $$\begin{array}{c}\hfill {\left.\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=0,{\left.G^{(1)}\right|}_{y_2=0}=0,\end{array}$$(2.4) ). The solution is constructed in the following form [14]: (3.3a) & (3.3b) $$\begin{array}{cc}\hfill G^{(1)}(y_1,y_2)& =\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{L}^{(1)}(s,y_2)e^{-is(y_1-{\mathit{\xi }}_1)}\mathrm{d}s,\hfill \\ \hfill L^{(1)}(s,y_2)& =\frac{sinh[\mathit{\gamma }(h-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[\mathit{\gamma }({\mathit{\xi }}_2+y_2-h)]}{2\mathit{\gamma }cosh(\mathit{\gamma }h)}.\hfill \end{array}$$(3.3a) & (3.3b) Notice that function $L^{(1)}(s,x_2)$ is meromorphic with respect to s. The (simple) poles of function $L^{(1)}(s,x_2)$ coincide with the zeros of the denominator:(3.4) $$\begin{array}{c}\hfill s_m=\pm i{\mathit{\sigma }}_m,{\mathit{\sigma }}_m=\sqrt{{\left(\frac{{\mathit{\pi }}_m}{h}\right)}^2-k^2},{\mathit{\pi }}_m=\mathit{\pi }(m+1/2),m=0,1,2,...\end{array}$$(3.4)

It is clear that for sufficiently low frequencies (small values of the wave number k) real-valued poles are absent, then all singularities of the integrand in (3.3a(3.3a) & (3.3b) $$\begin{array}{cc}\hfill G^{(1)}(y_1,y_2)& =\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{L}^{(1)}(s,y_2)e^{-is(y_1-{\mathit{\xi }}_1)}\mathrm{d}s,\hfill \\ \hfill L^{(1)}(s,y_2)& =\frac{sinh[\mathit{\gamma }(h-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[\mathit{\gamma }({\mathit{\xi }}_2+y_2-h)]}{2\mathit{\gamma }cosh(\mathit{\gamma }h)}.\hfill \end{array}$$(3.3a) & (3.3b) ) are located outside the contour of integration. This is the case we call ‘the one-mode regime’. However, with frequency increasing more and more poles cross the real axis $\mathrm{Im}(s)=0$. It can be proved that Sommerfeld’s principle of radiation condition predetermines that the integration path should bend a little around the positive poles from below and the negative ones from above. In practice, in the demonstrated examples of object’s identification, see Section 6 below, the value of the wave number is chosen at the top part of the considered one-mode frequency range.

By applying the Jordan lemma [17], the integral (3.3(3.3a) & (3.3b) $$\begin{array}{cc}\hfill G^{(1)}(y_1,y_2)& =\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{L}^{(1)}(s,y_2)e^{-is(y_1-{\mathit{\xi }}_1)}\mathrm{d}s,\hfill \\ \hfill L^{(1)}(s,y_2)& =\frac{sinh[\mathit{\gamma }(h-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[\mathit{\gamma }({\mathit{\xi }}_2+y_2-h)]}{2\mathit{\gamma }cosh(\mathit{\gamma }h)}.\hfill \end{array}$$(3.3a) & (3.3b) ) can be calculated as a series by residues at simple poles:(3.5) $$\begin{array}{c}\hfill G^{(1)}(y_1,y_2)=\sum _{m=0}^{\infty }\left\{cos\left[\frac{{\mathit{\pi }}_m(y_2-{\mathit{\xi }}_2)}{h}\right]-cos\left[\frac{{\mathit{\pi }}_m(y_2+{\mathit{\xi }}_2)}{h}\right]\right\}\frac{e^{-{\mathit{\sigma }}_m\left|y_1-x_1\right|}}{2h{\mathit{\sigma }}_m}.\end{array}$$(3.5)

3.2. Problem 2

Let us pass to the Green’s functions for Problem 2, which is to be determined as a solution of the boundary value problem (2.11(2.11) $$\begin{array}{c}\hfill {\displaystyle G^{(2)}(y_1,y_2)=\left\{\begin{array}{c}{G}_1^{(2)}(y_1,y_2),h_2<y<h,\\ G_2^{(2)}(y_1,y_2),0<y<h_2.\end{array}\right.}\end{array}$$(2.11) )–(2.13(2.13) $$\begin{array}{c}\hfill {\left.G_1^{(2)}\right|}_{y_2=h_2}={\left.G_2^{(2)}\right|}_{y_2=h_2},{\mathit{\mu }}_1{\left.\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}={\mathit{\mu }}_2{\left.\frac{\mathit{\partial }{G}_2^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}.\end{array}$$(2.13) ). By analogy to Problem 1 we apply the Fourier transform along variable $y_1$, which reduces the governing equations to a pair of ordinary differential equations, with respect to variable $y_2$ whose solution in Fourier images is:(3.6) $$\begin{array}{cc}\hfill {\displaystyle {\tilde{G}}_1^{(2)}(s,y_2)}& =A_1(s)e^{{\mathit{\gamma }}_1{y}_2}+B_1(s)e^{-{\mathit{\gamma }}_1{y}_2}+\frac{e^{is{\mathit{\xi }}_1}{e}^{-{\mathit{\gamma }}_1\left|{\mathit{\xi }}_2-y_2\right|}}{2{\mathit{\gamma }}_1},\hfill \\ \hfill {\displaystyle {\tilde{G}}_2^{(2)}(s,y_2)}& =A_2(s)e^{{\mathit{\gamma }}_2{y}_2}+B_2(s)e^{-{\mathit{\gamma }}_2{y}_2},{\mathit{\gamma }}_{1,2}={\mathit{\gamma }}_{1,2}(s)=\sqrt{s^2-k_{1,2}^2}.\hfill \end{array}$$(3.6)

Then the full representation for the Green’s functions under consideration is given by the inverse Fourier transform, as follows:(3.7) $$\begin{array}{c}\hfill G_{1,2}^{(2)}(y,\mathit{\xi })=\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{L}_{1,2}^{(2)}(s,y_2)e^{-is(y_1-{\mathit{\xi }}_1)}\mathrm{d}s.\end{array}$$(3.7)

The substitution of representations (3.6(3.6) $$\begin{array}{cc}\hfill {\displaystyle {\tilde{G}}_1^{(2)}(s,y_2)}& =A_1(s)e^{{\mathit{\gamma }}_1{y}_2}+B_1(s)e^{-{\mathit{\gamma }}_1{y}_2}+\frac{e^{is{\mathit{\xi }}_1}{e}^{-{\mathit{\gamma }}_1\left|{\mathit{\xi }}_2-y_2\right|}}{2{\mathit{\gamma }}_1},\hfill \\ \hfill {\displaystyle {\tilde{G}}_2^{(2)}(s,y_2)}& =A_2(s)e^{{\mathit{\gamma }}_2{y}_2}+B_2(s)e^{-{\mathit{\gamma }}_2{y}_2},{\mathit{\gamma }}_{1,2}={\mathit{\gamma }}_{1,2}(s)=\sqrt{s^2-k_{1,2}^2}.\hfill \end{array}$$(3.6) ) into the corresponding boundary conditions (2.12(2.12a) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_y{G}_1^{(2)}(y,\mathit{\xi })+k_1^2{G}_1^{(2)}(y,\mathit{\xi })=-{\mathit{\delta }(y}_1-{\mathit{\xi }}_1){\mathit{\delta }(y}_2-{\mathit{\xi }}_2),{\left.\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h}=0,\end{array}$$(2.12a) ), (2.13(2.13) $$\begin{array}{c}\hfill {\left.G_1^{(2)}\right|}_{y_2=h_2}={\left.G_2^{(2)}\right|}_{y_2=h_2},{\mathit{\mu }}_1{\left.\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}={\mathit{\mu }}_2{\left.\frac{\mathit{\partial }{G}_2^{(2)}}{\mathit{\partial }{y}_2}\right|}_{y_2=h_2}.\end{array}$$(2.13) ) leads to a $4\times 4$ SLAE for the four unknown quantities $A_1,B_1,A_2,B_2$. Its solution determines functions $L_{1,2}^{(2)}(s,y_2)$ in the following form ($h=h_1+h_2,H=h_1+2h_2$): (3.8a) & (3.8b) $$\begin{array}{cc}\hfill {\displaystyle L_1^{(2)}}& =\frac{{\mathit{\mu }}_{1}sinh({\mathit{\gamma }}_2{h}_2)\left\{cosh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}+\hfill \\ \hfill {\displaystyle }& +\frac{{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_2{h}_2)\left\{sinh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2{\mathit{\gamma }}_1\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}},\hfill \\ \hfill L_2^{(2)}& =\frac{{\mathit{\mu }}_{1}cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2-h)]sinh({\mathit{\gamma }}_2{y}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)}.\hfill \end{array}$$(3.8a) & (3.8b)

Obviously, Green’s functions for Problem 2 should degenerate to the one in the Problem 1. One can easily control that the two following limiting cases take place:

(i)

If the thickness of the second layer tends to zero then $h_2\to 0$, $h_1=h$, ${\mathit{\gamma }}_1=\mathit{\gamma }$, and representation (3.8a(3.8a) & (3.8b) $$\begin{array}{cc}\hfill {\displaystyle L_1^{(2)}}& =\frac{{\mathit{\mu }}_{1}sinh({\mathit{\gamma }}_2{h}_2)\left\{cosh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}+\hfill \\ \hfill {\displaystyle }& +\frac{{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_2{h}_2)\left\{sinh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2{\mathit{\gamma }}_1\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}},\hfill \\ \hfill L_2^{(2)}& =\frac{{\mathit{\mu }}_{1}cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2-h)]sinh({\mathit{\gamma }}_2{y}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)}.\hfill \end{array}$$(3.8a) & (3.8b) ) degenerates to (3.3b(3.8a) & (3.8b) $$\begin{array}{cc}\hfill {\displaystyle L_1^{(2)}}& =\frac{{\mathit{\mu }}_{1}sinh({\mathit{\gamma }}_2{h}_2)\left\{cosh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}+\hfill \\ \hfill {\displaystyle }& +\frac{{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_2{h}_2)\left\{sinh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2{\mathit{\gamma }}_1\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}},\hfill \\ \hfill L_2^{(2)}& =\frac{{\mathit{\mu }}_{1}cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2-h)]sinh({\mathit{\gamma }}_2{y}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)}.\hfill \end{array}$$(3.8a) & (3.8b) ).

(ii)

If the materials of the two elastic layers are identical then ${\mathit{\mu }}_1={\mathit{\mu }}_2=\mathit{\mu }$, ${\mathit{\gamma }}_1={\mathit{\gamma }}_2=\mathit{\gamma }$. Let point $\mathit{\xi }$ belong to the upper layer (I) and point y belong to the lower layer (II), then obviously $y_2<{\mathit{\xi }}_2$ in Equation (3.8a(3.8a) & (3.8b) $$\begin{array}{cc}\hfill {\displaystyle L_1^{(2)}}& =\frac{{\mathit{\mu }}_{1}sinh({\mathit{\gamma }}_2{h}_2)\left\{cosh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}+\hfill \\ \hfill {\displaystyle }& +\frac{{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_2{h}_2)\left\{sinh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2{\mathit{\gamma }}_1\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}},\hfill \\ \hfill L_2^{(2)}& =\frac{{\mathit{\mu }}_{1}cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2-h)]sinh({\mathit{\gamma }}_2{y}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)}.\hfill \end{array}$$(3.8a) & (3.8b) ). It is evident that in this case functions (3.8a(3.8a) & (3.8b) $$\begin{array}{cc}\hfill {\displaystyle L_1^{(2)}}& =\frac{{\mathit{\mu }}_{1}sinh({\mathit{\gamma }}_2{h}_2)\left\{cosh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}+\hfill \\ \hfill {\displaystyle }& +\frac{{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_2{h}_2)\left\{sinh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2{\mathit{\gamma }}_1\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}},\hfill \\ \hfill L_2^{(2)}& =\frac{{\mathit{\mu }}_{1}cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2-h)]sinh({\mathit{\gamma }}_2{y}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)}.\hfill \end{array}$$(3.8a) & (3.8b) ) and (3.8b(3.8a) & (3.8b) $$\begin{array}{cc}\hfill {\displaystyle L_1^{(2)}}& =\frac{{\mathit{\mu }}_{1}sinh({\mathit{\gamma }}_2{h}_2)\left\{cosh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}+\hfill \\ \hfill {\displaystyle }& +\frac{{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_2{h}_2)\left\{sinh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2{\mathit{\gamma }}_1\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}},\hfill \\ \hfill L_2^{(2)}& =\frac{{\mathit{\mu }}_{1}cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2-h)]sinh({\mathit{\gamma }}_2{y}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)}.\hfill \end{array}$$(3.8a) & (3.8b) ) should coincide with each other. It can easily be proved that in this case (3.9) $$\begin{array}{c}\hfill L_1^{(2)}(s,y_2)=L_2^{(2)}(s,y_2)=\frac{cosh[\mathit{\gamma }(h-{\mathit{\xi }}_2)]sinh(\mathit{\gamma }{y}_2)}{\mathit{\gamma }cosh(\mathit{\gamma }h)},\end{array}$$(3.9) that is to be proved.

Like in Problem 1, functions $L_{1,2}^{(2)}(s,y_2)$ are meromorphic with respect to variable s. However, it is difficult to find the zeros of the denominator in Equations (3.8(3.8a) & (3.8b) $$\begin{array}{cc}\hfill {\displaystyle L_1^{(2)}}& =\frac{{\mathit{\mu }}_{1}sinh({\mathit{\gamma }}_2{h}_2)\left\{cosh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}+\hfill \\ \hfill {\displaystyle }& +\frac{{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_2{h}_2)\left\{sinh[{\mathit{\gamma }}_1(h_1-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[{\mathit{\gamma }}_1({\mathit{\xi }}_2+y_2-H)]\right\}}{2{\mathit{\gamma }}_1\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}},\hfill \\ \hfill L_2^{(2)}& =\frac{{\mathit{\mu }}_{1}cosh[{\mathit{\gamma }}_1({\mathit{\xi }}_2-h)]sinh({\mathit{\gamma }}_2{y}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)}.\hfill \end{array}$$(3.8a) & (3.8b) ) analytically. Therefore, there is no efficient series representation of integral (3.7(3.7) $$\begin{array}{c}\hfill G_{1,2}^{(2)}(y,\mathit{\xi })=\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{L}_{1,2}^{(2)}(s,y_2)e^{-is(y_1-{\mathit{\xi }}_1)}\mathrm{d}s.\end{array}$$(3.7) ), the latter should be treated directly by an appropriate numerical quadrature formula.

3.3. Problem 3

At last, let us develop the Green’s tensor for Problem 3. By writing all four functions in terms of the inverse Fourier transform(3.10) $$\begin{array}{c}\hfill G_{\mathit{mj}}(y,\mathit{\xi })=\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{L}_{\mathit{mj}}(s,y_2)e^{-is(y_1-{\mathit{\xi }}_1)}\mathrm{d}s,(m,j=1,2),\end{array}$$(3.10)

it can be easily seen, on the basis of the governing equations, that the first two of them coincide with respective functions of Problem 2. Namely, one obtains: $L_{11}(s,y_2)=L_1^{(2)}(s,y_2)$, $L_{12}(s,y_2)=L_2^{(2)}(s,y_2)$. The other two components are determined in the following way: (3.11a) & (3.11b) $$\begin{array}{cc}\hfill L_{21}& =\frac{{\mathit{\mu }}_{2}cosh[{\mathit{\gamma }}_1(h-y_2)]cosh({\mathit{\gamma }}_2{\mathit{\xi }}_2)}{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)},\hfill \\ \hfill {\displaystyle L_{22}}& =\frac{{\mathit{\mu }}_{2}cosh({\mathit{\gamma }}_1{h}_1)\left\{cosh[{\mathit{\gamma }}_2({\mathit{\xi }}_2+y_2-h_2)]-cosh[{\mathit{\gamma }}_2(h_2-\left|{\mathit{\xi }}_2-y_2\right|)]\right\}}{2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}\hfill \\ \hfill & {\displaystyle -\frac{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)\left\{sinh[{\mathit{\gamma }}_2({\mathit{\xi }}_2+y_2-h_2)]+sinh[{\mathit{\gamma }}_2(h_2-\left|{\mathit{\xi }}_2-y_2\right|]\right\}}{2{\mathit{\gamma }}_2\left\{{\mathit{\mu }}_1{\mathit{\gamma }}_{1}sinh({\mathit{\gamma }}_1{h}_1)sinh({\mathit{\gamma }}_2{h}_2)+{\mathit{\mu }}_2{\mathit{\gamma }}_{2}cosh({\mathit{\gamma }}_1{h}_1)cosh({\mathit{\gamma }}_2{h}_2)\right\}}.}\hfill \end{array}$$(3.11a) & (3.11b)

4. Numerical treatment of the BIE

In order to solve the direct diffraction Problems 1–3  formulated above, one should solve integral Equations (2.7(2.7) $$\begin{array}{cc}\hfill {\displaystyle \frac{w(x)}{2}-{\int }_{l}K(y,x)w(y)\mathrm{d}{l}_y}& =f(x),x\in l,\hfill \\ \hfill {\displaystyle K(y,x)}& =\frac{\mathit{\partial }{G}^{(1)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}^{(1)}·n_y\}=n_y^1\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}.\hfill \end{array}$$(2.7) ), (2.15(2.15) $$\begin{array}{cc}& {\displaystyle \frac{w_1(x)}{2}-{\int }_{l}K(y,x)w_1(y)\mathrm{d}{l}_y=f(x),x\in l,}\hfill \\ \hfill & {\displaystyle K(y,x)=\frac{\mathit{\partial }{G}_1^{(2)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}_1^{(2)}.n_y\}=n_y^1\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}.}\hfill \end{array}$$(2.15) ), (2.20(2.20a) & (2.20b) $$\begin{array}{cc}& {\displaystyle x\in l_1:\frac{{\mathit{\mu }}_1}{2}{w}_1(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{11}(x),f_{11}(\mathit{\xi })=p_0{G}_{11}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),}\hfill \\ \hfill & {\displaystyle x\in l_2:\frac{{\mathit{\mu }}_2}{2}{w}_2(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{22}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{21}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{21}(x),f_{21}(\mathit{\xi })=p_0{G}_{21}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2).}\hfill \end{array}$$(2.20a) & (2.20b) ), respectively. Let us apply a collocation numerical technique in its simplest form. This implies to split the full closed contour l to N small chords of a polyline, then to apply a simple quadrature formula with a constant integrand over each small chord $l_j,j=1,2,\cdots ,N$ of length$\mathrm{\Delta }{l}_j$. Note that point $x=(x_1,x_2)$ runs over the same set of nodes as point $y=(y_1,y_2)$: $x_j=y_j\in l_j$  , $(j=1,2,\cdots ,N)$. As noticed above, the arising BIEs of the second kind in the case of smooth boundary contour l have no singularity as $y\to x$, hence the collocation method can be treated like in the regular cases with continuous kernels [18]. With the simplest treatment it is demonstrated here for Problem 1, Equation (2.7(2.7) $$\begin{array}{cc}\hfill {\displaystyle \frac{w(x)}{2}-{\int }_{l}K(y,x)w(y)\mathrm{d}{l}_y}& =f(x),x\in l,\hfill \\ \hfill {\displaystyle K(y,x)}& =\frac{\mathit{\partial }{G}^{(1)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}^{(1)}·n_y\}=n_y^1\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}.\hfill \end{array}$$(2.7) ):(4.1) $$\begin{array}{c}\hfill \frac{w_m}{2}-\sum _{j=1}^{N}K(x_j,x_m)\mathrm{\Delta }{l}_j{w}_j=f_m,m=1,\cdots ,N,\end{array}$$(4.1)

where $w_m=w(x_m),f_m=f(x_m)$. The case of Equation (2.15(2.15) $$\begin{array}{cc}& {\displaystyle \frac{w_1(x)}{2}-{\int }_{l}K(y,x)w_1(y)\mathrm{d}{l}_y=f(x),x\in l,}\hfill \\ \hfill & {\displaystyle K(y,x)=\frac{\mathit{\partial }{G}_1^{(2)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}_1^{(2)}.n_y\}=n_y^1\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}.}\hfill \end{array}$$(2.15) ), as well as the system of integral Equations (2.20(2.20a) & (2.20b) $$\begin{array}{cc}& {\displaystyle x\in l_1:\frac{{\mathit{\mu }}_1}{2}{w}_1(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{11}(x),f_{11}(\mathit{\xi })=p_0{G}_{11}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),}\hfill \\ \hfill & {\displaystyle x\in l_2:\frac{{\mathit{\mu }}_2}{2}{w}_2(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{22}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{21}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{21}(x),f_{21}(\mathit{\xi })=p_0{G}_{21}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2).}\hfill \end{array}$$(2.20a) & (2.20b) ), can be treated by analogy.

The only unclear question here is connected with the known property that a vertex of an angle, where two different elastic materials meet, may contain a singularity of the stress. In fact, this may happen in Problem 3 at two points where the boundary curve l intersects the interface line $y_2=h_2$, the situation which is locally like a vertex of the angle mentioned above.

Let us show that in the anti-plane problem this singularity is certainly absent at least in the case when the boundary contour l is locally orthogonal to the interface line $y_2=h_2$. In this case, we have two quarter-planes joined by one face, each of them with the second face free of load. The proof is based on the property that dynamic elastic problems possess the same singularities as the respective static problems. In the static case the governing equation is the Laplace equation, which should hold inside the upper (I) and the lower (II) quarter-plane. In the local polar coordinate system $(r,\mathit{\theta })$, $0<r<\infty ,-\mathit{\pi }/2<\mathit{\theta }<0$ for domain (II) and $0<r<\infty ,0<\mathit{\theta }<\mathit{\pi }/2$ for domain (I), this means:(4.2) $$\begin{array}{c}\hfill \frac{{\mathit{\partial }}^2{w}_1}{\mathit{\partial }{r}^2}+\frac{1}{r}\frac{\mathit{\partial }{w}_1}{\mathit{\partial }r}+\frac{1}{r^2}\frac{{\mathit{\partial }}^2{w}_1}{\mathit{\partial }{\mathit{\theta }}^2}=0,\frac{{\mathit{\partial }}^2{w}_2}{\mathit{\partial }{r}^2}+\frac{1}{r}\frac{\mathit{\partial }{w}_2}{\mathit{\partial }r}+\frac{1}{r^2}\frac{{\mathit{\partial }}^2{w}_2}{\mathit{\partial }{\mathit{\theta }}^2}=0.\end{array}$$(4.2)

Let us seek the solution to the ‘pair-of-quarter-planes’ problem in the following form:(4.3) $$\begin{array}{c}\hfill w_1=r^{{\mathit{\alpha }}_1}{f}_1(\mathit{\theta }),w_2=r^{{\mathit{\alpha }}_2}{f}_2(\mathit{\theta }),{\mathit{\alpha }}_1={\mathit{\alpha }}_2=\mathit{\alpha },\end{array}$$(4.3)

with a certain possible singularity as $r\to 0$, in the case if $\mathit{\alpha }$ is not a positive integer. Here, due to a symmetry (the opening angles in both the quarter-planes are equal to $\mathit{\pi }/2$), we set equal singularity power: ${\mathit{\alpha }}_1={\mathit{\alpha }}_2$. The substitution of (3.3(3.3a) & (3.3b) $$\begin{array}{cc}\hfill G^{(1)}(y_1,y_2)& =\frac{1}{2\mathit{\pi }}{\int }_{-\infty }^{\infty }{L}^{(1)}(s,y_2)e^{-is(y_1-{\mathit{\xi }}_1)}\mathrm{d}s,\hfill \\ \hfill L^{(1)}(s,y_2)& =\frac{sinh[\mathit{\gamma }(h-\left|y_2-{\mathit{\xi }}_2\right|)]+sinh[\mathit{\gamma }({\mathit{\xi }}_2+y_2-h)]}{2\mathit{\gamma }cosh(\mathit{\gamma }h)}.\hfill \end{array}$$(3.3a) & (3.3b) ) into (3.2(3.2) $$\begin{array}{c}\hfill {\tilde{G}}^{(1)}(s,y_2)=A(s)e^{\mathit{\gamma }{y}_2}+B(s)e^{-\mathit{\gamma }{y}_2}+\frac{e^{is{\mathit{\xi }}_1}{e}^{-\mathit{\gamma }\left|{\mathit{\xi }}_2-y_2\right|}}{2\mathit{\gamma }},\mathit{\gamma }(s)=\sqrt{s^2-k^2},\end{array}$$(3.2) ) leads to the ordinary differential equations $f_{1,2}^{″}(\mathit{\theta })+{\mathit{\alpha }}^2{f}_{1,2}=0$, whose general solution is:(4.4) $$\begin{array}{c}\hfill f_1(\mathit{\theta })=A_{1}cos(\mathit{\alpha }\mathit{\theta })+A_{2}sin(\mathit{\alpha }\mathit{\theta }),f_2(\mathit{\theta })=A_{3}cos(\mathit{\alpha }\mathit{\theta })+A_{4}sin(\mathit{\alpha }\mathit{\theta }).\end{array}$$(4.4)

Now, the boundary conditions on the interface $\mathit{\theta }=0:w_1=w_2,{\mathit{\mu }}_1\mathit{\partial }{w}_1/\mathit{\partial }\mathit{\theta }={\mathit{\mu }}_2\mathit{\partial }{w}_2/\mathit{\partial }\mathit{\theta }$, as well as over free side surfaces, $\mathit{\theta }=\mathit{\pi }/2:\mathit{\partial }{w}_1/\mathit{\partial }\mathit{\theta }=0$; $\mathit{\theta }=-\mathit{\pi }/2:\mathit{\partial }{w}_2/\mathit{\partial }\mathit{\theta }=0$, lead to a $4\times 4$ homogeneous SLAE for the unknown coefficients $A_1-A_4$. Since we seek a non-trivial solution of this system, its determinant should be trivial, which results in the simple transcendental equation: $sin(\mathit{\pi }\mathit{\alpha })=0$. Therefore, it is clear that appropriate values of this transcendental equation are $\mathit{\alpha }=0,1,2,3,\cdots$, since negative integers break the finite-energy condition. It follows from the obtained result that there is no singularity of the solution for such a simple geometry.

It should be noted that the case of ‘no singularity’ studied above is automatically attained if the centre of the elliptic object is located on the interface line $y_2=h_2$, since in this case the tangent lines at the two intersection points are locally orthogonal to the interface, i.e. just the case studied above in this section.

All examples presented below are related to the outer harmonic force of a unit amplitude $p_0/\mathit{\mu }$. Figure 4 demonstrates an example of the solution to the direct Problem 1 for the considered single-layered medium, which contains an elliptic cavity as a buried object. This is calculated for various inclination angles $\mathit{\theta }$ between the principal axis of the ellipses and the horizontal axis $x_1$. The presented diagrams are related to the displacement field over the upper free surface of the layer: $w(\mathit{\zeta },h)=\overline{w}(\mathit{\zeta })$. The variation of the mesh dimension from $N=20$ to $N=70$ shows insignificant variation of the diagrams, invisible in the figure. Similar diagrams for a round object are presented in Figure 5, for various horizontal coordinates of the object’s central point, quantity $x_1^c$. Figures 6 and 7 show in the direct Problem 3 the diagram of the displacement over the boundary contour l for an elliptic object of a number of different aspect ratio. This is plotted separately for the upper part of contour l, which is located fully in the upper layer (I), and for the lower part of contour l, which is fully located in the lower layer (II). The relation of elastic parameters between the two elastic media is the same as copper (medium I) and steel (medium II).

Figure 4. Direct Problem 1, $h=2.5$, $k=0.5$, $x_0=2h$. The displacement over the upper surface, elliptic void: semi-axes $a=0.2h,b=0.1h$, central point $(x_1^c,x_2^c)=(-0.4h,0.4h)$, $\mathit{\theta }$ is the slope of the principal axis with respect to $x_1$.

Figure 5. Direct Problem 1, $h=2.5$, $k=0.5$, $x_0=2h$. The displacement over the upper surface, round void: $a=b=0.1h$, $x_2^c=0.4h$.

Figure 6. Direct Problem 3, $h=1.0$, $h_1=h_2=h/2$, $k_1=1.0$, $k_2=0.4$, $x_0=h/2$. The displacement $w_1$ over the upper part of the boundary contour l versus polar angle of the current boundary point on the interval $(0,\mathit{\pi })$, elliptic void: semi-axes $a=0.05,b$, central point $(x_1^c,x_2^c)=(h/2,h/2)$.

Figure 7. Direct Problem 3, $h=1.0$, $h_1=h_2=h/2$, $k_1=1.0$, $k_2=0.4$, $x_0=h/2$. The displacement $w_2$ over the lower part of the boundary contour l versus polar angle of the current boundary point on the interval $(\mathit{\pi },2\mathit{\pi })$, elliptic void: semi-axes $a,b=0.05$, central point $(x_1^c,x_2^c)=(h/2,h/2)$.

Figure 8. Direct Problem 1, $h=4.0$, $k=0.36$, $x_0=h/2$, $n=4$. The displacement w over the boundary line versus polar angle, three different rectangular voids with filleted corners: $(x_1^c,x_2^c)=(h/2,h/2)$.

Additionally, in order to demonstrate the possibilities of the proposed method, we plot in Figure 8 the displacement w at the points of the boundary contour l, direct problem 1, for three different rectangular voids with filleted corners. The defining equation is as follows:(4.5) $$\begin{array}{cc}\hfill {\displaystyle x_1(\mathit{\psi })}& =x_1^c+\mathit{\rho }(\mathit{\psi })cos\mathit{\psi },x_2(\mathit{\psi })=x_2^c+\mathit{\rho }(\mathit{\psi })sin\mathit{\psi },\hfill \\ \hfill {\displaystyle \mathit{\rho }(\mathit{\psi })}& ={\left[{\left(\frac{cos\mathit{\psi }}{a}\right)}^n+{\left(\frac{sin\mathit{\psi }}{b}\right)}^n\right]}^{-1/n},0\le \mathit{\psi }\le 2\mathit{\pi }.\hfill \end{array}$$(4.5)

5. Formulation of the inverse Problems 1–3

The inverse problem for any respective direct Problem 1–3 is to restore, or to identify, the location and the geometry of the buried object, under assumption that this is of an elliptic form. The formulation is as follows. Under condition of the harmonic oscillation of the layered medium generated by the outer point force applied over the upper free surface $y_2=h$ at point $y_1=x_0$, the information about the amplitude of the free surface oscillations, function $w(\mathit{\zeta },h)=\overline{w}(\mathit{\zeta })$ is known at a fixed frequency, over a certain interval $\mathit{\zeta }\in [c,d],y_2=h$. It is assumed that $x_0\notin [c,d]$. In order to solve the inverse problem, the input data, function $\overline{w}(\mathit{\zeta })$, is taken from the solution of the respective direct problem, perturbed by a random noise of a certain amplitude. As can be seen from the previous sections, function $\overline{w}(\mathit{\zeta })$ is easily calculated for any problem if respective BIEs are solved numerically by the collocation method. More precisely,

In Problem 1:(5.1) $$\begin{array}{c}\hfill \overline{w}(\mathit{\zeta })={\int }_{l}w(y)\frac{\mathit{\partial }G(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.1)

In Problem 2:(5.2) $$\begin{array}{c}\hfill \overline{w}(\mathit{\zeta })={\int }_l{w}_1(y)\frac{\mathit{\partial }{G}_1^{(2)}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.2)

In Problem 3:(5.3) $$\begin{array}{c}\hfill {\mathit{\mu }}_1\overline{w}(\mathit{\zeta })={\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y+{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.3)

Here the terms f in Equations (2.6(2.6) $$\begin{array}{cc}\hfill {\displaystyle w(\mathit{\xi })}& =f(\mathit{\xi })+{\int }_{l}w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y,\hfill \\ \hfill {\displaystyle f(\mathit{\xi })}& =\frac{p_0}{\mathit{\mu }}{\left.G^{(1)}(y,\mathit{\xi })\right|}_{y_1=x_0,y_2=h}=\frac{p_0}{\mathit{\mu }}{G}^{(1)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),\hfill \end{array}$$(2.6) ) and (2.14(2.14) $$\begin{array}{c}\hfill {\displaystyle w_1(\mathit{\xi })=f(\mathit{\xi })+{\int }_l{w}_1(y)\frac{\mathit{\partial }{G}_1^{(2)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}l,f(\mathit{\xi })=\frac{p_0}{{\mathit{\mu }}_1}{G}_1^{(2)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),}\end{array}$$(2.14) ), as well as terms $f_{11},f_{21}$ in Equation (2.20(2.20a) & (2.20b) $$\begin{array}{cc}& {\displaystyle x\in l_1:\frac{{\mathit{\mu }}_1}{2}{w}_1(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{11}(x),f_{11}(\mathit{\xi })=p_0{G}_{11}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),}\hfill \\ \hfill & {\displaystyle x\in l_2:\frac{{\mathit{\mu }}_2}{2}{w}_2(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{22}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{21}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{21}(x),f_{21}(\mathit{\xi })=p_0{G}_{21}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2).}\hfill \end{array}$$(2.20a) & (2.20b) ), all have been omitted. This means that the amplitude of the free-surface oscillations, function $\overline{w}(\mathit{\zeta })$ is taken in its relative form, when compared with the analogous amplitude for the medium without any buried object.

The inverse problem on object’s image identification is to find contour l from Equations (5.1(5.1) $$\begin{array}{c}\hfill \overline{w}(\mathit{\zeta })={\int }_{l}w(y)\frac{\mathit{\partial }G(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.1) )–(5.3(5.3) $$\begin{array}{c}\hfill {\mathit{\mu }}_1\overline{w}(\mathit{\zeta })={\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y+{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.3) ), if the input data i.e. function $\overline{w}(\mathit{\zeta })$ is known, at least with a certain random noise, over a chosen interval $\mathit{\zeta }\in [c,d]$. It should be noted that not only contour l is unknown in Equations (5.1(5.1) $$\begin{array}{c}\hfill \overline{w}(\mathit{\zeta })={\int }_{l}w(y)\frac{\mathit{\partial }G(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.1) )–(5.3(5.3) $$\begin{array}{c}\hfill {\mathit{\mu }}_1\overline{w}(\mathit{\zeta })={\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y+{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.3) ) but also functions $w(y),w_1(y),[w_1(y),w_2(y)]$ are unknown on the unknown contour l. The so-formulated mathematically inverse problem is obviously non-linear with respect to a constitutive equation of this contour l. Generally, this problem is nonlinear and ill-posed. However, as indicated above, in the case of elliptic object this problem is solved in a finite-dimensional space, hence this becomes well-posed. The questions of existence and uniqueness are out of the scope of the present study.

Let us rewrite the basic BIEs (2.7(2.7) $$\begin{array}{cc}\hfill {\displaystyle \frac{w(x)}{2}-{\int }_{l}K(y,x)w(y)\mathrm{d}{l}_y}& =f(x),x\in l,\hfill \\ \hfill {\displaystyle K(y,x)}& =\frac{\mathit{\partial }{G}^{(1)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}^{(1)}·n_y\}=n_y^1\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}^{(1)}}{\mathit{\partial }{y}_2}.\hfill \end{array}$$(2.7) ), (2.15(2.15) $$\begin{array}{cc}& {\displaystyle \frac{w_1(x)}{2}-{\int }_{l}K(y,x)w_1(y)\mathrm{d}{l}_y=f(x),x\in l,}\hfill \\ \hfill & {\displaystyle K(y,x)=\frac{\mathit{\partial }{G}_1^{(2)}(y,x)}{\mathit{\partial }{n}_y}=\{\mathrm{\nabla }{G}_1^{(2)}.n_y\}=n_y^1\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_1}+n_y^2\frac{\mathit{\partial }{G}_1^{(2)}}{\mathit{\partial }{y}_2}.}\hfill \end{array}$$(2.15) ), (2.20(2.20a) & (2.20b) $$\begin{array}{cc}& {\displaystyle x\in l_1:\frac{{\mathit{\mu }}_1}{2}{w}_1(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{11}(x),f_{11}(\mathit{\xi })=p_0{G}_{11}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),}\hfill \\ \hfill & {\displaystyle x\in l_2:\frac{{\mathit{\mu }}_2}{2}{w}_2(x)-{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{22}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y}\hfill \\ \hfill & {\displaystyle -{\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{21}(y,x)}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y=f_{21}(x),f_{21}(\mathit{\xi })=p_0{G}_{21}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2).}\hfill \end{array}$$(2.20a) & (2.20b) ), for any direct Problem 1 – 3, symbolically as an operator equation:(5.4) $$\begin{array}{c}\hfill \mathbf{A}g=f,\end{array}$$(5.4)

where in Problem 1: ${g=w|}_l$; in Problem 2: $g=w_1{|}_l$; in Problem 3: $g=g_1\cup g_2=(w_1{|}_{l_1})\cup (w_2{|}_{l_2})$. The numerical collocation technique, developed in the previous section, permits efficient solution of respective BIEs written in a certain discrete form; this implies the numerical inversion of operator $\mathbf{A}$, if the boundary contour l of the void is known:(5.5) $$\begin{array}{c}\hfill g={\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f.\end{array}$$(5.5)

Here we have indicated explicitly that the inverse operator ${\mathbf{A}}^{-1}$ depends upon the geometry of the elliptic contour l, which is defined at most by the five geometric parameters: two Cartesian coordinates of its center $(x_1^c,x_2^c)$, the lengths of its two semi-axes a and b, and the slope angle $\mathit{\theta }$ between the principal axis and the horizontal coordinate axis $x_1$.

Further, let us rewrite symbolically relations (5.1(5.1) $$\begin{array}{c}\hfill \overline{w}(\mathit{\zeta })={\int }_{l}w(y)\frac{\mathit{\partial }G(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.1) )–(5.3(5.3) $$\begin{array}{c}\hfill {\mathit{\mu }}_1\overline{w}(\mathit{\zeta })={\mathit{\mu }}_1{\int }_{l_1}{w}_1(y)\frac{\mathit{\partial }{G}_{11}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y+{\mathit{\mu }}_2{\int }_{l_2}{w}_2(y)\frac{\mathit{\partial }{G}_{12}(y,\mathit{\zeta })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y.\end{array}$$(5.3) ) in a unified form, as a certain operator relation having the input data, function $\overline{w}(\mathit{\zeta })$, in the right-hand side:(5.6) $$\begin{array}{c}\hfill [\mathbf{K}(x_1^c,x_2^c,a,b,\mathit{\theta })g](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d],\end{array}$$(5.6)

where again operator $\mathbf{K}$ depends on the geometry of contour l which is equivalent to its dependence again upon the five parameters $(x_1^c,x_2^c,a,b,\mathit{\theta })$. By substitution of (5.5(5.5) $$\begin{array}{c}\hfill g={\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f.\end{array}$$(5.5) ) into (5.6(5.6) $$\begin{array}{c}\hfill [\mathbf{K}(x_1^c,x_2^c,a,b,\mathit{\theta })g](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d],\end{array}$$(5.6) ), one obtains:(5.7) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d].\end{array}$$(5.7)

In Equation (5.7(5.7) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d].\end{array}$$(5.7) ) all included functions and operators are known. The only unknown quantities are five geometric parameters defining the geometry of the object’s boundary curve l. If the amplitude of vibrations of the upper surface ($y_2=h$) is known exactly as a solution to the direct problem, then the solution of the operator Equation (5.7(5.7) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d].\end{array}$$(5.7) ) should coincide with the real contour l which is used for the direct problem. However, the construction of a stable numerical algorithm, to determine five parameters $(x_1^c,x_2^c,a,b,\mathit{\theta })$ from operator Equation (5.7(5.7) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d].\end{array}$$(5.7) ) is not simple, is one does not known a priori the contour l of the direct problem. This question is much more complex in the case when the input data, function $\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d]$ is extracted from the solution of the respective direct problem, further perturbed by a random noise. Therefore, the proposed algorithm for object identification should be stable with respect to small perturbations of the input data.

6. Numerical algorithm and examples of the reconstruction of elliptical voids

Let us assume that function $\overline{w}(\mathit{\zeta })$ in Equation (5.7(5.7) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d].\end{array}$$(5.7) ) is known in a discrete set of nodes ${\mathit{\zeta }}_q\in [c,d],(q=1,\cdots ,Q)$. Then Equation (5.7(5.7) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f](\mathit{\zeta })=\overline{w}(\mathit{\zeta }),\mathit{\zeta }\in [c,d].\end{array}$$(5.7) ) in the discrete form reads as follows:(6.1) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f]({\mathit{\zeta }}_q)=\overline{w}({\mathit{\zeta }}_q),q=1,\cdots ,Q.\end{array}$$(6.1)

The proposed numerical algorithm is based on a minimization of the discrepancy functional $\mathrm{\Omega }$ [19] of Equation (6.1(6.1) $$\begin{array}{c}\hfill [\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f]({\mathit{\zeta }}_q)=\overline{w}({\mathit{\zeta }}_q),q=1,\cdots ,Q.\end{array}$$(6.1) ):(6.2) $$\begin{array}{cc}\hfill {\displaystyle \mathrm{\Omega }(x_1^c,x_2^c,a,b,\mathit{\theta })}& ={\left|\left|[\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f]({\mathit{\zeta }}_q)-\overline{w}({\mathit{\zeta }}_q)\right|\right|}^2\hfill \\ \hfill {\displaystyle }& =\sum _{q=1}^Q{\left\{[\mathbf{K}{\mathbf{A}}^{-1}(x_1^c,x_2^c,a,b,\mathit{\theta })f]({\mathit{\zeta }}_q)-\overline{w}({\mathit{\zeta }}_q)\right\}}^2,\hfill \end{array}$$(6.2)

which in the non-linear problem is generally non-quadratic.

It should be noted that, in spite of the formulated inverse problems are theoretically finite-dimensional, hence well-posed, in the numerical implementation this may lead to a practical instability in the case of the noisy input data. By this reason, we used in our numerical experiments the steepest descent iteration method which is known to provide a sort of iterative regularization [20]. The method applied thus brings an additional stability to the identification algorithm proposed.

Some examples of the reconstruction are demonstrated here in tables. For all examples, we used $Q=20$, to form the array of the input data, such that ${\mathit{\zeta }}_q=c+m{h}^{\ast },h^{\ast }=(d-c)/Q$ , $q=0,1,...,Q$. It should be noted that for all tables we choose the following criterion, to stop the iteration process on the nth step: ${\mathit{\epsilon }}_n=\left|[\mathrm{\Omega }(x^{(n+1)})-\mathrm{\Omega }(x^{(n)})]/\mathrm{\Omega }(x^{(n)})\right|\sim {10}^{-4}$ .

Table 1. Examples of the reconstruction, Problem 1: $k=1.5$, $h=2.0$, $x_0=6.0$, ${(x_1^c,x_2^c,a,b,\mathit{\theta })}^0=(2/3,2/3,0.25,0.25,\mathit{\pi }/3)$, Noise level = 0 %, bold font – exact geometry.

${\mathit{\epsilon }}_n$	No. iter.	$x_1^c$	$x_2^c$	a	b	$\mathit{\theta }$	$[c,d]$
		1.25	1.2	0.2	0.3	$\mathit{\pi }\mathbf{/}\mathbf{4}\mathbf{=}\mathbf{0.785}$
$5.2·{10}^{-4}$	10	1.354	1.138	0.336	0.236	1.068	[$-1,1$]
$2.1·{10}^{-4}$	22	1.289	1.233	0.228	0.294	0.918	[$-2,2$]
$1.7·{10}^{-4}$	33	1.274	1.185	0.224	0.276	0.802	[$-3,3$]
$1.3·{10}^{-4}$	12	1.260	1.190	0.203	0.297	0.775	[$-4,4$]
		1.0	1.0	0.25	0.2	$\mathit{\pi }\mathbf{/}\mathbf{6}\mathbf{=}\mathbf{0.524}$
$3.0·{10}^{-4}$	23	1.017	1.100	0.320	0.195	0.631	[$-1,1$]
$1.1·{10}^{-4}$	7	1.042	1.0360	0.202	0.298	0.652	[$-2,2$]
$1.5·{10}^{-4}$	17	1.008	1.014	0.196	0.212	0.695	[$-3,3$]
$1.0·{10}^{-4}$	25	0.983	1.055	0.236	0.180	0.540	[$-4,4$]

Table 1 reflects some identification results for Problem 1. This is principally aimed at the study of identification precision versus the length of the sensing interval $[c,d]$.

Table 2 demonstrates object identification for Problem 2. Here only horizontal principal axis of the elliptic object has been used in our numerical experiments. The number of the geometrical parameters to be reconstructed here is four since in this case the slope angle $\mathit{\theta }$ is excluded from the array of standard five parameters.

Table 2. Examples of the reconstruction, Problem 2: $k_1=0.62$, $k_2=0.22$, $\mathit{\theta }=0$, $x_0=4.25$, $[c,d]=[-3,3]$, bold font – exact geometry.

${(x_1^c,x_2^c,a,b)}^0$	Noise (%)	${\mathit{\epsilon }}_n$	No. iter.	$x_1^c$	$x_2^c$	a	b
$h_1=h_2=1.0$				1.35	1.2	0.15	0.1
(1.5, 1.5, 0.25, 0.25)	0	$0.91·{10}^{-4}$	12	1.445	1.174	0.140	0.095
	10	$0.42·{10}^{-4}$	15	1.444	1.163	0.138	0.095
	20	$0.21·{10}^{-4}$	23	1.443	1.170	0.141	0.086
$h_1=h_2=1.0$				1.75	1.6	0.1	0.3
(1.5, 1.5, 0.25, 0.25)	0	$0.01·{10}^{-4}$	17	1.728	1.583	0.111	0.274
	10	$0.32·{10}^{-4}$	25	1.694	1.542	0.120	0.254
	20	$0.16·{10}^{-4}$	27	1.722	1.571	0.115	0.265
$h_1=2.0,h_2=1.0$				1.5	1.5	0.25	0.25
(2.0, 2.0, 0.5, 0.5)	0	$0.74·{10}^{-4}$	26	1.403	1.602	0.278	0.232
	10	$0.40·{10}^{-4}$	32	1.399	1.663	0.315	0.205
	20	$0.23·{10}^{-4}$	32	1.402	1.634	0.250	0.227
$h_1=2.0,h_2=1.0$				2.25	1.75	0.4	0.1
(2.0, 2.0, 0.5, 0.5)	0	$0.11·{10}^{-4}$	11	2.217	1.718	0.385	0.118
	10	$0.14·{10}^{-4}$	11	2.181	1.680	0.361	0.120
	20	$0.35·{10}^{-4}$	12	2.202	1.653	0.374	0.120

Table 3. Examples of the reconstruction, Problem 3: $k_1=0.62$, $k_2=0.22$, $h_1=h_2=0.75$, $x_0=3.25$, ${(x_1^c,a,b)}^0=(1.0,1.0,2.0)$, $[c,d]=[-3,3]$, bold font – exact geometry.

Noise (%)	No. iter.	$x_1^c$	a	b
		0.75	0.25	0.25
0	8	0.750	0.250	0.250
10	9	0.748	0.250	0.245
20	8	0.739	0.267	0.235
		1.0	0.1	0.35
0	10	0.984	0.116	0.362
10	6	0.976	0.124	0.382
20	12	0.961	0.135	0.380

Table 3 demonstrates the precision of the object’s image reconstruction for Problem 3. Only three geometrical parameters are to be restored from the list of five parameters in the general case. In fact, the centre of the void is placed just on the interface straight line $y=h_2$, hence $\mathit{\theta }=0$ and $x_2^c=h_2$ in this case.

Table 4. Examples of the reconstruction with different positions of void’s central point, Problem 1, a rectangular void approximated by an elliptic defect: $k=0.38$, $h=4.0$, $a=0.375h$, $b=0.25h$, $x_0=6.0$, $\mathit{\theta }=0$, $n=4$, $[c,d]=[-5,5]$,  ${(x_1^c,x_2^c,a,b)}^0=(1.25,1.25,1.25,1.25)$, Noise level = 10 %, bold font – exact configuration.

$x_1^c$	$x_2^c$	Relative error of void’s area (%)
1.0	1.0
1.035	0.962	1.4
1.15	1.0
1.282	0.961	1.7
1.1	0.9
1.161	0.862	3.2
0.85	1.2
0.885	1.261	3.9

Moreover, we have performed a numerical experiment on reconstruction of a non-elliptic object, more exactly – a rectangular void with filleted corners, by approximating it by an ellipse, see Table 4. For this aim, we first have solved the direct problem for the exact object (rectangle) as shown in Figure 8. The so obtained input data, the calculated displacement over the upper surface of the layer, is used to solve the inverse problem in the class of elliptic defects by the algorithm presented above. As can be seen from Table 4, the relative error in the identification of defect’s area does not exceed 4%. The maximum relative error in the identification of the central point coordinates is around 6%.

At last, let us compare our results with some others known in literature. The authors could not find any published paper, for a direct comparison. However, Pompei et al. [21] solve a static anti-plane problem for identification of elliptic voids in an elastic half-space. Their results can be obtained from the present work, Problem 1, as a particular case when thickness h of our layer tends to infinity: $h\to \infty$, and wave number k becomes very small: $k\to 0$. Besides, we should reverse the direction of the vertical axis $y_2$ to be directed downward. This requires to modify slightly the boundary conditions, namely, the harmonic load to act on the upper boundary $y_2=0$, and the lower boundary $y_2=h$ to be fixed. With so doing, the representation for Green’s function may directly be taken from [14], chapter 3. The only change is that function $f(\mathit{\xi })$ in Equation (2.6(2.6) $$\begin{array}{cc}\hfill {\displaystyle w(\mathit{\xi })}& =f(\mathit{\xi })+{\int }_{l}w(y)\frac{\mathit{\partial }{G}^{(1)}(y,\mathit{\xi })}{\mathit{\partial }{n}_y}\mathrm{d}{l}_y,\hfill \\ \hfill {\displaystyle f(\mathit{\xi })}& =\frac{p_0}{\mathit{\mu }}{\left.G^{(1)}(y,\mathit{\xi })\right|}_{y_1=x_0,y_2=h}=\frac{p_0}{\mathit{\mu }}{G}^{(1)}(x_0,h,{\mathit{\xi }}_1,{\mathit{\xi }}_2),\hfill \end{array}$$(2.6) ) should be taken with opposite sign.

The comparison of our results with [21], on a particular example is demonstrated in Table 5.

Table 5. Comparison with [21], Problem 1: $k=1·{10}^{-4}$, $x_0=-6.0$, ${(x_1^c,x_2^c,a,b,\mathit{\theta })}^0=(1.0,1.0,0.5,0.5,\mathit{\pi }/6)$, Noise level = 0 %, $[c,d]=[-5,5]$.

${\mathit{\epsilon }}_n$	No. iter.	$x_1^c$	$x_2^c$	a	b	$\mathit{\theta }$
Exact		1.5	0.0	1.0	0.5	0.0
[21]		1.520	0.025	0.999	0.527	0.020
$h=20$	63	1.545	0.018	1.048	0.709	0.213
$h=40$	60	1.494	0.005	1.039	0.604	0.015
Exact		2.5	1.0	1.3	0.5	$\mathit{\pi }\mathbf{/}\mathbf{2}\mathbf{=}\mathbf{1.571}$
[21]		2.470	1.001	0.493	1.293	3.140
$h=20$	242	2.533	1.071	0.491	1.278	0.003
$h=40$	206	2.490	1.013	1.340	0.492	1.515

Note: The bold values in Table 5 are significant, to disctinct exact values from the reconstructed ones.

7. Conclusions

The basic results and the principal conclusions are listed here below:

(1)	The paper proposes some efficient methods to solve the anti-plane dynamic problem of the linear isotropic elasticity theory, for layered elastic media with a buried object located inside the media.
(2)	On the basis of the constructed Greens’ functions, the direct problem is reduced to a single BIE or to a system of BIEs, holding over the boundary contour of the buried object.
(3)	There is developed an efficient numerical algorithm for the system of integral equations on the basic of collocation method, which solves the formulated direct diffraction problems.
(4)	There is developed a numerical algorithm, to solve the inverse object identification problems in the layered media, which is founded on the minimization of the discrepancy functional. Some examples demonstrating the solution to the direct and the inverse problems are presented as diagrams and tables. The solution of the inverse problems is adapted to the case of elliptical voids.
(5)	The precision of the reconstruction weakly depends on the noise level. This is much more sensitive to the geometry of the object, as well as to the values of other geometric parameters, such as the length of the interval $[c,d]$.
(6)	It is clearly seen from Table 1 that object’s identification becomes more precise if the interval of the ‘measurement’ $[c,d]$, to create the input data base, becomes wider.
(7)	Table 2 confirms that the elongated objects are identified with a less precision that is quite natural from the heuristic point of view.
(8)	Table 3 demonstrates that the identification with less number of parameters (only three geometrical parameters) is performed with a better precision than in the general case of five parameters.
(9)	Table 4 shows a good precision in the identification of the defects of non-elliptic shapes. It is very important that the area of the defect is identified precisely. This fact, of course, is of a great practical significance in the ultrasonic NDT.
(10)	Table 5 demonstrates an acceptable precision when our dynamic results of the identification in the layer are compared with a static problem in the half-plane. This comparison is admissible for low frequencies and deep layer. It is also interesting that some results of the reconstruction in Table 5, which initially seem to be far from the true geometry, in fact describe a void very close to the exact geometry. The difference becomes inessential there since the reconstructed object is turned through right angle, and simultaneously the bigger axis of the ellipse becomes its smaller one.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The second author is grateful to the Russian Science Foundation (RSCF), for its support by Project [number 15-19-10008].