Boundary integral equation methods for the elastic and thermoelastic waves in three dimensions

TL;DR

This paper develops new regularized boundary integral equation formulations for elastic and thermoelastic wave problems in three dimensions, enabling more accurate and efficient numerical solutions using boundary element methods.

Contribution

It introduces novel regularized formulations for hyper-singular boundary integral operators in 3D elastic and thermoelastic wave equations, simplifying numerical computations.

Findings

Regularized formulations reduce singularity complexity.

Numerical examples demonstrate high accuracy.

Effective boundary element method implementation.

Abstract

In this paper, we consider the boundary integral equation (BIE) method for solving the exterior Neumann boundary value problems of elastic and thermoelastic waves in three dimensions based on the Fredholm integral equations of the first kind. The innovative contribution of this work lies in the proposal of the new regularized formulations for the hyper-singular boundary integral operators (BIO) associated with the time-harmonic elastic and thermoelastic wave equations. With the help of the new regularized formulations, we only need to compute the integrals with weak singularities at most in the corresponding variational forms of the boundary integral equations. The accuracy of the regularized formulations is demonstrated through numerical examples using the Galerkin boundary element method (BEM).

Figures

Tables

Table 1. Table 1: Numerical errors Error with respect to the meshsize h ℎ h .

$h$	Error	Order
0.4880	1.46E-3	–
0.3871	7.51E-4	2.87
0.2668	2.95E-4	2.51
0.1913	1.33E-4	2.39
0.1005	3.01E-5	2.31

Table 2. Table 2: Numerical errors Error ~ ~ Error \widetilde{\mbox{Error}} with respect to the meshsize h ℎ h .

$h$	$\widetilde{\text{Error}}$	Order
0.4880	5.32E-4	–
0.3871	2.22E-4	3.77
0.2668	1.07E-4	1.96
0.1913	3.59E-5	3.28
0.1005	2.61E-6	4.07

Equations

$${\Delta }^{\ast }u+\rho {\omega }^{2}u$$

$$T(\partial ,\nu )u$$

$$\underset{r\to \infty }{\lim }r\left(\frac{\partial u_t}{\partial r}-i{k}_t{u}_t\right)=0,r=|x|,t=p,s,$$

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

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

$$u_p=-\frac{1}{k_p^2}\text{mbox}grad\text{mbox}divu,u_s=\frac{1}{k_s^2}\text{mbox}curl\text{mbox}curlu,$$

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

$$c_p:=\sqrt{\mu /\rho },c_s:=\sqrt{(\lambda +2\mu )/\rho }.$$

$$\gamma =(\lambda +\frac{2}{3}\mu )\alpha ,\eta =\frac{T_0\gamma }{{\lambda }_0},$$

$${\Delta }^{\ast }u+\rho {\omega }^{2}u-\gamma \nabla p$$

$$\Delta p+qp+i\omega \eta \nabla \cdot u$$

$$LU=0,L:=\left[\begin{array}{cc}(\mu \Delta +\rho {\omega }^2)I_3+(\lambda +\mu )\nabla \nabla \cdot & -\gamma \nabla \\ q\eta \kappa \nabla \cdot & \Delta +q\end{array}\right],U=(u^{\top },p{)}^{\top }.$$

$$\tilde{T}(\partial ,\nu )U:=\left[\begin{array}{cc}T(\partial ,\nu )& -\gamma \nu \\ 0& {\partial }_{\nu }\end{array}\right]U=F.$$

$$u=u^1+u^2+u^s,p=p^1+p^2,$$

$$\Delta u^i+k_i^2{u}^i=0,i=1,2\text{mbox}and\Delta u^s+k_s^2{u}^s=0$$

$$\text{mbox}curl{u}^i=0i=1,2\text{mbox}and\text{mbox}div{u}^s=0,$$

$$\Delta p^i+k_i^2{p}^i=0,i=1,2.$$

$$k_1^2+k_2^2=q(1+ϵ)+k_p^2,k_1^2{k}_2^2=q{k}_p^2,$$

$$k_1$$

$$k_2$$

$$C_{\pm }=i(1+ϵ)c_p^2\pm (1+i)c_p\sqrt{2\omega \kappa }.$$

