Weighted integral solvers for elastic scattering by open arcs in two dimensions

TL;DR

This paper introduces a new numerical method for elastic scattering problems involving open arcs in two dimensions, combining weighted operators, a novel Calderón relation, and spectral quadrature for high accuracy and efficiency.

Contribution

It develops a novel approach using weighted operators and a new Calderón relation for elastic scattering by open arcs, with demonstrated high accuracy and efficiency.

Findings

High-accuracy results with few iterations

Effective for both low and high frequencies

Numerical examples confirm efficiency

Abstract

We present a novel approach for the numerical solution of problems of elastic scattering by open arcs in two dimensions. Our methodology relies on the composition of weighted versions of the classical operators associated with Dirichlet and Neumann boundary conditions in conjunction with a certain "open-arc elastic Calder\'on relation" whose validity is demonstrated in this paper on the basis of numerical experiments, but whose rigorous mathematical proof is left for future work. Using this Calder\'on relation in conjunction with spectrally accurate quadrature rules and the Krylov-subspace linear algebra solver GMRES, the proposed overall open-arc elastic solver produces results of high accuracy in small number of iterations---for low and high frequencies alike. A variety of numerical examples in this paper demonstrate the accuracy and efficiency of the proposed methodology.

Tables

Table 1. Table 1: Near-field errors for elastic scattering by a unit circle. GMRES tol: 10 − 12 superscript 10 12 10^{-12} .

$\omega$	$N$	$\text{Dir}\left(S\right)$	$\text{Dir}\left(NS\right)$	$\text{Dir}\left(\tilde{N}S\right)$	$\text{Neu}\left(N\right)$	$\text{Neu}\left(NS\right)$
	30	$6.59\times {10}^{-4}$	$6.59\times {10}^{-4}$	$6.59\times {10}^{-4}$	$9.53\times {10}^{-3}$	$9.53\times {10}^{-3}$
10	40	$1.55\times {10}^{-6}$	$1.55\times {10}^{-6}$	$1.55\times {10}^{-6}$	$6.02\times {10}^{-5}$	$6.02\times {10}^{-5}$
	60	$2.54\times {10}^{-12}$	$2.50\times {10}^{-12}$	$2.68\times {10}^{-12}$	$2.13\times {10}^{-11}$	$2.16\times {10}^{-11}$
	150	$1.33\times {10}^{-7}$	$1.33\times {10}^{-7}$	$1.33\times {10}^{-7}$	$1.58\times {10}^{-7}$	$1.58\times {10}^{-7}$
50	160	$2.51\times {10}^{-10}$	$2.51\times {10}^{-10}$	$2.54\times {10}^{-10}$	$1.27\times {10}^{-9}$	$1.27\times {10}^{-9}$
	200	$7.12\times {10}^{-13}$	$1.58\times {10}^{-12}$	$9.53\times {10}^{-12}$	$2.22\times {10}^{-12}$	$1.65\times {10}^{-12}$

Table 2. Table 2: Near-field errors for elastic scattering by a Spiral-Shaped Arc. GMRES tol: 10 − 8 superscript 10 8 10^{-8} .

$\omega$	$N$	$\text{Dir}\left(S^w\right)$	$\text{Dir}\left(N^w{S}^w\right)$	$\text{Dir}\left({\tilde{N}}^w{S}^w\right)$	$\text{Neu}\left(N^w\right)$	$\text{Neu}\left(N^w{S}^w\right)$
	100	$2.97\times {10}^{-2}$	$2.97\times {10}^{-2}$	$2.97\times {10}^{-2}$	$2.22\times {10}^{-1}$	$2.22\times {10}^{-1}$
10	150	$7.41\times {10}^{-5}$	$7.41\times {10}^{-5}$	$7.41\times {10}^{-5}$	$4.33\times {10}^{-3}$	$4.33\times {10}^{-3}$
	200	$2.06\times {10}^{-11}$	$2.39\times {10}^{-11}$	$2.10\times {10}^{-11}$	$7.42\times {10}^{-10}$	$7.43\times {10}^{-10}$
	600	$6.11\times {10}^{-3}$	$1.95\times {10}^{-4}$	$6.11\times {10}^{-3}$	$1.49\times {10}^{-1}$	$4.22\times {10}^{-2}$
50	800	$4.76\times {10}^{-10}$	$8.51\times {10}^{-8}$	$1.19\times {10}^{-7}$	$1.46\times {10}^{-6}$	$1.86\times {10}^{-6}$
	1000	$2.01\times {10}^{-14}$	$1.33\times {10}^{-7}$	$1.97\times {10}^{-12}$	$1.57\times {10}^{-8}$	$9.36\times {10}^{-9}$

Table 3. Table 3: Iterations for elastic scattering by a Flat Strip. GMRES tol: 10 − 5 superscript 10 5 10^{-5} .