$$u^j=o(r^{-1}),{\partial }_{x_i}{u}^j=O(r^{-2}),$$

$$p^j=o(r^{-1}),{\partial }_{x_i}{p}^j=O(r^{-2}),$$

$$u^s=o(r^{-1}),r({\partial }_r{u}^s-i{k}_s{u}^s)=O(r^{-1}).$$

$$u(x)=({\mathcal{D}}_s\phi )(x):={\int }_{\Gamma }(T({\partial }_y,{\nu }_y)E(x,y){)}^{\top }\phi (y)d{s}_y,\forall x\in {\Omega }^c,$$

$$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],x\ne y.$$

$${\gamma }_{k_t}(x,y)=\frac{\exp (i{k}_t|x-y|)}{4\pi |x-y|},x\ne y,t=p,s.$$

$$W_s\phi (x)=-f,x\in \Gamma ,$$

$$W_{s}u(x):=-\underset{z\to x\in \Gamma ,z\notin \Gamma }{\lim }T({\partial }_z,{\nu }_x){\int }_{\Gamma }(T({\partial }_y,{\nu }_y)E(z,y){)}^{\top }u(y)d{s}_y,x\in \Gamma .$$

$$⟨W_s\phi ,v⟩=-⟨f,v⟩\text{mbox}forallv\in H^{1/2}(\Gamma {)}^3.$$

$$U(x)=(\tilde{\mathcal{D}}\Psi )(x):={\int }_{\Gamma }({\tilde{T}}^{\ast }({\partial }_y,{\nu }_y){\tilde{E}}^{\top }(x,y){)}^{\top }\Psi (y)d{s}_y,\forall x\in {\Omega }^c,$$

$$\tilde{E}(x,y)=\left[\begin{array}{cc}{E}_{11}(x,y)& E_{12}(x,y)\\ E_{21}^{\top }(x,y)& E_{22}(x,y)\end{array}\right],$$

$$E_{11}(x,y)$$

Full text

Keywords: Elastic wave, thermoelastic wave, hyper-singular boundary integral operator, regularized formulation, boundary element method

1 Introduction

In this paper, we apply the BIE method to solve the three dimensional time-harmonic elastic and thermoelastic scattering problems that are of great importance in many fields of applications such as geophysics, seismology, non-destructive testing and material sciences, to name a few. We are interested in the wave scattering by a bounded impenetrable obstacle immersed in an infinite isotropic solid. Based upon the assumption that any deformation of the elastic medium occurs under a constant temperature, the displacement in the solid can be modeled by the Navier equation together with an appropriate radiation condition at infinity ([22]). If considering the temperature fluctuations caused by the dynamic deformation, one arrives at the Biot system of equations ([5, 22, 33]) describing the interaction of the temperature and displacement fields.

When considering the wave propagation in an unbounded domain, one could use the so-called Dirichlet-to-Neumann (DtN) map (or non-reflecting boundary condition) on a closed artificial boundary which decomposes the exterior region into two parts. After that, the original scattering problem could be solved on the bounded region. The DtN map for the elastic wave has been used for numerical simulations in the open literatures ([14, 17]) and some properties of the DtN map have been investigated in [1, 24, 25]. For the thermoelastic wave, the explicit formulation of the DtN map is still unknown. We refer to [8, 9, 12, 13] for the mathematical analysis of the thermoelastic wave. The BIE method is another conventional numerical method for solving the scattering problems, and it has been widely used in acoustics, electromagnetics, elastodynamics and thermoelastodynamics ([4, 2, 7, 10, 11, 18, 19, 20, 21, 23, 27, 26, 28, 31, 34, 35, 36, 37, 38]). The BIE method takes some advantages over domain discretization methods, including that the boundary integral representation of the solution fulfills the radiation condition naturally, the dimension of the computational domain is reduced by one, and to name a few. Various numerical techniques, including the Galerkin scheme, the Nyström method, the fast multipole method, and the spectral method etc., have been developed for the efficient transformation of the BIE into a linear system in the past decades. In this work, we will use the Galerkin scheme ([18, 19, 34]) for the numerical solutions, and its advantages include the availability of mathematical convergence analysis allowing $h$-$p$ approximations, and particularly the strength on dealing with the hyper-singularities in the boundary integrals.

During the application of the BIE method, there needs for the use of hyper-singular BIO in many situations, containing the removal of the pollution of eigenfrequencies for the BIE ([6]), and the solution of the Neumann boundary value problem using the Fredholm integral equation of the first kind ([15, 32]), etc.. Theoretical analysis indicates that the hyper-singular BIO is equivalent to the Hadamard finite part of a hyper-singular integral ([19]) which is usually difficult to be calculated accurately. One usually needs additional treatments in numerics for the correct evaluation of the classically non-integrable boundary integrals arising from the hyper-singular BIO. There are already existing many works ([15, 26, 16, 19, 29, 30, 27, 38, 2]) on this issue, and the main idea consists of rewriting the hyper-singular BIO in terms of a composition of differentiation and weakly singular operators for the Laplace equation, the Helmholtz equation, the time-harmonic Navier equation and the Lamé equation. This composition, in fact, is a regularization procedure ([19, 38]) of the hyper-singular distribution, and is useful for the variational formulation and related computational procedures.

In this paper, we consider the three dimensional elastic and thermoelastic scattering problems with Neumann boundary conditions. For each problem, we apply the double-layer potential to represent the solution, and then the original boundary value problem is reduced to a Fredholm boundary integral equation of the first kind with the corresponding hyper-singular BIO. Following the work in [38, 2] for the two-dimensional case, and utilizing the tangential Günter derivative, we derive the new and analytically accurate regularized formulations for the hyper-singular BIO associated with the three dimensional time-harmonic Navier equation and Biot system of linearized thermoelasticity, respectively. As a result, in the corresponding weak forms, all involved integrals are at most weakly-singular. This work is an extension of our work only considering two dimensional elastic waves, and the extension is in fact non-trivial. In the numerical implementations, applying the special local coordinate system given in [34] and the Gauss quadrature rules on the triangle element, we present a semi-analytic strategy to evaluate all the weakly-singular integrals effectively. Although we only consider $C^{2}$ boundary in our numerical tests, the theoretical results actually could be extended to the Lipschitz case in terms of the properties of Günter derivative given in [3]. The convergence analysis of the numerical scheme could be obtained following the standard techniques in [19], and we will omit it in this work.

The rest of this paper is organized as follows. In Section 2, we introduce the exterior elastic and thermoelastic scattering problems, and then describe the BIE and the Galerkin BEM in Section 3. In section 4, we propose the new and analytically accurate regularized formulations for the hyper-singular BIO in three dimensions. Finally, we discuss a semi-analytic strategy for the numerical implementation of the Galerkin scheme and present numerical results of several examples in section 5.

2 Mathematical problems

Let $\Omega\subset{\mathbb{R}}^{3}$ be a bounded, simply connected and impenetrable body with $C^{2}$ boundary $\Gamma=\partial\Omega$. The exterior complement of $\Omega$ is denoted by $\Omega^{c}={\mathbb{R}}^{3}\setminus\overline{\Omega}\subset{\mathbb{R}}^{3}$. Assume that $\Omega^{c}$ is occupied by a linear and isotropic elastic solid characterized by the Lamé constants $\lambda$ and $\mu$ ($\mu>0$, $3\lambda+2\mu>0$) and mass density $\rho>0$. Let $\omega>0$ be the frequency of propagating waves.

2.1 Elastic scattering problem (ESP)

Assume that the temperature is always a constant and suppress the time-harmonic dependence $e^{-i\omega t}$. Then the displacement field $u$ in the solid can be modeled by the following exterior ESP: Given $f\in H^{-1/2}(\Gamma)^{3}$, determine $u=(u_{1},u_{2},u_{3})^{\top}\in H_{loc}^{1}(\Omega^{c})^{3}$ satisfying

[TABLE]

and the Kupradze radiation condition ([22])

[TABLE]

uniformly with respect to all $\hat{x}=x/|x|\in\mathbb{S}^{2}:=\{x\in{\mathbb{R}}^{3}:|x|=1\}$. Here, $\Delta^{*}$ is the Lamé operator defined by