$\omega$	$N$	$\text{Dir}\left(S^w\right)$	$\text{Dir}\left(N^w{S}^w\right)$	$\text{Dir}\left({\tilde{N}}^w{S}^w\right)$	$\text{Neu}\left(N^w\right)$	$\text{Neu}\left(N^w{S}^w\right)$
10	160	21	12	12	31	13
30	480	33	23	16	88	24
50	800	41	30	19	137	34
80	1280	48	43	21	204	47
100	1600	51	50	23	246	58

Table 4. Table 4: Computing time (seconds) required for elastic scattering by a Flat Strip. GMRES tol: 10 − 5 superscript 10 5 10^{-5} .

$\omega$	$N$	$\text{Dir}\left(S^w\right)$	$\text{Dir}\left(N^w{S}^w\right)$	$\text{Dir}\left({\tilde{N}}^w{S}^w\right)$	$\text{Neu}\left(N^w\right)$	$\text{Neu}\left(N^w{S}^w\right)$
10	160	0.01	0.06	0.06	0.17	0.06
30	480	0.07	0.41	0.29	2.9	0.53
50	800	0.22	1.3	0.73	11.2	1.9
80	1280	0.60	4.4	2.1	40.4	8.9
100	1600	0.97	7.7	3.3	76.3	15.5

Table 5. Table 5: Iterations for elastic scattering by a Spiral-Shaped Arc. GMRES tol: 10 − 5 superscript 10 5 10^{-5} .

$\omega$	$N$	$\text{Dir}\left(S^w\right)$	$\text{Dir}\left(N^w{S}^w\right)$	$\text{Dir}\left({\tilde{N}}^w{S}^w\right)$	$\text{Neu}\left(N^w\right)$	$\text{Neu}\left(N^w{S}^w\right)$
10	160	73	75	61	198	76
30	480	109	139	94	444	145
50	800	126	196	113	651	204
80	1280	147	274	139	985	291
100	1600	160	330	157	1222	348

Table 6. Table 6: Computing time (seconds) for elastic scattering by a Spiral-Shaped Arc. GMRES tol: 10 − 5 superscript 10 5 10^{-5} .

$\omega$	$N$	$\text{Dir}\left(S^w\right)$	$\text{Dir}\left(N^w{S}^w\right)$	$\text{Dir}\left({\tilde{N}}^w{S}^w\right)$	$\text{Neu}\left(N^w\right)$	$\text{Neu}\left(N^w{S}^w\right)$
10	160	0.08	0.25	0.18	0.8	0.23
30	480	0.36	2.7	1.6	10.7	3.7
50	800	0.93	15.9	4.6	41.9	15.7
80	1280	2.5	47.6	13.7	183.2	51.2
100	1600	3.7	53.9	23.7	347.8	93.8

Equations

$${\Delta }^{\ast }u+\rho {\omega }^{2}u=0\text{mbox}in{\mathbb{R}}^2\backslash \Gamma ,$$

$$u=F\text{mbox}on\Gamma$$

$$T(\partial ,\nu )u=G\text{mbox}on\Gamma .$$

$${\Delta }^{\ast }=\mu \text{mbox}div\text{mbox}grad+(\lambda +\mu )\text{mbox}grad\text{mbox}div,$$

$$Tu:=2\mu {\partial }_{\nu }u+\lambda \nu \mathrm{div}u+\mu \nu \times \mathrm{curl}u,\nu =({\nu }^1,{\nu }^2){}^{\top },$$

$$\tilde{T}u:=(\mu +\tilde{\mu }){\partial }_{\nu }u+\tilde{\lambda }\nu \mathrm{div}u+\tilde{\mu }\nu \times \mathrm{curl}u,$$

$$\tilde{\mu }=\frac{\mu (\lambda +\mu )}{\lambda +3\mu },\tilde{\lambda }=\lambda +\mu -\tilde{\mu },$$

$$u(x)=(\mathcal{S}\Psi )(x):={\int }_{\Gamma }E(x,y))\varphi (y)d{s}_y,\forall x\in {\mathbb{R}}^2\backslash \Gamma ,$$

$$u(x)=(\mathcal{D}\Psi )(x):={\int }_{\Gamma }(T({\partial }_y,{\nu }_y)E(x,y){)}^{\top }\psi (y)d{s}_y,\forall x\in {\mathbb{R}}^2\backslash \Gamma ,$$

$${\gamma }_{k_t}(x,y)=\frac{i}{4}{H}_0^{(1)}(k_t|x-y|),x\ne y,$$

$$k_s:=\omega /c_p,k_p=\omega /c_s$$

$$c_p=\sqrt{\mu /\rho }\text{mbox}and{c}_s=\sqrt{(\lambda +2\mu )/\rho },$$

$$E(x,y)=\frac{1}{\mu }{\gamma }_{k_s}(x,y)I+\frac{1}{\rho {\omega }^2}{\nabla }_x{\nabla }_x^{\top }\left[{\gamma }_{k_s}(x,y)-{\gamma }_{k_p}(x,y)\right].$$

$$S[\varphi ](x)$$

$$N[\psi ](x)$$

$$S[\varphi ]=F,N[\psi ]=G\text{mbox}on\Gamma .$$