[TABLE]

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

[TABLE]

where $\nu$ is the outward unit normal to the boundary $\Gamma$ and $\partial_{\nu}:=\nu\cdot\mbox{grad}$ is the normal derivative. In (2.3), $u_{p}$ and $u_{s}$ are referred as the compressional wave and the shear wave, respectively, and they are given by

[TABLE]

where the wave numbers $k_{s},k_{p}$ are defined as

[TABLE]

with

[TABLE]

For the uniqueness of the ESP (2.1)-(2.3), we refer to [22, 1].

2.2 Thermoelastic scattering problem (TESP)

Now we consider the temperature fluctuations caused by the dynamic deformation. In this case, the elastic medium $\Omega^{c}$ is additionally characterized by the coefficient of thermal diffusity $\kappa$ and the coupling constants $\gamma$, $\eta$ given by

[TABLE]

respectively, where $\alpha$ is the volumetric thermal expansion coefficient, $T_{0}$ is a reference value of the absolute temperature and $\lambda_{0}$ is the coefficient of thermal conductivity. Denote by $\epsilon:=\gamma\eta\kappa/(\lambda+2\mu)$ the dimensionless thermoelastic coupling constant which assumes ’small’ positive values for most thermoelastic media and $q=i\omega/\kappa$. Suppressing the time-harmonic dependence $e^{-i\omega t}$, the displacement field $u$ and the temperature variation field $p$ can be modeled by the following Biot system of linearized thermoelasticity

[TABLE]

Rewriting (2.4)-(2.5) into a matrix form, we obtain

[TABLE]

On the boundary of the scatterer, we assume the Neumann boundary condition

[TABLE]

It follows ([22]) that the wave field $U$ admits the decomposition

[TABLE]

where the vector displacement fields $u^{1},u^{2},u^{s}$ satisfy the vectorial Helmholtz equations

[TABLE]

with

[TABLE]

and the scalar temperature fields $p^{1}$ and $p^{2}$ satisfy the following scalar Helmholtz equations

[TABLE]

Here, the wave numbers $k_{1},k_{2}$, corresponding to the elastothermal and thermoelastic waves respectively, are the roots of the characteristic system

[TABLE]

for which $\mbox{Im}\,k_{i}\geq 0,i=1,2$. In particular,

[TABLE]

where

[TABLE]

We assume that the scattered field $U$ satisfies the following Kupradze radiation conditions as $r=|x|\rightarrow\infty$ for $i=1,2,3$ and $j=1,2$

[TABLE]

The direct TESP to be considered in this paper is to determine the displacement field $u$ and the temperature variation field $p$ satisfying (2.6), the boundary condition (2.7) and the Kupradze radiation conditions. For given $F\in(H^{-1/2}(\Gamma))^{4}$, we refer to [22, 7] for the uniqueness of the direct problem.

3 Numerical method

In this section, we derive the BIE for solving the ESP and TESP, respectively and give a brief introduction to the Galerkin BEM for the discretization of the derived BIE.

3.1 BIE for ESP

For the ESP, it follows from the potential theory ([22]) that the unknown function $u$ can be represented as

[TABLE]

where $\mathcal{D}_{s}$ is referred as the double-layer potential and $E(x,y)$ is the fundamental displacement tensor of the time-harmonic Navier equation (2.1) in ${\mathbb{R}}^{3}$ taking the form

[TABLE]

In (3.2) and the following, $I$ denotes the $3\times 3$ identity matrix, and $\gamma_{k_{t}}(x,y)$ is the fundamental solution of the Helmholtz equation in ${\mathbb{R}}^{3}$ with wave number $k_{t}$ and takes the form

[TABLE]

Operating with the traction operator on (3.1), taking the limits as $x\to\Gamma$ and applying the jump relations and the boundary condition (2.2) , we arrive at the BIE on $\Gamma$

[TABLE]

where $W_{s}:H^{s}(\Gamma)^{3}\rightarrow H^{s-1}(\Gamma)^{3}(s\geq 1/2)$ is called the hyper-singular BIO defined by

[TABLE]

The standard weak formulation of (3.4) reads: Given $f\in H^{-1/2}(\Gamma)^{3}$, find $\varphi\in H^{1/2}(\Gamma)^{3}$ such that

[TABLE]

Here and in the sequel, $\langle\cdot,\cdot\rangle$ denotes the $L^{2}$ duality pairing between $H^{-1/2}(\Gamma)^{d}$ and $H^{1/2}(\Gamma)^{d}$ for $d\in{\mathbb{Z}}^{+}$.

3.2 BIE for TESP

For the TESP, it follows from the potential theory ([22, 7]) that the unknown function $U$ can be represented as

[TABLE]

where $\widetilde{E}(x,y)$ is the fundamental solution of the Biot system (2.6) in ${\mathbb{R}}^{3}$ given by

[TABLE]

with

[TABLE]

and $\widetilde{T}^{*}(\partial,\nu)$ is the corresponding Neumann operator of the adjoint problem of (2.6) taking the form

[TABLE]

Operating with the operator $\widetilde{T}$ on (3.7), taking the limits as $x\to\Gamma$ and applying the boundary condition (2.7) , we obtain the BIE on $\Gamma$

[TABLE]

where the hyper-singular BIO $\widetilde{W}_{s}$ is defined by

[TABLE]

The standard weak formulation of (3.8) reads: Given $F\in H^{-1/2}(\Gamma)^{4}$, find $\Psi\in H^{1/2}(\Gamma)^{4}$ such that

[TABLE]

Remark 3.1**.**

For the wellposedness of the variational equations (3.6) and (3.10), we refer the readers to [2, 7, 22]. The pollution of eigenfrequencies on the uniqueness can be removed by applying the so-called Burton-Miller formulation, see [2, 6] for example.

3.3 Galerkin BEM

We only propose the Galerkin Scheme for solving ESP, and the corresponding procedure and formulas for TESP are quite similar and will be neglected. Let $\mathcal{H}_{h}$ be a finite dimensional subspace of $H^{1/2}(\Gamma)$. Then the Galerkin approximation of (3.4) reads: Given $f$, find $\varphi_{h}\in\mathcal{H}_{h}^{3}$ satisfying

[TABLE]

In the following, we describe briefly the reduction of the Galerkin equation (3.11) into its discrete linear system of equations.

Let $\Gamma_{h}=\cup_{i=1}^{N}\overline{\tau_{i}}$ be a uniform boundary element mesh of $\Gamma$ where each $\tau_{i}$ is a plane triangle with vertex $x_{i_{1}},x_{i_{2}},x_{i_{3}}$ ordered counter clockwise. Let $\{x_{j}\}_{j=1}^{M}$ be the set of all nodes of the triangulation. Using the reference element

[TABLE]

the point $x\in\tau_{i}$ can be parameterized as

[TABLE]

Let $\{\psi_{j}\}_{j=1}^{M}$ be the set of piecewise linear basis functions. We seek the approximate solution

[TABLE]

where $u_{j}\in{\mathbb{C}}^{3},j=1,\cdots,M$ are unknown nodal values of $u_{h}$ at $x_{j}$. For the boundary value $f$, we interpolated it as

[TABLE]

where $x_{i_{*}}=(x_{i_{1}}+x_{i_{2}}+x_{i_{3}})/3$ is the mid point of the element $\tau_{i}$ and $\{\phi_{i}\}_{i=1}^{N}$ is the set of piecewise constant basis functions. Substituting the interpolation forms into (3.11) and setting $\psi_{j},j=1,\cdots,M$ as test functions, we arrive at a linear system of equations

[TABLE]

where for $k,j=1,\cdots,M$, $i=1,\cdots,N$,

[TABLE]

[TABLE]

[TABLE]

4 Regularized formulation for the hyper-singular BIO

In this section, we derive the new regularized formulations for the hyper-singular BIOs $W_{s}$ and $\widetilde{W}_{s}$ in three dimensions such that the coefficient matrix ${\bf A}_{h}$ ( or $\widetilde{\bf A}_{h}$) can be evaluated in a more effective and accurate way. More precisely, using the derived regularized formulations, only classically integrable and weakly-singular integrals are involved in the weak forms of $W_{s}$ and $\widetilde{W}_{s}$. Before doing this, we introduce the hyper-singular BIO associated with the Helmholtz equation and the Günter derivatives.