$$\varphi \sim \frac{{\xi }_1}{\sqrt{d}}+{\eta }_1,\phi \sim {\xi }_2\sqrt{d}+{\eta }_2,$$

$$\varphi =\frac{\alpha }{\sqrt{d}},\phi =\beta \sqrt{d},$$

$$S^w[\alpha ]=S\left[\frac{\alpha }{w}\right],N^w[\beta ]=N[\beta w],$$

$$S^w[\alpha ]=F,N^w[\beta ]=G\text{mbox}on\Gamma ,$$

$$\displaystyle{\color[rgb]{0.000,0.000,0.000}N_{a,c}S_{a,c}=-\frac{I}{4}+(K_{a,c}^{*})^{2}},$$

$$\displaystyle{\color[rgb]{0.000,0.000,0.000}K^{*}}[\psi](x)$$

$$\displaystyle{\color[rgb]{0.000,0.000,0.000}\widetilde{K}^{*}}[\psi](x)$$

$$\tilde{N}[\psi ](x)$$

$$C_{\lambda ,\mu }=\frac{\mu }{2(\lambda +2\mu )}.$$

$$(K^{\ast }{)}^2-C_{\lambda ,\mu }^{2}I=(K^{\ast }-K_0^{\ast })(K^{\ast }+K_0^{\ast })+(K_0^{\ast }{)}^2-C_{\lambda ,\mu }^{2}I$$

$$NS=-\left(\frac{1}{4}-C_{\lambda ,\mu }^2\right)I+K^{(1)},\tilde{N}S=-\frac{1}{4}I+K^{(2)},$$

$$N^w{S}^w[\alpha ]=N^w[F],{\tilde{N}}^w{S}^w[\alpha ]={\tilde{N}}^w[F],N^w{S}^w[\beta ]=G\text{mbox}on\Gamma .$$

$$\tilde{T}(\partial ,\nu )u(x)=(\lambda +\mu )\nu (\nabla \cdot u)+\mu {\partial }_{\nu }u+\tilde{\mu }M(\partial ,\nu )u$$

$$M(\partial ,\nu )u(x)={\partial }_{\nu }u-\nu (\nabla \cdot u)+\nu \times \mathrm{curl}u.$$

$$M(\partial ,\nu )u(x)=A\frac{du(x)}{d{s}_x},A=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right).$$

Full text

1 Introduction

We consider the problem of numerical evaluation of elastic waves diffracted by infinitely thin open surfaces [24, 27, 28, 29, 31]. This problem plays fundamental roles in a number of important applications in science and engineering, including non-destructive testing of materials, characterization of fractures, energy production from natural gas and geothermal resources, mining, etc. These problems present considerable mathematical and computational challenges in view of the highly oscillatory character of the associated time-harmonic elastic fields and, for solvers based on use of volumetric discretizations, the unbounded character of the physical domains that must be considered in connection with the aforementioned applications.

As is the case for wave scattering problems in acoustics and electromagnetics [15, 25], the boundary integral equation methods in elasticity require discretization of domains of lower dimensionality than those required by volumetric discretization methods (such as finite difference or finite element methods [7, 29]). These equations can generally be treated effectively even for high frequencies [12, 8, 13, 9]. As is well known, however, the classical boundary integral equations for open surfaces (or, in two dimensions, open arcs) are not second-kind Fredholm integral equations. This setting presents some difficulties. On one hand, as the eigenvalues of the left-hand side operators accumulate around zero and/or infinity, solution of these problems by means of Krylov-subspace iterative solvers such as GMRES (which are commonly-used in conjucntion with accelerated integral-equation solvers) often requires large numbers of iterations for convergence, and thus, large computing costs. Additionally, the evaluation of the Hadamard finite part of the hypersingular integral operator associated to the elastic Neumann problem [2] has also remained a significant challenges in this context [13].

These difficulties were recently addressed, for the 2D and 3D acoustic and 2D electromagnetic contexts, in the contributions [9, 10, 22]—which, in particular, introduced Fredholm integral equations of the second-kind, and associated numerical algorithms, for problems of acoustic scattering by open arcs and surfaces. The second-kind equations were obtained in these articles by utilizing compositions of appropriately modified versions {\color[rgb]{0.000,0.000,0.000}S_{a}^{w}} and {\color[rgb]{0.000,0.000,0.000}N_{a}^{w}} (which incorporate explicitly the singular character of the integral equation densities) of the acoustic single-layer and hypersingular integral operators {\color[rgb]{0.000,0.000,0.000}S_{a}} and {\color[rgb]{0.000,0.000,0.000}N_{a}}. As shown in [9, 22], a generalization of the classical closed-surface Calderón formulas holds for the {\color[rgb]{0.000,0.000,0.000}S_{a}^{w}} and {\color[rgb]{0.000,0.000,0.000}N_{a}^{w}} operators in the open arc case, which, in particular, gives rise to the aforementioned second-kind integral formulations for the acoustic open-arc problems. It was verified numerically in these contributions that, as predicted by theory, the eigenvalues of the proposed second-kind operators remain bounded away from zero and infinity, even for problems of high frequency.