4.1 Hyper-singular BIO for acoustic scattering problem

Consider the Helmholtz equation

[TABLE]

with wave number $k>0$. Denote by $V_{f}:H^{s-1}(\Gamma)\rightarrow H^{s}(\Gamma)$ and $W_{f}:H^{s}(\Gamma)\rightarrow H^{s-1}(\Gamma)$, $s\geq 1/2$ the single-layer and hyper-singular BIO defined by

[TABLE]

respectively. It follows from Lemma 1.2.2 in [19] that the hyper-singular BIO $W_{f}$ can be expressed as

[TABLE]

4.2 Günter derivatives

Now we describe the Günter derivatives that play essential roles in the proof of our main results. Define the operator $M(\partial,\nu)$, whose elements are also called the Günter derivatives, as

[TABLE]

Then the traction operator can be rewritten as

[TABLE]

A direct calculation yields

[TABLE]

Then

[TABLE]

which implies that

[TABLE]

The properties of the operator $M(\partial,\nu)$ ([3, 19, 22]) shows that for any scaler fields $p,q$, vector fields $u,v$ and tensor field $E$, there hold the Stokes formulas

[TABLE]

and

[TABLE]

4.3 Hyper-singular BIO for ESP

We now investigate the operator $W_{s}$. Following the results in [27, 38], we have for $x\neq y$,

[TABLE]

and

[TABLE]

Using the Günter derivatives, the hyper-singular operator $W_{s}$ can be rewritten as

[TABLE]

where

[TABLE]

Theorem 4.1**.**

The hyper-singular BIO $W_{s}$ in three dimensions can be expressed alternatively as

[TABLE]

*where $J_{\nu_{x},\nu_{y}}=\nu_{y}\nu_{x}^{\top}-\nu_{x}\nu_{y}^{\top}$. *

Proof.

See A. ∎

4.4 Hyper-singular BIO for TESP

Now we consider the operator $\widetilde{W}_{s}$. Note that the hyper-singular kernel of $\widetilde{W}_{s}$ is

[TABLE]

where

[TABLE]

For $U=(u^{\top},p)$ and $V=(v^{\top},q)$, we have

[TABLE]

where

[TABLE]

Lemma 4.2**.**

For $x\neq y$, it follows that

[TABLE]

and

[TABLE]

Proof.

See B. ∎

We first investigate the term $W_{1}u$. Observe that the first term in $W_{11}$ is consistent with the the kernel of hyper-singular BIO $W_{s}$. For $W_{1}u$, we have the following regularized formulation.

Theorem 4.3**.**

The hyper-singular operator $W_{1}$ can be expressed alternatively as

[TABLE]

Here, the constants $C_{i},i=1,2,3$ are given by

[TABLE]

and

[TABLE]

Proof.

See C. ∎

Next we investigate the terms $W_{2}p$ and $W_{3}u$. We have

Theorem 4.4**.**

The hyper-singular operators $W_{2}$ and $W_{3}$ can be expressed as

[TABLE]

and

[TABLE]

*respectively. *

Proof.

It follows that

[TABLE]

and

[TABLE]

which further implies (4.15) by the Stokes formulas (4.6). The proof of (4.16) is similar and we omit it here. ∎

Finally, we investigate the term $W_{4}p$. From the results for acoustic wave (4.1), we immediately conclude that

Theorem 4.5**.**

The hyper-singular operator $W_{4}$ can be expressed as

[TABLE]

Remark 4.6**.**

It can be easily verified from the Stokes formulas of Günter derivatives that using the proposed regularized formulations, all the integrals in the corresponding weak forms of $W_{s}u$ and $\widetilde{W}_{s}U$ are at most weakly-singular.

5 Numerical tests

In this section, we present several numerical examples to demonstrate the accuracy of the proposed scheme solving the exterior ESP and TESP. We now take the ESP as the model to describe the method for numerical implementations.

5.1 Numerical implementations

Using the weak form of the regularized formulation (4.10), it follows that (3.13) can be retreated as

[TABLE]