Multiple challenges arise as extensions of these methods to elastic open-arc problems are attempted. At a basic level, elastic Calderón formulas have not been studied even in the closed-surface case—likely, on account of the fact that, unlike the acoustic wave case, the classical Neumann-Poincaré double-layer operator $K$ and its adjoint $K^{*}$, which play important roles in the Calderón relations, are not compact in the elastic case [4, 21]. (Reference [17] mentions all the elastic operators relevant to the closed-surface Calderón calculus but it does not utilize the Calderón projectors as regularization tools). But recent results [4] have established that the closed surface elastostatic ($\omega=0$) double-layer integral operator is polynomially compact, which suggests that the composition $NS$ of closed-surface single-layer and hypersingular integral operators may be a Fredholm operator of second-kind even in the case $\omega\neq 0$. This in fact established in Section 3.2 below, where it is additionally shown that the eigenvalues of the operator $NS$ accumulate at a certain point which depends on the elastic Lamé parameters. Moreover, for the Dirichlet problem particulary, an artificial traction operator $\widetilde{T}$ can be introduced (see (2.3) for which the resulting double-layer operator $\widetilde{K}$ is compact [18, 21]—which results in an alternative second-kind Calderón-like formula for the composition $\widetilde{N}S$. But this formula is only applicable for Dirichlet problem since this artificial traction operator does not correspond to an actual physical traction.

Relying on the newly studied property of the closed-surface Calderón formula for elastic wave, it is naturally to consider the extension of the existing second-kind Fredholm integral equations for acoustic open-arcs [9, 22] to the elastic case. In view of these contributions, and unlike the approaches [3, 14, 19, 23, 26, 32], we consider the composition of weighted versions of the classical single-layer and hypersingular integral operators [2]. We find that the benefits of this approach are two-fold. On one hand, the new method enjoys high-order accuracy: the weighted versions $S^{w}$, $N^{w}$ of the single-layer and hypersingular operators extract the solutions’ edge singularity [16] explicitly and the applied quadrature rules provide spectral convergence. On the other hand, the method gives rise to well-behaved linear algebra: as numerically demonstrated in Section 3.3, the eigenvalues of $N^{w}S^{w}$ and $\widetilde{N}^{w}S^{w}$ are bounded away from zero and infinity. (A theoretical proof of this fact is left for future work.). And, as desired, the composite operator $N^{w}S^{w}$ or $\widetilde{N}^{w}S^{w}$ requires small number of iterations when used in conjunction with the linear iterative solver GMRES. The new formulation for the Neumann problem is especially beneficial, as it give rise to order-of-magnitude improvements in computing times over the corresponding hypersingular formulation. Such gains do not occur for the Dirichlet problem, although the new formulation $\widetilde{N}^{w}S^{w}$ requires fewer iterations than $S^{w}$, since the application of the operator $S^{w}$ is significantly less expensive than the application of the operator $\widetilde{N}^{w}$.

A number of additional techniques are proposed to accurately evaluate the elastic hypersingular integral operators $N^{w}$ and $\widetilde{N}^{w}$. For closed surface scattering problems in elasticity, the novel and exact regularized formulation presented in [6, 33] show that the hypersingular operator in two dimensions can be transformed into a composition of weakly-singular integrals and tangential-derivative operators that involve the Günter derivative and integration-by-parts. For the weighted hypersingular operators $N^{w}$ and $\widetilde{N}^{w}$ in the open-arc case, thanks to the edge-vanishing weight function $w$, the results in the closed-surface case can be extended to the open arc case since all singular terms arising from the integration-by-parts calculation can be eliminated. Additionally, the tangential derivative evaluations, which are approximated in our method by means of FFTs, has proven more efficient than the alternative treatment [20].

The remainder of this paper is organized as follows. Section 2 describes the Dirichlet and Neumann problems of elastic scattering by open surfaces, along with the associated classical boundary integral formulations and including a brief discussion of certain associated challenges. Section 3.1 introduces new weighted operators that explicitly account for the known singular character of the solutions, and Section 3.2 presents our novel investigation of the closed-surface Calderón formula, including numerical verifications of a polynomial-compactness theoretical result. Section 3.3 then presents a numerical examination of the Calderón relation for open arcs in light of the eigenvalue distributions obtained for non-trivial open arc problems. An exact regularized formulation for the hyper-singular operator is presented in Section 3.4. The high order quadrature rules we use for evaluation of the new integral operators are described in Section 4.1. Numerical demonstrations presented in Section 4.2, for both low and high frequencies and for various geometries, demonstrate the high-accuracy and high-order of convergence enjoyed by the proposed approach, as well as the reduced numbers of GMRES linear-algebra iterations required by the algorithm for convergence.

2 Preliminaries

2.1 The elastic scattering problem

Let $\Gamma$ denote a smooth open arc in the plane $\mathbb{R}^{2}$. The complement of $\Gamma$ in $\mathbb{R}^{2}$ is occupied by a linear isotropic and homogeneous elastic medium characterized by the Lamé constants $\lambda,\mu$ satisfying $\mu>0$, $\lambda+\mu>0$, and by the mass density $\rho>0$. Suppressing the time-harmonic dependence $e^{-i\omega t}$ in which $\omega>0$ is the frequency, the displacement field is the solution of the time-harmonic Navier equation

[TABLE]

together the appropriate Kupradze radiation condition [21] at infinity. On $\Gamma$ the solution is assumed to satisfy either the Dirichlet boundary condition

[TABLE]

or the Neumann boundary condition

[TABLE]

Here $\Delta^{*}$ is the Lamé operator defined by

[TABLE]

and $T=T(\partial,\nu)$ is the traction (or stress) operator on the boundary defined as

[TABLE]

in which $\nu$ is the unit normal to the boundary $\Gamma$ and $\partial_{\nu}:=\nu\cdot\nabla$ is the normal derivative.

Along with (2.2), we define the modified (unphysical) traction operator

[TABLE]

where $\widetilde{\lambda}+\widetilde{\mu}=\lambda+\mu$. It easily follows that $\widetilde{T}=T$ when $\widetilde{\lambda}=\lambda,\widetilde{\mu}=\mu$. The selection

[TABLE]

which is made throughout this article, is justified in 3.2 and subsequent sections.

2.2 Boundary integral equations

It follows from potential theory that the solution of (2.1) under Dirichlet and Neumann boundary conditions can be expressed in terms of single- and double-layer potentials,

[TABLE]

and

[TABLE]

respectively. Here, calling $\gamma_{k_{t}}(x,y)$ the fundamental solution of the Helmholtz equation in ${\mathbb{R}}^{2}$ with wave number $k_{t}$,

[TABLE]

and letting

[TABLE]

denote the wave number of the compressional and shear waves for isotropic dymanic elasticity, respectively, where

[TABLE]

$E(x,y)$ denotes the fundamental displacement tensor for the Navier equation in ${\mathbb{R}}^{2}$:

[TABLE]

As is known, using the single-layer and hypersingular operators

[TABLE]

and

[TABLE]

the Dirichlet and Neumann problems reduce to the boundary integral equations

[TABLE]

It is well known that, as demonstrated in Fig.1, the eigenvalues of the single-layer and hypersingular integral operators in equations (2.10) accumulate at zero and infinity, respectively. As a result (and as illustrated in numerical tests in section 4.2), the solution of these equations by means of Krylov-subspace iterative solvers such as GMRES generally requires large numbers of iterations. In addition, as discussed in section 3.1, the solutions $\phi$ and $\psi$ of equations (2.10) are not smooth at the end-points of $\Gamma$, and, thus, they give rise to low order convergence (and require high discretization of the densities for a given accuracy) unless such singularities are appropriately treated.

3 Weighted operators and new integral solvers

3.1 Regularity and singular behavior at the edge

The singular character of the solutions of scattering by open arcs is well documented [16]. In particular, it is known that $\phi,\varphi$ can be expressed in the forms

[TABLE]

where $d$ denotes the distance to the edge, $\xi_{1}$ and $\xi_{2}$ are smooth cut-off functions, and where the functions $\eta_{1}$ and $\eta_{2}$ are somewhat smoother than $\phi$ and $\varphi$. If the curve itself and the boundary functions $F$ and $G$ are infinitely differentiable, it follows that

[TABLE]

where $\alpha$ and $\beta$ are infinitely differentiable functions throughout $\Gamma$, up to and including the endpoints. Thus, the singular character of these solutions is fully characterized by the factors $d^{1/2}$ and $d^{-1/2}$ in equation (3.2).

In view of the regularity results (3.2), we introduce a positive integral weight $w(x)>0$ with asymptotic behavior $w\sim\sqrt{d}$ at the edge (by which it is implied that the quotient $w/\sqrt{d}$ is infinitely differentiable up to the edge), and we define the weighted operators

[TABLE]

so that equations (2.10) for the Dirichlet and Neumann problems may be expressed in the forms

[TABLE]

respectively. In what follows a smooth parameterizations $x(t)=(x_{1}(t),x_{2}(t))$ defined in the interval $[-1,1]$ is used for the curve $\Gamma$, and a canonical choice is made for the weight $w$: $w(x(t))=\sqrt{1-t^{2}}$.

3.2 Calderón relation for closed surfaces

It is well known that for acoustic scattering by closed smooth surfaces there holds the Calderón formula

[TABLE]

where {\color[rgb]{0.000,0.000,0.000}S_{a,c},K_{a,c}^{*}} and {\color[rgb]{0.000,0.000,0.000}N_{a,c}} are the corresponding single-layer, adjoint of double-layer and hyper-singular boundary integral operators of acoustic scattering by closed surfaces, respectively, and {\color[rgb]{0.000,0.000,0.000}K_{a,c}^{*}} is compact. For elastic scattering by closed surfaces, we can obtain a similar Calderón formula for the single-layer and hyper-singular boundary integral operators $S,N$. However, as discussed below, the adjoint $K^{*}$ of the double-layer operator is non-compact—a fact that makes the elastic scattering problem more challenging than its acoustic counterpart. But as indicated in what follows, $NS$ can still be viewed as a compact perturbation of a multiple of the identity operator.

The operator {\color[rgb]{0.000,0.000,0.000}K^{*}} is given by

[TABLE]

Utilizing the artificial traction operator (2.3), we additionally define the operators

[TABLE]

For the special choice (2.4) of the constants $\widetilde{\lambda}$ and $\widetilde{\mu}$, it is known [18, 21] that the kernel of $\widetilde{K}^{*}$ is weakly-singular, which implies that that $\widetilde{K}^{*}$ is a compact operator. However, the kernel of $K^{*}$ is strongly singular and, consequently, the operator $K^{*}$ itself is not compact. As indicated above, this presents a singificant difficulty in our context which can, however, be resolved by appealing to a recent result [4]in the elasto-static context—as indicated in what follows.

Let $K_{0}^{*}$ denote the adjoint of the elastic double-layer operator in the zero-frequency case $\omega=0$. The spectral properties of this kind of operator are best known in the electrostatic case (which involves the Laplace equation), where they relate to plasmonics and cloaking by anomalous localized resonances, which occur at eigenvalues and at the accumulation point of eigenvalues, respectively. The recent contribution [4] shows that cloaking by anomalous localized resonance also occurs in the elastic case. A fundamental result in that paper states that $(K_{0}^{*})^{2}-C_{\lambda,\mu}^{2}I$ is a compact operator (and, in particular, that $K_{0}^{*}$ is polynomially compact), where $C_{\lambda,\mu}$ is a constant that depends on the Lamé parameters:

[TABLE]

This result is useful in our context. Indeed, noting that $K^{*}-K_{0}^{*}$ has a weakly-singular kernel it follows that $K^{*}-K_{0}^{*}$ is a compact operator, and we obtain

[TABLE]

is compact. We thus obtain the elastic Calderón formulae for closed smooth surfaces

[TABLE]

where $K^{(1)},K^{(2)}$ are compact operators.

To verify the above results numerically, we consider the problem of elastic scattering by a circular scatterer of radius one, and we choose $\rho=1$, $\mu=1$, $\omega=50$. Fig. 2 displays the eigenvalue distribution of $NS$ and $\widetilde{N}S$ for various values of $\lambda$. (The eigenvalue computation was based on the exact regularized formulations for the hypersingular integral operator given in section 3.4, implemented by means of the high-order Nyström methodology [15] together with FFT for evaluation of tangential derivatives.) Choosing a large enough number of descritization points the eigenvalues of $NS$ and $\widetilde{N}S$ are seen to accumulate at $(-1/4+C_{\lambda,\mu}^{2},0)$ and $(-1/4,0)$, respectively, as predicted by the result embodied in Equation (3.5).

3.3 Numerical study of the Calderón relation for open arcs

Given the smoothness of the solutions of the equations arising from the weighted operators $S^{w}$, $N^{w}$ and $\widetilde{N}^{w}$, and in light of equation (3.5), it is reasonable to consider the composite operator $N^{w}S^{w}$ as a possible basis for solution of open-arc problems. The operator $\widetilde{N}^{w}S^{w}$ could also employed in the Dirichlet case. Figure (3) presents the eigenvalue distribution for discrete versions of the operators $N^{w}S^{w}$ and $\widetilde{N}^{w}S^{w}$—which were obtained numerically for a spiral-shaped open arc using the spectrally accurate quadrature rules described in section 4.1. It can be see from Fig. 3 that the eigenvalues of $N^{w}S^{w}$ and $\widetilde{N}^{w}S^{w}$ are bounded away from zero and infinity. A theoretical analysis of the Calderón relation for acoustic diffraction by open arcs is presented in [9, 22]. The corresponding analysis for elastic problems is left for future work.

In view of this numerical evidence and theoretical background on related problems we suggest that the open arc problems can be solved effectively by means of the integral equations

[TABLE]

The smooth solutions $\alpha$ and $\beta$ of these equations are related to the singular solutions of (2.10) via the relations $\phi=\alpha/w$ and $\psi=\beta\cdot w$.

3.4 Regularized formulation for hyper-singular operator

In this section, we derive an accurate regularized formulation for the hyper-singular boundary integral operator $\widetilde{N}^{w}$ in light of the techniques presented in [6, 33] for the closed-surface case. The regularized formulation for the hyper-singular boundary integral operator $N^{w}$ can be obtained directly by letting $\widetilde{\mu}=\mu$.

3.4.1 Closed-surface case

The artificial traction operator (2.3) can be expressed in the form

[TABLE]

where the operator $M(\partial,\nu)$, whose elements are also called the Günter derivatives, is defined as

[TABLE]

In fact, the operator $M$ in 2D is the multiplication of the tangential derivative and a constant matrix, i.e.,

[TABLE]

Then following the technique in [33] which is part of the conversion of traction and normal derivatives into tangential derivatives, it can be established that

[TABLE]

where

[TABLE]

We omit the proof here.

3.4.2 Open-arc case

For the open-arc case, due to the smooth boundary-vanishing weight $w$ in $\widetilde{N}^{w}$, the integration-by-parts formula

[TABLE]

holds. We thus, can obtain from the steps in [33] the regularized formulation for the hyper-singular boundary integral operator $\widetilde{N}^{w}$ as follows

[TABLE]

4 Numerical experiments

4.1 Numerical implementations

In this section, we present spectral quadrature rules for the operators $S^{w}$, $N^{w}$ and $\widetilde{N}^{w}$ which give rise to an efficient and accurate solver for the general elastic and thermoelastic open arc diffraction problems.

As indicated in Section 3.1, without loss of generality we may use a smooth parameterization $x(t)=(x_{1}(t),x_{2}(t)),t\in[-1,1]$ of $\Gamma$ satisfying $|x^{\prime}(t)|=|dx(t)/dt|>0$, we choose the weight $w$ as $w(x(t))=\sqrt{1-t^{2}}$. Hence, the operator $S^{w}$ give rise to

[TABLE]

Introducing the change of variables $t=\cos\theta$ and $\tau=\cos\theta^{\prime}$ and defining $\nu_{\theta}=\nu_{x(\cos\theta)}$, we obtain the periodic weighted single-layer operator (using the same notation)

[TABLE]

It follows from the definition of the fundamental solution and the series expansions of Bessel functions [1, 15], we can obtain the decomposition

[TABLE]

where

[TABLE]

with $r(t,\tau)=|x(t)-x(\tau)|$ and

[TABLE]

When $t=\tau$, we have

[TABLE]

and

[TABLE]

Use of the Chebyshev points $\left\{\theta_{n}=\frac{\pi(2n+1)}{2N}\right\},n=0,1,\cdots,N-1$ gives rise to a spectrally convergent cosine representation for smooth, $\pi$-periodic and even function $\widetilde{\phi}$ as

[TABLE]

It is known from the diagonal property of Symm’s operator in the cosine basis $e_{n}(\theta)=\cos(n\theta)$ [11, 32] that

[TABLE]

Applying equation 4.2 to each term of expansion (4.1) we can obtain the well-known spectral quadrature rule for the logarithmic kernel

[TABLE]

where

[TABLE]

Together with the trapezoidal integration for smooth function we therefore obtain the spectrally quadrature approximation of the operator $S_{o}$ that

[TABLE]

Then we can evaluate $S^{w}[\phi]$ in the sets of quadrature points $\{\theta_{n},n=0,\cdots,N-1\}$ by means of a matrix-vector multiplication involving the matrix ${\bf S}^{w}$ whose elements are defined by

[TABLE]

in which the quantities $R_{j}^{(N)}(\theta_{i})$ can be evaluated through

[TABLE]

where

[TABLE]

while this expression can significantly reduce the computational cost of ${\bf S}_{o}^{w}$.

In order to evaluate $N^{w}$ and $\widetilde{N}^{w}$, we use the regularized formulations given in section 3.4. In (3.8), the first term on the right hand side can be evaluated by means of a rule analogous to the single-layer operator of acoustic [9]). The other terms take one of the following forms

[TABLE]

where $S_{j},j=1,2,3$, whose kernels are of logarithmic type, can be evaluated analogous to $S^{w}$, and

[TABLE]

Then we approximate the quantity $T_{0}[\widetilde{\phi}]$ by means of term per term differentiation of the sine expansion of the function $\widetilde{\phi}(\theta)\sin\theta$ (which can itself be produced efficiently by means of an FFT). The quantity $D_{0}[\widetilde{\phi}]$ can be evaluated by invoking classical FFT-based Chebyshev differentiation rules [30]).

4.2 Numerical examples

The numerical examples presented in what follows were obtained by means of a Matlab implementation of the quadrature rules introduced in section 4.1 for numerical evaluation of the operators $S^{w}$, $N^{w}$ and $\widetilde{N}^{w}$ in conjunction with the iterative linear algebra solver GMRES. In all cases the (maximum) errors reported were evalueated by comparisons with exact or highly-resolved numerical solutions. Select the elastic parameters as $\rho=1$, $\mu=1$ and $\lambda=2$. We call the equations in (3.3) and (3.6) as $\mbox{Dir}(S^{w})$, $\mbox{Neu}(N^{w})$ and $\mbox{Dir}(N^{w}S^{w})$, $\mbox{Dir}(\widetilde{N}^{w}S^{w})$, $\mbox{Neu}(N^{w}S^{w})$, respectively.