in which the entries of the matrix $M_{\psi_{j},\psi_{k}}$ are given by

[TABLE]

In (5.1), all the integrals are at most weakly-singular. It can be obtained from decompositions that the weakly-singular kernels in (5.1) are of types

[TABLE]

To compute the weakly singular integrals efficiently, we apply the special local coordinate system given in [34] to the boundary element $\tau$ with vertex $x_{(1)},x_{(2)},x_{(3)}$. Define the unit vector $r_{\tau,2}=(x_{(3)}-x_{(2)})/(|x_{(3)}-x_{(2)}|)$. Set $q_{1}=(x_{(2)}-x_{(1)})\cdot r_{\tau,2}$ and $q_{2}=(x_{(3)}-x_{(1)})\cdot r_{\tau,2}$. Then the intersection point $x^{*}$ can be determined by $x^{*}=x_{(3)}-q_{2}r_{\tau,2}$ or equivalently, $x^{*}=x_{(2)}-q_{1}r_{\tau,2}$. Define another unit vector $r_{\tau,1}=(x^{*}-x_{(1)})/(p_{\tau})$, $p_{\tau}=|x^{*}-x_{(1)}|$ and the proportional coefficients $\alpha_{i}=q_{i}/s_{\tau},\quad i=1,2$. Thus, the boundary element $\tau$ can be parameterized as

[TABLE]

Then for $x,y\in\tau$, $|x-y|^{2}=(p_{x}-p_{y})^{2}+(q_{x}-q_{y})^{2}$. The outward unit norma $\nu_{\tau}$ to the element $\tau$ is determined by $\nu_{\tau}=r_{\tau,1}\times r_{\tau,2}$. Moreover, the piecewise linear basis function $\psi_{(1)}$ for the vertex $x_{(1)}$ on $\tau$ can be formulated as $\psi_{(1)}(x)=\psi_{(1)}(x(p_{x},q_{x}))=(p_{\tau}-p_{x})/(p_{\tau}),x\in\tau$, and $\nabla_{x}\psi_{(1)}(x)=-r_{\tau,1}/p_{\tau}$. In addition, using the above parametrisation we have

[TABLE]

or

[TABLE]

Now we present the main computing strategy of the numerical implementation. Set $\tau=\tau_{i}$. Corresponding to the piecewise linear basis function $\psi_{i_{m}}(x),m=1,2,3$ on $\tau$, set

[TABLE]

Using this reorder strategy and the local coordinate system, $\psi_{i_{m}}(x)=\psi_{(1)}(x(p_{x},q_{x}))$ on $\tau$. Therefore,

[TABLE]

are all constants.

The nonsingular integrals involved in (5.1) can be approximated by Gaussian quadrature for triangular elements and we only need to consider the following integrals

[TABLE]

which can be numerically computed following the steps described in [34] in a semi-analytic sense.

5.2 Numerical examples

In the numerical tests, the direct solver ’$\setminus$’ in Matlab is employed for solutions of the linear system (3.12). The impenetrable obstacle $\Omega$ is set to be a unit ball (see Figure 2 (a)) or star-like (see Figure 2 (b)) with radial function

[TABLE]

For these two obstacles, the origin $O$ is in $\Omega$. In our numerical tests, we first compute the unknown potentials $\varphi_{h}$ and $\Psi_{h}$ on $\Gamma_{h}$ by solving the variational equations (3.6) and (3.10), respectively and then put them into the solution representations (3.1) and (3.7) to get the numerical solutions $u_{h}$ and $U_{h}$ in $\Omega^{c}$, i.e.,

[TABLE]

5.2.1 Numerical examples for ESP

Set $\omega=1$, $\rho=1$, $\lambda=2$, $\mu=1$. Let the exact solution be

[TABLE]

Denote $\Gamma_{m}:=\{x=(x_{1},x_{2},x_{3})^{\top}\in{\mathbb{R}}^{3}:x_{1}=2\cos\theta,x_{2}=2,x_{3}=1.5\cos\theta,\theta\in[0,2\pi]\}$. Define the numerical error

[TABLE]