We first demonstrate the accuracy of the regularized formulation given in section 3.4. To see this, let $\Gamma$ be a unit circle and the exact solutions for elastic and thermoelastic problems are given by $u=\nabla_{x}H_{0}^{(1)}(k_{p}|x-z_{0}|)$ and $u=E_{12}(x,z_{0})$, $p=E_{22}(x,z_{0})$, respectively where $z_{0}=(0,0.5)$. We apply the Nyström method [15] for the discretization of the operators $S$, $N$ and $\widetilde{N}$. Table 1 shows the spectral (exponentially-fast) convergence which further demonstrate the accuracy of our regularized formulations.

Now we consider the elastic scattering by open arcs. Let $u^{inc}$ be a plane incident pressure wave $u^{inc}=d_{inc}\exp(ik_{p}x\cdot d_{inc})$ where $d_{inc}=(\cos\theta_{inc},\sin\theta_{inc})$ and $\theta_{inc}$ is the incident direction. This means that the boundary data is given by

[TABLE]

To demonstrated the high-order character of the algorithm we consider the elastic scattering by the ’Spiral-Shaped Arc’ characterized by the parameterization

[TABLE]

Choose $\theta_{inc}=\pi/4$. We present the numerical errors of the solution produced by means of the operators $S^{w}$, $N^{w}$, $\widetilde{N}^{w}S^{w}$ and $N^{w}S^{w}$ for the elastic scattering problems which demonstrate the spectral (exponentially-fast) convergence.

In Table 3 and 4, we present the iteration numbers required to achieve the GMRES tolerance $10^{-5}$ for the scattering by a flat strip $[-1,1]$ and the exponential spiral mentioned above. It can be seen that for Neumann problems, $\mbox{Neu}(N^{w})$ requires very large number of iterations as the frequency grows and, thus, the computing times required by the low-iteration equation $\mbox{Neu}(N^{w}S^{w})$ are significantly lower than those required by $\mbox{Neu}(N^{w})$. But for the Dirichlet problem, $\mbox{Dir}(N^{w}S^{w})$ requires almost the same number of iterations as $\mbox{Neu}(N^{w}S^{w})$ which is significantly larger that $\mbox{Dir}(S^{w})$ for the spiral-shaped arc. But the corresponding equation $\mbox{Dir}(\widetilde{N}^{w}S^{w})$ requires fewer iterations than $\mbox{Dir}(S^{w})$ for both flat strip and spiral-shaped arc. We present the corresponding total time required by the solver to find the solution once the needed matrixes for iteration, $S^{w}$, $S_{1}$, $S_{2}$ and $S_{3}$ mentioned in section ?, are stored. It should be pointed out that the total computing cost of the equations $\mbox{Dir}(\widetilde{N}^{w}S^{w})$ and $\mbox{Dir}(N^{w}S^{w})$ is generally higher than $\mbox{Dir}(S^{w})$ since the application of the operator in $\mbox{Dir}(S^{w})$ is significantly less expensive than the application of operator in $\mbox{Dir}(\widetilde{N}^{w}S^{w})$ and $\mbox{Dir}(N^{w}S^{w})$. In addition, from the regularized formulation, the cost for the generation of the matrixes $\widetilde{N}^{w}$ and $N^{w}$ (except the computing of tangential derivatives) is about 5 times of that for $S^{w}$.