For simplicity, we use ’RP’ and ’IP’ to stand for ’real part’ and ’imaginary part’, respectively. The exact and numerical solutions on $\Gamma_{m}$ are plotted in Figure 3 for Obstacle I with $h=0.1005$. We observe that the numerical solutions are in a perfect agreement with the exact ones from the qualitative point of view. In Table 1, we present the numerical errors Error with respect to the meshsize $h$ which indicate the asymptotic convergence order $O(h^{2})$. These results verify the accuracy of the regularized formulation for hyper-singular BIO $W_{s}$.

Next, we consider the scattering of an incident plane wave $u^{in}$ taking the form

[TABLE]

by Obstacle II where $(\theta^{in},\phi^{in})$ is the incident direction. In this case, $f=-T(\partial,\nu)u^{in}$ on $\Gamma$. We choose $\theta^{in}=\pi/2$ and $\phi^{in}=0$. The real and imaginary parts of the numerical solution $u_{h}$ on four unit spheres surrounding the obstacle is presented in Figure 4.

5.2.2 Numerical examples for TESP

Choose $\omega=1$, $\rho=2$, $\lambda=1$, $\mu=1$, $\kappa=1$, $\eta=0.2$ and $\gamma=0.1$. The exact solution is set to be

[TABLE]

and $z=(0,1,0.3,0.2)^{\top}\in\Omega$. Define the numerical error

[TABLE]

We plot the exact and numerical solutions on $\Gamma_{m}$ in Figure 5 for Obstacle I with $h=0.1005$. The numerical solutions are in a perfect agreement with the exact ones from the qualitative point of view. In Table 2, we present the numerical errors $\widetilde{\mbox{Error}}$ with respect to the meshsize $h$ which also indicate the convergence. These results verify the accuracy of the regularized formulation for hyper-singular BIO $\widetilde{W}_{s}$.

Finally, we consider the scattering of an incident point source $U^{in}=({u^{in}}^{\top},p^{in})^{\top}$ taking the form

[TABLE]

by Obstacle II where $z$ is the location of point source. In this case, $F=-\widetilde{T}(\partial,\nu)U^{in}$ on $\Gamma$. We choose $z=(0,0,2)^{\top}$. The real and imaginary parts of the numerical solutions $U_{h}$ on four unit spheres surrounding the obstacle is presented in Figure 6.

Appendix A Proof of Theorem 4.1

We know from (4.9) that

[TABLE]

where

[TABLE]

Note that

[TABLE]

where

[TABLE]

We obtain from (4.3) that

[TABLE]

From (4.1) we can obtain that

[TABLE]

For $g_{3}(z)$, we know from (4.8) that

[TABLE]

Therefore, (1.2)-(1.4) yields

[TABLE]

where

[TABLE]

Note that for $i,j=1,2,3$,

[TABLE]

We conclude that

[TABLE]

We obtain from the Stokes formula (4.4) that

[TABLE]

On the other hand,

[TABLE]

Thus,

[TABLE]

Finally, since

[TABLE]

we have

[TABLE]

We complete the proof of (4.10) by a combination of (1.1) and (1.5)-(1.7).

Appendix B Proof of Lemma 4.2

For some matrix $A$ or vector $B$, we denote $(A)_{ij}$ and $(B)_{i}$ their Cartesian components, respectively. Let

[TABLE]

Then we have

[TABLE]

[TABLE]

and

[TABLE]

Therefore, from (4.2) and (2.1)-(2.3) we have

[TABLE]

Note that

[TABLE]

Hence,

[TABLE]

which completes the proof of (4.12). The proof of (4.13) follows in similar way, and we skip it here.

Appendix C Proof of Theorem 4.3

Following the same steps in A we can obtain that

[TABLE]

On the other hand, we have

[TABLE]

and

[TABLE]

Then we have

[TABLE]

Then (4.14) can be proved by combining (3.1) and (3.2).

Acknowledgments

The work of G. Bao is supported in part by a NSFC Innovative Group Fund (No.11621101), an Integrated Project of the Major Research Plan of NSFC (No. 91630309), and an NSFC A3 Project (No. 11421110002). The work of L. Xu is partially supported by a Key Project of the Major Research Plan of NSFC (No. 91630205), and a NSFC Grant (No. 11771068).