In order to obtain an indication of the manner in which an open arc problem can be viewed as a limit of closed-curve problems and provide an independent verification of the validity of our solvers, we consider a test case in which the flat strip $[-1,1]$ is viewed as the limit as $a\rightarrow 0$ of the family of closed curves $x(t)=(\cos t,a\sin t)$. It is known that the scattered field $u$ admits the asymptotic behavior

[TABLE]

Then for the Dirichlet problem, the P-part $u_{p}^{\infty}$ and S-part $u_{p}^{\infty}$ of the far-field pattern of $u$ are given by

[TABLE]

Using the closed-curve Nyström method we evaluate the scattered field of Dirichlet problems for values of $a$ approaching 0 and present the P-part and S-part of the far-field pattern side-by-side the corresponding far-field pattern for the limiting open arc as produced by the $S^{w}$-based open-arc solver. Clearly, the closed-curve and open-arc solutions are quite close to each other.

5 Conclusion

We have introduced new integral solvers and associated numerical algorithms for the elastic scattering by open arcs with Dirichlet or Neumann boundary condition in two dimensions. The new methods enjoy spectral convergence and reduce the number of GMRES iterations consistently across various geometries and frequency regimes. In particular, the new formulation is highly beneficial for the Neumann problem, giving rise to order-of-magnitude improvements in computing time over the original hypersingular formulation. Theoretical investigation of the Calderón formula in open-arc case and the generalization of the method enabling efficient solution of problems of elastic and thermoelastic scattering by open-surfaces in three dimensions are left for future work.

Acknowledgments

OB gratefully acknowledge support by NSF, AFOSR and DARPA through contracts DMS- 1411876 and FA9550-15-1-0043 and HR00111720035, and the NSSEFF Vannevar Bush Fellowship under contract number N00014-16-1-2808. LX is partially supported by a Key Project of the Major Research Plan of NSFC (No. 91630205), and a NSFC Grant (No. 11771068).