Time-harmonic acoustic scattering from a non-locally perturbed trapezoidal surface

TL;DR

This paper introduces a new radiation condition for acoustic scattering from a non-locally perturbed trapezoidal surface, proving uniqueness and existence of solutions, and develops a numerical mode matching method validated by experiments.

Contribution

It presents a novel radiation condition for non-local perturbations, proving theoretical results and creating an effective numerical method for acoustic scattering problems.

Findings

The new radiation condition is sharper than the Angular Spectrum Representation.

Existence and uniqueness of solutions are established under the new condition.

Numerical experiments confirm the validity and effectiveness of the proposed method.

Abstract

This paper is concerned with acoustic scattering from a sound-soft trapezoidal surface in two dimensions. The trapezoidal surface is supposed to consist of two horizontal half-lines pointing oppositely, and a single finite vertical line segment connecting their endpoints, which can be regarded as a non-local perturbation of a straight line. For incident plane waves, we enforce that the scattered wave, post-subtracting reflected plane waves by the two half lines of the scattering surface in certain two regions respectively, satisfies an integral form of Sommerfeld radiation condition at infinity. With this new radiation condition, we prove uniqueness and existence of weak solutions by a coupling scheme between finite element and integral equation methods. This consequently indicates that our new radiation condition is sharper than the Angular Spectrum Representation, and has generalized the radiation condition for scattering problems in a locally perturbed half-plane. Furthermore, we develop a numerical mode matching method based on this new radiation condition. A perfectly matched layer is setup to absorb outgoing waves at infinity. Since the medium composes of two horizontally uniform regions, we expand, in either uniform region, the scattered wave in terms of eigenmodes and match the mode expansions on the common interface between the two uniform regions, which in turn gives rise to numerical solutions to our problem. Numerical experiments are carried out to validate the new radiation condition and to show the performance of our numerical method.

Equations

$$\Gamma =\{(x_1,0)|x_1\le 0\}\cup V_h\cup \{(x_1,-h)|x_1\ge 0\}$$

$$(x_1,x_2)\in {\Omega }_{\Gamma }\to (x_1,x_2+s)\in {\Omega }_{\Gamma }\text{mbox}foralls>0,$$

$$\Delta u^{tot}+k^2{u}^{tot}=0$$

$$u^{tot}=0$$

$$\displaystyle\|u\|_{H^{1}_{\varrho}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}})}=\left[\int_{{\color[rgb]{0.000,0.000,0.000}\Omega_{a}}}\bigg{(}\big{|}(1+|x_{1}|^{2})^{\varrho/2}u\big{|}^{2}+\left|\nabla\left[(1+|x_{1}|^{2})^{\varrho/2}u\right]\right|^{2}\bigg{)}dx\right]^{1/2}\,.$$

$$\displaystyle||u||^{\prime}:=\left[\int_{{\color[rgb]{0.000,0.000,0.000}\Omega_{a}}}(1+|x_{1}|^{2})^{\varrho}\bigg{(}\big{|}u\big{|}^{2}+\big{|}\nabla u|^{2}\bigg{)}dx\right]^{1/2}\,,\quad u\in H^{1}_{\varrho}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}}).$$

$$u^{sc}(x)=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}\exp (i[(x_2-h)\sqrt{k^2-{\xi }^2}+x_1\xi ]){\hat{u}}_a^{sc}(\xi )d\xi ,x_1>a,$$

$${\hat{u}}_a^{sc}(\xi )=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}\exp (-i{x}_1\xi )u^{sc}(x_1,a)d{x}_1.$$

$$\underset{r\to \infty }{\lim }{\int }_{S_r}|{\partial }_{r}v-ikv{|}^{2}ds\to 0,\underset{r>0}{\sup }{\int }_{S_r}|v{|}^{2}ds<\infty ,$$

$${\Omega }^L:=\{x\in {\Omega }_{\Gamma }:x_2>x_1\tan \theta \},$$

$${\Omega }^R:=\{x\in {\Omega }_{\Gamma }:x_2<x_1\tan \theta -h/{\cos }^2\theta \},$$

$${\Omega }^M:=\{x\in {\Omega }_{\Gamma }:x_1\tan \theta -h/{\cos }^2\theta <x_2<x_1\tan \theta \}.$$

$$\mathcal{L}:=\{x=s(\cos \theta ,\sin \theta ):s>0\},$$

$${\mathcal{L}}^{\prime }:=\{x=s(\cos \theta ,\sin \theta )+(0,-h/{\cos }^2\theta ):s>0\}.$$

$$\displaystyle u^{tot}_{L}={\color[rgb]{0.000,0.000,0.000}u^{in}+u^{re}}=e^{i(\alpha x_{1}-\beta x_{2})}-e^{i(\alpha x_{1}+\beta x_{2})},$$

$$\displaystyle u^{tot}_{R}={\color[rgb]{0.000,0.000,0.000}u^{in}+u^{re}_{h}}=e^{i(\alpha x_{1}-\beta x_{2})}-c_{h}\,e^{i(\alpha x_{1}+\beta x_{2})},$$

$$\displaystyle v:=\left\{\begin{array}[]{lll}u^{tot}-u^{tot}_{L}&&\mbox{in}\quad\Omega^{L},\\ u^{tot}-u^{tot}_{R}&&\mbox{in}\quad\Omega^{M}\cup\Omega^{R};\end{array}\right.$$

$$\displaystyle v^{\prime}:=\left\{\begin{array}[]{lll}u^{tot}-u^{tot}_{L}&&\mbox{in}\quad\Omega^{L}\cup\Omega^{M},\\ u^{tot}-u^{tot}_{R}&&\mbox{in}\quad\Omega^{R}.\end{array}\right.$$

$$\displaystyle\begin{split}v^{+}-v^{-}=(1-c_{h})e^{i(\alpha x_{1}+\beta x_{2})}\quad&&\mbox{on}\quad{\mathcal{L}},\\ (v^{\prime})^{+}-(v^{\prime})^{-}=(1-c_{h})e^{i(\alpha x_{1}+\beta x_{2})}\quad&&\mbox{on}\quad{\mathcal{L}}^{\prime}.\end{split}$$

$$\displaystyle\begin{split}\partial_{\nu}^{+}v-\partial_{\nu}^{-}v=0\quad&&\mbox{on}\quad{\mathcal{L}},\\ \partial_{\nu}^{+}(v^{\prime})-\partial_{\nu}^{-}(v^{\prime})=0\quad&&\mbox{on}\quad{\mathcal{L}}^{\prime}.\end{split}$$

$${\int }_{S_r}|{\partial }_{r}v-ikv{|}^{2}ds$$

$${\int }_{S_r}|v{|}^{2}ds$$

$${\int }_{S_r^M}|{\partial }_{r}w-ikw{|}^{2}ds\to 0,{\int }_{S_r^M}|w{|}^{2}ds=O(1)\text{mbox}asr\to 0.$$

$${\int }_{S_r^M}|{\partial }_{r}w-ikw{|}^{2}ds$$

$${\int }_{S_r^M}|w{|}^{2}ds=(c_h-1{)}^{2}r{\int }_{{\theta }^{\ast }(r)}^{\theta }d\phi =O(1)$$

$$u^{tot}=u_L^{tot}+v\text{mbox}in{\Omega }^L,u^{tot}=u_R^{tot}+v\text{mbox}in{\Omega }^M\cup {\Omega }^R,$$

$$\underset{r\to \infty }{\lim }{\int }_{S_r\cap \{x_2\ge a\}}|{\partial }_r\Phi -ik\Phi {|}^{2}ds\to 0,\underset{r>0}{\sup }{\int }_{S_r\cap \{x_2\ge a\}}|\Phi {|}^{2}ds<\infty ,r=|x|.$$

$$\sqrt{r}({\partial }_r\Phi -ik\Phi )\to 0\text{mbox}asr\to \infty ,x_2\ge a.$$

$$v(x)$$

$${\int }_{S_{R_1}}v(y){\partial }_{\nu }\Phi (x;y)-{\partial }_{\nu }v(y)\Phi (x;y)ds(y)$$

$$v(x)={\int }_{S_R}\left[v(y){\partial }_{\nu }\Phi (x;y)-{\partial }_{\nu }v(y)\Phi (x;y)\right]ds(y)+f(x)$$

$$\displaystyle f(x):={\color[rgb]{0.000,0.000,0.000}(1-c_{h})}\lim_{R_{1}\rightarrow\infty}\int_{{\mathcal{L}}_{R_{1}}}e^{i(\alpha y_{1}+\beta y_{2})}\partial_{\nu}\Phi(x;y)\;ds(y).$$

Full text

Furthermore, we develop a numerical mode matching method based on this new radiation condition. A perfectly matched layer is setup to absorb outgoing waves at infinity. Since the medium composes of two horizontally uniform regions, we expand, in either uniform region, the scattered wave in terms of eigenmodes and match the mode expansions on the common interface between the two uniform regions, which in turn gives rise to numerical solutions to our problem. Numerical experiments are carried out to validate the new radiation condition and to show the performance of our numerical method.

Keywords: Helmholtz equation, trapezoidal surface, Sommerfeld radiation condition, half plane, non-local perturbation, mode matching method, perfectly matched layer.

1 Introduction

Wave scattering in a layered medium and in a half plane has numerous applications in scientific and engineering areas [14]. For the purpose of proving well-posedness of the scattering model, understanding the physical radiation behavior of the wave field at infinity is critical. A rigorous description of the asymptotics not only helps to design a proper radiation condition for the wave field at infinity, but also helps to truncate the unbounded domain with an accurate boundary condition for numerically solving the problem.

Certainly, the radiation condition is structure-related. For instances, if a bounded obstacle is embedded into a homogeneous background medium, the scattered wave is purely outgoing at infinity and satisfies the classic Sommerfeld radiation condition (SRC) [17]. When the structure is filled in by a two-layered medium with a locally perturbed planar surface, the perturbed wave field due to the local perturbation, is outgoing at infinity and satisfies the SRC [36, 13, 2, 28]; see also [1, 40, 24, 39] for studies on impenetrable locally perturbed surfaces. However, in the case of a globally perturbed rough surface, by which we mean a non-local perturbation of a planar surface such that the surface lies within a finite distance of the original plane, one in general cannot explicitly extract an outgoing wave from the scattered wave to meet the SRC. The Angular Spectrum Representation (ASR) radiation condition (see [18, 9, 10]) or the Upward Propagating Radiation Condition (UPRC) (see [41, 42, 12]) can be viewed as a rigorous formulation of a radiation condition to show the well-posedness of the problem in both two and three dimensions. The radiation condition relies also on the type of incoming waves in rough surface scattering problems. The authors in [10] proved that the scattered field incited by a plane wave decays slower than that for a point source wave in the horizontal direction. It was recently proved in [22] that the scattered field due to a point source wave fulfills the Sommerfeld radiation condition in a half plane, which however does not hold true for plane wave incidence.

In this paper, we shall study a special class of two-dimensional (2D) rough surface scattering problems, where the globally perturbed surface is assumed to consist of two horizontal half lines pointing oppositely and one single vertical line segment connecting their endpoints; we shall propose a novel SRC-type radiation condition that is sharper than ASR. Relying on this new radiation condition, we shall prove existence and uniqueness of weak solutions and design a numerical method for the scattering problem.

Let $\Gamma$ be a sound-soft rough surface in two dimensions and suppose that the region above $\Gamma$, which we denote by $\Omega_{\Gamma}$, is occupied by a homogeneous and isotropic medium. Consider a time-harmonic plane wave $u^{in}=e^{ikx\cdot d}$ incident onto the rough surface $\Gamma\subset{\mathbb{R}}^{2}$ from above. Here, $k>0$ denotes the wave number, $d=(\cos\theta,-\sin\theta)$ stands for the incident direction and $\theta\in(0,\pi)$ is the incident angle with the positive $x_{1}$-axis. In this paper, we suppose that $\Gamma$ consists of three parts (see Fig. 1):

[TABLE]

where $V_{h}=\{(0,x_{2})|-h\leq x_{2}\leq 0\}$ is a finite vertical line segment with the height $h>0$.

Such kind of trapezoidal surfaces is a non-local perturbation of the straight line $\{x_{2}=0\}$, which can be regarded as a special kind of globally perturbed rough surfaces. Since the region $\Omega_{\Gamma}$ fulfills the geometrical condition

[TABLE]

by [10] there exists a unique solution {\color[rgb]{0.000,0.000,0.000}u^{tol}}\in H_{\varrho}^{1}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}}) for any $a>0$ and $\varrho\in(-1,-1/2)$ such that

[TABLE]

Here, {\color[rgb]{0.000,0.000,0.000}\Omega_{a}}=\{x\in\Omega_{\Gamma},x_{2}<a\} is the unbounded strip between $\Gamma$ and $x_{2}=a$ (see Figure 1), and H^{1}_{\varrho}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}}) stands for weighted Sobolev space defined as

[TABLE]

One can also employ the following equivalent norm to ||\cdot||_{H^{1}_{\varrho}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}})}:

[TABLE]

Further, the scattered field $u^{sc}=u^{tot}-u^{in}$ in $\Omega_{\Gamma}$ is required to fulfill the Upward Angular Spectrum Representation (UASR), which can be written as

[TABLE]

where $\hat{u}^{sc}_{a}(\xi)$ denotes the Fourier transform of $u^{sc}(x_{1},a)$ in $x_{1}$, given by

[TABLE]

In this equation $\sqrt{k^{2}-\xi^{2}}=i\sqrt{\xi^{2}-k^{2}}$ when $\xi^{2}>k^{2}$. The representation of $u^{sc}$ in the integral (1.6) can be interpreted as a formal radiation condition in the physics and engineering literature on rough surface scattering (see e.g. [9, 18]).

If $\Gamma$ is a local perturbation of the original straight line $\{x_{2}=0\}$ , it was proved in a series of papers (see [1, 2, 24, 28, 39, 40]) that $u^{sc}$ can be decomposed into two parts: $u^{sc}=u^{re}+v$ where $u^{re}$ is the scattered field corresponding to unperturbed surface $\{x_{2}=0\}$ and $v\in H^{1}(\Omega_{\Gamma})$ is caused by the local perturbation satisfying the half-plane Sommerfeld radiation condition

[TABLE]

where $S_{r}:=\{x\in\Omega_{\Gamma}:|x|=r\}$. Physically, $u^{re}=-e^{ikx\cdot d^{\prime}}$ with $d^{\prime}:=(d_{1},-d_{2})$ is uniquely determined by Snell’s law and (1.7) indicates that $v$ approximately becomes purely outgoing at infinity. We note that both $u^{re}$ and $v$ fulfill the ASR (1.6), so that $u^{sc}$ also satisfies the ASR. The decomposition of $u^{sc}$ into the sum of $u^{re}$ and $v$ provides a deep insight into the wave phenomenon in a locally perturbed half plane, which however cannot apply to our scattering problem in a non-locally perturbed half-plane. In this paper, we shall propose a new radiation condition, which not only satisfies the ASR but also allows us to generalize the above decomposition for the non-locally perturbed surface $\Gamma$ under consideration.

Since $\Gamma$ in (1.1) is no longer a local perturbation of $\{x_{2}=0\}$, subtracting $u^{re}$, defined in the last paragraph, from $u^{sc}$ no longer yields an outgoing wave at infinity; clearly, $u^{sc}$ contains a reflected plane wave $u^{re}_{h}=-e^{ikx^{\prime}_{h}\cdot d}$ with $x_{h}^{\prime}:=(x_{1},-2h-x_{2})$, due to the half line $\{x_{2}=-h\}$, which is in general different from $u^{re}$ for $h\neq 0$. It turns out that $u^{re}_{h}$ and $u^{re}$ reside in two non-overlapping regions, respectively, which occupy the whole medium above $\Gamma$ but are separated by a straight line ${\mathcal{L}}$ parallel to $d^{\prime}$, the propagating direction of $u^{re}$. Consequently, we enforce that subtracting $u^{re}$ in the region left to ${\mathcal{L}}$ and then $u^{re}_{h}$ in the region right to ${\mathcal{L}}$ from $u^{sc}$, yields an outgoing wave at infinity satisfying (1.7). Based on this new radiation condition, we prove uniqueness and existence of weak solutions by a coupling scheme between finite element and integral equation methods.

To numerically compute $u^{sc}$, we place a perfectly matched layer (PML) [3] to absorb the outgoing waves extracted from $u^{sc}$. Since $u^{sc}$ differs from outgoing waves by only reflected plane waves propagating parallel to $d^{\prime}$, we impose a homogeneous Robin-type boundary condition to eliminate $u^{re}$ and $u^{re}_{h}$ on the boundary of the PML by following [32]. As the medium structure can be decomposed into two $x_{1}$-uniform segment, the numerical mode matching (NMM) method of [32] will be adapted to our scattering problem. More precisely, we expand $u^{sc}$ in terms of eigenmodes in each uniform segment by separating variables in the governing equations, and then match the two mode expansions on the common interface separating the two regions. This in turn gives rise to a linear system of the unknown Fourier coefficients in the expansions. Solving the linear system yields a numerical solution to our scattering problem. Numerical experiments shall be carried out to show the performance of our numerical methods based on the new radiation condition.

The NMM method [14, 32], a.k.a mode expansion method or modal methods [6, 29, 37], and its numerous numerical invariants [16, 20, 21, 25, 27, 30, 33, 35, 38, 19, 4, 5] are applicable when the structure can be divided into a number of segments, where the medium becomes uniform along one spatial variable. The classical mode matching method solves the eigenmodes analytically while the NMM methods solve the eigenmodes by numerical methods, and they are easier to implement and applicable to more general structures. The mode matching method and its variants have the advantage of avoiding discretizing one spatial variable. They are widely used in engineering applications, since many designed structures are indeed piecewise uniform.

For numerical simulations of waves, the perfectly matched layer (PML) [3, 15] is an important technique for truncating unbounded domains. It is widely used with standard numerical methods, such as finite element methods [34] and spectral methods, etc., that discretize the whole computational domain. For structures that are piecewise homogeneous, the boundary integral equation (BIE) methods [8, 7, 26] are popular since they can automatically take care of radiation conditions at infinity while discretizing only interfaces of the structure. For scattering problems in layered media, PML can also be incorporated with BIE methods to efficiently truncate interfaces that extend to infinity [31].

The remaining of this paper is organized as follows. In the subsequent section 2, we enforce a new radiation condition on the scattered field and prove uniqueness and existence of solutions to our globally perturbed scattering problem. In section 3, we develop a NMM method to solve the scattering problem. Section 4 is devoted to numerical examples.

2 Well-posedness

This section is devoted to existence and uniqueness of the 2D wave scattering problem in an upper-half space with a trapezoidal sound-soft boundary $\Gamma=\{(x_{1},0)|x_{1}\leq 0\}\cup\{(0,x_{2})|-h\leq x_{2}\leq 0\}\cup\{(x_{1},-h)|x_{1}\geq 0\}$ for $h>0$. A new radiation condition will be proposed in subsection 2.1 and our main result, Theorem 2.4, will be proved in subsection 2.2.

2.1 Radiation condition

Without loss of generality, we suppose that the incident angle $\theta\in(0,\pi/2)$. This means that the incoming wave is incident onto $\Gamma$ from the left hand side of the upper half plane. We divide the region $\Omega_{\Gamma}$ into three parts $\Omega_{\Gamma}=\Omega^{L}\cup\Omega^{M}\cup\Omega^{R}$ with (see Figure 2)

[TABLE]

These domains are separated by the following two rays

[TABLE]

Denote by $S_{r}^{\eta}$ ($\eta=L,M,R$) the restriction of the circular curve $S_{r}$ to $\Omega^{\eta}$, that is, $S_{r}^{\eta}=\{x:x\in S_{r}\cap\Omega^{\eta}\}$. Let $u^{tot}_{L}$, $u^{tot}_{R}$ be the uniquely determined total fields incited by the plane wave $u^{in}$ incident on the sound-soft straight lines $\{x_{2}=0\}$ and $\{x_{2}=-h\}$, respectively. We refer to [9, 10, 41, 42] for uniqueness and existence of the solution in Hölder continuous spaces or in weighted Sobolev spaces using integral equation or variational methods. Mathematically, they are given explicitly by

[TABLE]

where \alpha:=k{\color[rgb]{0.000,0.000,0.000}\cos\theta},\beta:=k{\color[rgb]{0.000,0.000,0.000}\sin\theta}, $c_{h}:=e^{2i\beta h}$. Denote by $u^{tot}$ the total field to our scattering problem. To introduce our radiation condition, we need to define two functions

[TABLE]

Obviously, $v$ coincides with $v^{\prime}$ in $\Omega^{L}\cup\Omega^{R}$ and $v$ differs from $v^{\prime}$ only over the region $\Omega^{M}$. We remark that, since $u^{tot}\in H^{1}_{loc}(\Omega_{\Gamma})$, the function $v$ is discontinuous on {\color[rgb]{0.000,0.000,0.000}{\cal L}}, whereas $v^{\prime}$ is discontinuous on {\color[rgb]{0.000,0.000,0.000}{\cal L}}^{\prime}. More precisely, it holds that

[TABLE]

Here, the notation $(\cdot)^{\pm}$ denote respectively the limits taking from left and right hand sides. However, the normal derivatives of $v$ (resp. $v^{\prime}$) are continuous when getting across {\color[rgb]{0.000,0.000,0.000}{\cal L}} (resp. {\color[rgb]{0.000,0.000,0.000}{\cal L}}^{\prime}), that is,

[TABLE]

In the following lemma we show that the half-plane Sommerfeld radiation conditions of $v$ and $v^{\prime}$ are equivalent.

Lemma 2.1

The function $v$ satisfies the Sommerfeld radiation condition (1.7) if and only if $v^{\prime}$ satisfies (1.7).

Proof. We first note that the Sommerfeld radiation condition of $v$ are understood as

[TABLE]

respectively, as $r=|x|\rightarrow\infty$. Set $w=v-v^{\prime}$. Then w={\color[rgb]{0.000,0.000,0.000}u^{re}-u_{h}^{re}}=(c_{h}-1)e^{i(\alpha x_{1}+\beta x_{2})} in $\Omega^{M}$ and $w=0$ in $\Omega^{L}\cup\Omega^{R}$. Hence, we only need to prove that

[TABLE]

Let $(r,\varphi)$ be the polar coordinate of $x$, and let $(r,\theta^{*}(r))$ be the polar coordinate of the intersection point of $S_{r}$ and {\color[rgb]{0.000,0.000,0.000}\cal L}^{\prime}. Using the fact that $w(x)=w(r,\varphi)=(c_{h}-1)e^{ikr\cos(\theta-\varphi)}$ together with the mean value theorem, it is easy to see

[TABLE]

as $r\rightarrow\infty$, since $|\theta^{*}(r)-\theta|=O(1/r)$. Analogously, the second relation in (2.16) follows from the inequality

[TABLE]

as $r\rightarrow\infty$. The proof of the lemma is complete. $\Box$

By Lemma 2.1, our new radiation condition in $\Omega_{\Gamma}$ is defined as follows.

Definition 2.2

The scattered field $u^{tot}-u^{in}$ is said to be outgoing if the function $v$ or $v^{\prime}$ satisfies the half-plane Sommerfeld radiation condition (1.7).

Below we present several remarks concerning this new radiation solution.

Remark 2.3

• (i)

Obviously, the definition of our outgoing radiation condition depends on the incident angle $\theta$ and the height $h$ of the scattering surface $\Gamma$ in the vertical direction. If $h=0$ (i.e., $\Gamma$ is a local perturbation of the original straight line $\{x_{1}=0\}$), then we see $u_{L}^{tot}=u_{R}^{tot}$ and thus the proposed radiation condition could be reduced to the usual condition for scattering problems in a locally-perturbed half-plane (see e.g., [1, 2, 24, 28, 39, 40]). Moreover, we remark that $u-u^{in}$ still satisfies the Angular Spectrum Representation (1.6) for general rough surface problems.

• (ii)

*In the definition of $v$ (see (2.10)), the ray ${\mathcal{L}}$ can be replaced by another ray lying in $\Omega_{\Gamma}$ which starts from any point on $\Gamma$ with the direction $(\cos\theta,\sin\theta)$. By the proof of Lemma 2.1, the Sommerfeld radiation condition of $v$ is also equivalent to that of the new function defined in the modified domain. Such substitution also applies to $v^{\prime}$ and the ray ${\mathcal{L}}^{\prime}$. *

Suppose that $u-u^{in}$ is an outgoing radiation solution. By Definition 2.2, the total field $u^{tol}$ can be decomposed into

[TABLE]

where $u_{L}^{tot}$ and $u_{R}^{tot}$ are defined in (2.28) and the function $v$ fulfills the Sommerfeld radiation condition (1.7).

The main results of this section is stated in the following theorem. Its proof will be carried out in the subsequent section.

Theorem 2.4

Assume that the scattering interface $\Gamma$ is given by (1.1) and $u^{in}=e^{ikx\cdot d}$ is a plane wave. Then the scattering problem (1.3)-(1.4) admits a unique solution of the form (2.17). Moreover, it is the unique solution in the weighted Sobolev space H_{\varrho}^{1}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}}) for any $a>0$ and $\varrho\in(-1,-1/2)$.

2.2 Proof of Theorem 2.4

Let $G(x,z)=\frac{i}{4}H_{0}^{(1)}(k|x-y|)$ be the free-space Green’s function to the Helmholtz equation, where $H_{0}^{(1)}$ is the Hankel function of the first kind of order zero. Suppose that $u^{in}(x,y)=G(x,y)$ is an incoming point source wave emitting from the source position located by $y\in\Omega_{\Gamma}$. By [10], there admits a unique scattered field belonging to the space H_{\varrho}^{1}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}}) for any $a>0$ and $\varrho\in(-1,0)$, which satisfies the UASR (1.6) in $x_{2}>a$. Let $\Phi(x,y)$ ($y\in\Omega_{\Gamma}$) be the unique total field caused by the incoming point source wave $u^{in}(x,y)$. Obviously, $\Phi(x,y)$ can be regarded as the Green’s function to our scattering problem. The proof of Theorem 2.4 relies essentially on the following proposition.

Proposition 2.5

The Green’s function $\Phi(x,y)(y\in\Omega_{\Gamma})$ fulfills the half-plane Sommerfeld radiation condition (1.7).

When $\Gamma$ is the graph of a $C^{1,1}$-smooth function, the assertion of Proposition 2.5 for a compactly supported source term is already contained in [11, Theorem 5.1] but without a detailed proof. It was further proved in [22, Lemma 2.2] for general rough surface scattering problems that \Phi(\cdot,y)\in H_{\varrho}^{1}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}}\cap\{x\in{\mathbb{R}}^{2}:|x_{1}|>R\}) for any $a>\max_{x\in\Gamma}\{x_{2}\}$, $|\varrho|<1$ and $R>|y_{1}|$, and that

[TABLE]

Moreover, the above radiation conditions are proved to be equivalent to the classical Sommerfeld radiation condition that

[TABLE]

In our case of the non-locally perturbed surface $\Gamma$ given by (1.1), we may choose $a\geq 0$ because the solution is continuous up to the surface {$x_{2}=0$}. In the remaining part of this subsection we shall verify that the total field of our scattering problem for plane wave incidence can be decomposed into the form of (2.17). We will follow the coupling scheme of [2, 23, 28] between the finite element and boundary integral equation methods, but modified to be applicable to our case in a non-locally perturbed half plane.

Set $\Omega_{R}^{+}=\{x:|x|<R\}\cap\Omega_{\Gamma}$. Choose $R>h$ such that $k^{2}$ is not the eigenvalue of $-\Delta$ in $\Omega_{R}^{+}$; see Figure 3.

Choose $R_{1}>R$. Using Green’s formula, we may represent $v$ in $\Omega^{L}\cap\{x:R<|x|<R_{1}\}$ as the integral representation

[TABLE]

with ${\mathcal{L}}_{R_{1}}:=\{x\in{\mathcal{L}}:R<|x|<R_{1}\}$ for $R_{1}>R$. Here we have assumed that the normal vector $\nu$ on ${\mathcal{L}}$ is directed into the left hand side and that on $S_{R}$ into the exterior $|x|>R$; see Figure 3. Analogously, for $x\in(\Omega^{R}\cup\Omega^{M})\cap\{x:R<|x|<R_{1}\}$ we get

[TABLE]

Adding (2.18) and (2.19) together and making use of the jump conditions of $v$ on ${\mathcal{L}}$ (see (2.14) and (2.15)), we obtain

[TABLE]

for $x\in\Omega^{L}\cap\{x:R<|x|<R_{1}\}$. In view of the Sommerfeld radiations of $v$ and $G$, we may find that

[TABLE]

as $R_{1}\rightarrow\infty$. Hence,

[TABLE]

for all x\in\Omega_{\Gamma}{\color[rgb]{0.000,0.000,0.000}\backslash\overline{\Omega_{R}^{+}}}, where

[TABLE]

The proof of the existence of the limit on the right hand side of (2.21) will be given in the Appendix. Note that $f$ vanishes identically if $h=0$ or $u^{in}=0$, and that the integral on $S_{R}$ appearing in (2.20) is understood as the sum of those integrals over $S_{R}^{L}$ and $S_{R}^{M}\cup S_{R}^{R}$. Since the Green’s function $\Phi$ is weakly singular, the jump relation for double layer potentials gives (cf. (2.14))

[TABLE]

For notational simplicity, we write $S_{1}=S_{R}^{L}$, S_{2}=S_{R}^{M}\cup S_{R}^{{\color[rgb]{0.000,0.000,0.000}R}} and define

[TABLE]

Taking the limit $x\rightarrow S_{R}$, we get

[TABLE]

Here $I$ is the identify operator, $\mathcal{D}$ and $\mathcal{S}$ are defined by

[TABLE]

for $i,j=1,2$. We remark that the jump relations for $D_{ij}$ and $S_{ij}$ remain valid, since $\Phi(\cdot;y)-G(\cdot;y)$ is of $C^{\infty}$-smoothness. On the other hand, using integration by part we may find

[TABLE]

for all $\varphi\in H^{1}(\Omega_{R}^{+})$ such that $\varphi=0$ on $\Gamma\cap\{x:|x|<R\}$. This implies that

[TABLE]

where the function $w$ is defined by (cf. (2.10))

[TABLE]

Introduce the variational space $X=X_{0}\times X_{1}$, where

[TABLE]

Combining (2.23) and (2.22) gives the variational formulation for the unknown solution pair $(u^{tot},p)\in X$ with $p=(p_{1},p_{2})^{\top}\in X_{1}$ as follows

[TABLE]

for all $(\varphi,\chi)\in X$ with $\chi=(\chi_{1},\chi_{2})^{\top}$, where $\hat{f}:=(f|_{S_{1}},f|_{S_{2}})^{\top}$ and

[TABLE]

Note that $\hat{u}^{tot}:=({u}^{tot}|_{S_{1}},{u}^{tot}|_{S_{2}})^{\top}$. Here we have used the notation

[TABLE]

In comparison with the variational formulation for local perturbation scattering problems (see e.g., [22, 28]), we have an additional term $\hat{f}$ appearing on the right hand side of (2.27), due to the fact that $h>0$. In addition, the integral over $S_{R}$ is split into the sum of corresponding integrals over $S_{1}$ and $S_{2}$, because of the unknown functions $p_{1}$ and $p_{2}$. The operator $A$ takes the same form as the case of a locally perturbed half plane. Hence, using mapping properties of $D_{ij}$, $S_{ij}$ and arguing analogously to [22], one can prove that $A:X\rightarrow X^{*}$ is a Fredholm operator with index zero. We omit the details, since the proof of [22, 28] carries over to our case easily.

To prove uniqueness, we assume that $u^{in}=0$. This implies that $w=0$ and $f=0$; recall (2.21) and (2.26) for the definition. Hence, the right hand side of the variational formulation (2.27) vanishes and $u^{tot}=v$ fulfills the radiation condition (1.7). Since $k^{2}$ is not the Dirichlet eigenvalue of $-\Delta$, one can extend the solution $u^{tol}$ of the homogeneous boundary value problem A\left(({\color[rgb]{0.000,0.000,0.000}u^{tol}},p),(\varphi,\chi)\right)=0 for all $(\varphi,\chi)\in X$ from $\Omega_{R}^{+}$ to the whole half plane $\Omega_{\Gamma}$, which also satisfies the Sommerfeld radiation condition (1.7). The extended solution also fulfills the Angular Spectrum Representation. Hence, by uniqueness to rough surface scattering problem under the geometrical condition (1.2) (see [10]), we get $u^{tot}=0$, which proves the uniqueness of solutions to the problem (2.27). Existence follows straightforwardly from Fredholm alternative theory. This finishes the proof of Theorem 2.4. $\Box$

Once $u^{tot}\in H^{1}(\Omega_{R})$ (and thus $v$) is obtained from (2.27), the solution $u^{tot}$ can be extended from $\Omega_{R}^{+}$ to the region $|x|>R$ via $u^{tot}=w+v$, where $v$ is expressed by (2.20) in terms of the trace of $v$ on $S_{R}$ and the function $f$.

Remark 2.6

*Theorem 2.4 can be readily carried over to the following cases: *

• (i)

The incident angle $\theta\in(\pi/2,\pi)$, or the scattering interface is a local perturbation of $\Gamma$ defined by (1.1) (that is, the interface coincides with $\Gamma$ in the exterior of a compact set). The non-local surface shown in Figure 4 can be analogously treated as well. Note that our approach applies to a half-plane which satisfies the geometrical assumption (1.2). As indicated by Remark 2.3 (ii), we may choose a ray starting from any point on $\Gamma$ with the direction $(\cos\theta,\sin\theta)$ as ${\cal L}$ to define $\Omega_{L}$ as the domain on the left of ${\cal L}$. The ray ${\cal L}^{\prime}$ can be chosen as any ray parallel to and on the right hand side of ${\cal L}$. Then $\Omega_{R}$ is defined as the domain on the right of ${\cal L}^{\prime}$ and $\Omega_{M}$ as the middle domain between ${\cal L}$ and ${\cal L}^{\prime}$. One easily checks that Lemma 2.1 is still valid.

• (ii)

An inhomogeneous medium with compact contrast function is embedded into $\Omega_{\Gamma}$. In this case, the wave equation for the total field becomes $(\Delta+k^{2}q(x))u^{tot}=0$ in $\Omega_{\Gamma}$, where $q\in L^{\infty}(\Omega_{\Gamma})$ and $D:=\mbox{Supp}(1-q)$ is a compact set of $\Omega_{\Gamma}$. The local perturbation is understood as the scattering effect due to the compactly supported inhomogeneity $q$.

• (iii)

The incoming wave is a point source wave emitted from some source position located in $\Omega_{\Gamma}$, that is, $u^{in}(x)=\Phi(x,z)$ for some $z\in\Omega_{\Gamma}$. Then the total field can be decomposed into the form (2.17) with

[TABLE]

where $z^{\prime}=(z_{1},-z_{2})$, $z^{\prime}_{h}=(z_{1},-2h-z_{2})$ for $z=(z_{1},z_{2})$.

Remark 2.7

The fundamental differences between the arguments in [22] and the proof of Theorem 2.4 are summarised as follows. In [22], a similar coupling scheme was employed to establish well-posedness of time-harmonic acoustic scattering from a locally perturbed sound-soft periodic surface. The mathematical analysis there was mostly placed upon the justification of the Sommerfeld radiation condition of the Green’s function (that is, Proposition 2.5). However, in this paper we consider an acoustic scattering problem in a globally perturbed half plane. Compared to the case of local perturbation, essential difficulties for trapezoidal surfaces arise from the additional term $\hat{f}$ on the right hand side of the new variational formulation (2.27). We have to prove the convergence of the limit in the definition of the function $f$ (see (2.21)) , which requires ingenious analysis to estimate the asymptotic behavior of the Green function $\Phi(x,z)$ as $|x|\rightarrow\infty$ uniformly in all $z\in S_{R}$; see the lengthy arguments in the Appendix for the details.

The variational formulation established in this section is helpful to establish well-posedness of our scattering problem. However, it is hard to implement the resulting numerical scheme, due to the heavy computational cost on the background Green’s function. Instead, a mode matching method will be adopted in the subsequent section to get numerical solutions, where the decomposition form (2.17) will be used to truncate the unbounded domain with an accurate boundary condition.

3 A numerical mode matching method

Without the knowledge of a precise radiation behavior of the total wavefield $u^{tot}$ infinity, one cannot apply existing truncation techniques such as PML method, absorbing boundary condition (ABC) method, etc., to truncate $\Omega_{\Gamma}$ as we don’t know at all what boundary conditions should be imposed after the truncation, let alone developing further numerical methods to compute $u^{tot}$! In this section, we will propose a numerical mode matching (NMM) method to compute $u^{tot}$ utilizing the newly proposed radiation condition.

We remark that there are some major differences between the current work and [32]. In [32], an NMM method has been developed for the scattering problem in a two-layer medium with a stratified inhomogeneity, where an outgoing wavefield is much easier to extract in terms of directly subtracting the background reference wavefield from the total wavefield so that the aformentioned truncation techniques can be easily applied. Moreover, a hybrid Robin-Dirichlet boundary condition was proposed therein on the PML boundary to make the mode expansion procedure applicable, but this in fact is unnecessary in some standard methods like FEM methods since one can simply put zero Dirichlet boundary condition on the whole PML boundary to terminate the outgoing wavefield. Unlike [32], no uniform background reference wavefield is available for us to extract an outgoing wavefield for the scattering problem under consideration. Whichever one chooses, $u^{tot}_{L}$ or $u^{tot}_{R}$, as the background reference wavefield, (2.17) tells that neither $u^{tot}-u^{tot}_{L}$ nor $u^{tot}-u^{tot}_{R}$ is outgoing since the two difference wavefields contain plane waves parallel to $e^{i(\alpha x_{1}+\beta x_{2})}$ in $\Omega_{M}\cup\Omega_{R}$ and $\Omega_{L}$, respectively. Nevertheless, one can use PMLs and then zero Dirichlet boundary condition to directly terminate the outgoing wavefield $v$, but one definitely meets a “pseudointerface” dependent on the incident angle, i.e., ray ${\cal L}$, across which $v$ is discontinuous. To avoid such a moving pseudointerface, we have two approaches: (1). Directly compute the partially outgoing wavefield $u^{tot}-u^{tot}_{L}$ (or $u^{tot}-u^{tot}_{R}$) in the whole computational domain; we must impose a hybrid Robin-Dirichlet boundary condition on the PML boundary to eliminate the reflection of the plane wave $e^{i(\alpha x_{1}+\beta x_{2})}$; (2). Setup a fixed pseudointerface, e.g., $x_{1}=0$, and then compute $u^{tot}-u^{tot}_{L}$ on the left of the pseudointerface and $u^{tot}-u^{tot}_{R}$ on the right; like (1), we still need to impose a hybrid Robin-Dirichlet boundary condition on the PML boundary since one of the two still is partially outgoing. Consequently, Robin-Dirichlet boundary condition is necessary in any numerical methods unless one could tolerate the movement of pseudointerface as incident angle varies. Here, we propose to use the second approach and will develop an NMM method to compute $u^{tot}$.

For simplicity, we assume that $\Gamma$ is defined in (1.1) and there is no inhomogeneity above $\Gamma$; this NMM method is applicable as well for $\Gamma$ locally perturbed with multiple vertical and horizontal segments and with additional rectangular inhomogeneities above; see Remark 2.6 (i) for details.

As shown in Figure 1, the vertical $x_{2}$-axis splits $\Omega_{\Gamma}$ into two $x_{1}$-invariant mediums

[TABLE]

Introduce the following function

[TABLE]

notice that $u^{sc}$ in general is not outgoing. Then, $u^{sc}$ satisfies

[TABLE]

In $\Omega_{1}$, applying the method of separation of variables, we insert $u^{sc}(x_{1},x_{2})=\phi(x_{2})\psi(x_{1})$ into (3.32) and obtain the following eigenvalue equations for $\phi$

[TABLE]

and the associated equation for $\psi$

[TABLE]

Next, along $x_{2}$-axis, we use a complex-coordinate transformation

[TABLE]

where $\sigma_{2}(t)=0$ for $t<L$ and $\sigma_{2}(t)>0$ for $t\geq L$ and the half-plane $x_{2}\geq L$ with a nonzero absorbing function $\sigma_{2}$ is called the perfectly matched layer (PML), which is capable of absorbing outgoing waves quite efficiently [3, 15]. According to Eqs. (2.10) and (3.31), $u^{sc}$ and the outgoing wave $v$ in $\Omega_{1}$ differ by a multiple of the reflected plane wave $e^{i\alpha x_{1}+i\beta x_{2}}$ in $\Omega_{1}/\Omega_{L}$ only when $\Omega_{L}\subset\Omega_{1}$. If $\beta$ is close to [math], this plane wave can propagate nearly parallel to the PML entrance $x_{2}=L$, which numerically causes inefficiency of the PML absorption. To resolve this issue, our previous work [32] suggests to eliminate the plane wave by the following Robin-type relation,

[TABLE]

The Sommerfeld radiation condition (1.7) for $v$ implies that $v$ is approximately an outgoing wave $e^{ikr}$ at infinity so that it decays rapidly in the PML as $x_{2}\rightarrow\infty$. Consequently, since $\partial_{x_{2}}v-i\beta v$ also produces an outgoing wave at infinity,

[TABLE]

decays as well, implying

[TABLE]

approaches [math] for $x_{2}\rightarrow\infty$. Thus, by setting $\hat{\phi}(x_{2})=\phi(\hat{x}_{2})$ and by terminating the PML layer at $x_{2}=L+d$ for the PML thickness $d>0$, we get from (3.35-3.36) that

[TABLE]

with the following boundary conditions

[TABLE]

Employing the Chebyshev collocation method in [38] to solve the above eigenvalue problems Eqs. (3.38), (3.39), and (3.40), we obtain $N$ solutions of eigenpairs $\{\beta_{j}^{1},\hat{\phi}_{j}^{1}\}_{j=1}^{N}$ when $N$ collocation points $\{x_{2}^{j}\}_{j=1}^{N}$are used to discretize $x_{2}\in[0,L+d]$. According to [32], ${\rm Im}(\beta_{j}^{1})\geq 0$ so that ${\rm Im}(\sqrt{\beta_{j}^{1}})\geq 0$ and ${\rm Re}(\sqrt{\beta_{j}^{1}})\geq 0$.

Now, inserting each eigenpair into Eq. (3.37) yields two independent solutions $e^{-i\sqrt{\beta_{j}^{1}}x_{1}}$ and $e^{i\sqrt{\beta_{j}^{1}}x_{1}}$. We claim that $u^{sc}$ propagates only towards negative $x_{1}$-axis so that we choose $\psi_{j}=e^{-i\sqrt{\beta_{j}^{1}}x_{1}}$ that propagates towards negative $x_{1}$-axis. To show this, we distinguish two possible cases occurring here. If $\alpha\geq 0$, the given incident wave $u^{in}$ propagates towards positive $x_{1}$-axis so that one easily sees that $\Omega_{1}\subset\Omega_{L}$, which indicates that $u^{sc}=v$ is outgoing in $\Omega_{1}$. If otherwise $\alpha<0$, $u^{in}$ now propagates towards negative $x_{1}$-axis so that $\Omega_{L}\subset\Omega_{1}$, then $u^{sc}$ is an outgoing wave in $\Omega_{1}$ plus a multiple of reflected wave $e^{i\alpha x_{1}+i\beta x_{2}}$ in $\Omega_{1}/\Omega_{L}$ which still propagates towards negative $x_{1}$-axis. Consequently, $\psi_{j}=e^{-i\sqrt{\beta_{j}^{1}}x_{1}}$ such that we get $N$ eigenmodes to approximate $\hat{u}^{s}$,

[TABLE]

where $x_{2}$ is collocated at $\{x_{2}^{j}\}_{j=1}^{N}$.

Repeating the same procedure of variable separation in $\Omega_{2}$, one obtains $N+M$ eigenmodes to approximate $\hat{u}^{s}$,

[TABLE]

where $x_{2}$ is collocated at the common points $\{x_{2}^{j}\}_{j=1}^{N}$ in $[0,L+d]$ and also at $M$ extra points $\{x_{2}^{j}\}_{j=N+1}^{N+M}$ in $V_{h}=[-h,0]$. On $x_{1}=0$ separating $\Omega_{1}$ and $\Omega_{2}$, we have for $j=1,\ldots,N$ that

[TABLE]

by (3.31), and for $j=N+1,\ldots,N+M$ that

[TABLE]

by (3.33).

Eqs. (3.43-3.45) together with the expansions (3.41) and (3.42) give rise to a linear system of $2N+M$ equations for the unknowns $\{c_{j}\}_{j=1}^{N}$ and $\{d_{j}\}_{j=1}^{N+M}$. Solving this linear system, we get $c_{j}$ and $d_{j}$ so that $u^{sc}$ in $\Omega_{1}$ and $\Omega_{2}$ are obtained by (3.41) and (3.42), which eventually gives the total field $u^{tot}$ by (3.31).

As for incident cylindrical waves, the total wave field $\Phi$, post-subtracting the free-space Green’s function $G$, satisfies the Sommerfeld radiation condition (1.7) in the whole domain above $\Gamma$. Therefore, the mode matching procedure described above can be adapted with ease to solve the scattering problem for cylindrical incident waves; we omit the details here.

4 Numerical examples

In this section, we will carry out several numerical experiments to validate our newly proposed radiation condition. In all examples, we assume that the free-space wavelength $\lambda=1$ so that the free-space wavenumber $k_{0}=2\pi$, and the refractive index of the background medium above $\Gamma$ is $n=1$. In setting up the PML, we choose

[TABLE]

for a positive constant $\sigma$.

Example 1. In the first example, we directly analyze the scattering problem for $\Gamma$ in (1.1) with $h=\lambda$. We consider two different incident waves: (1) a plane incident wave with the incident angle $\theta=\frac{\pi}{6}$; (2) a cylindrical incident wave excited by the source point $(0.2,0.2)$. We observe the total wave field $u^{tot}$ in the domain $[-2.5,2.5]\times[-2.5,2.5]$ above $\Gamma$ so that we take $L=2.5$ in the PML.

To validate our numerical solutions, we use $\sigma=70$ and $d=1$ in the PML, and compute $N=280$ eigenmodes in $\Omega_{1}$, and $N+M=420$ eigenmodes in $\Omega_{2}$, i.e., $M=140$ points are used to discretize $V_{h}$, in the NMM method, to get a reference solution $u_{\rm ref}$ for each of the two incident waves, as shown in Figure 5 below.

To illustrate the absorption efficiency of our PML and to validate the newly proposed radiation condition (1.7), we compute the following relative error,

[TABLE]

for different values of $\sigma$ and $d$ in (4.46). The set $S=\{(x,y)|x=0,y=-1,0,2.5\}$ defines where numerical solutions and the reference solution are compared; this choice is typical since it contains all corners of $\Gamma$ and the interior boundary point of the PML.

Figure 6(a) shows the convergence curve of $E_{\rm rel}$ for $\sigma=70$ and for different values of $d$, ranging from $0.001$ to $1$; for a fixed $d=1$, Figure 6(b) shows the convergence curve of $E_{\rm rel}$ for different values of $\sigma$, ranging from $0.1$ to $70$. We observe from Figure 6 that $E_{\rm ref}$ decays exponentially with the PML parameters $\sigma$ and $d$ at the beginning when numerical discretization error is not dominant.

Example 2. In this example, we slightly modify the medium in Example 1 by attaching a penetrable medium of refractive index $n=2$ to the vertical segment $V_{h}$ of $\Gamma$, as shown in Figure 7(c). Again, we compute total wavefields in $[-2.5,2.5]\times[-2.5,2.5]$ for the same two incident waves used in Example 1. To apply the NMM method, we now need to split the medium above the PEC surface into three $x$-uniform regions $\Omega_{1}$, $\Omega_{2}$ and $\Omega_{3}$ separated by $x=0$ and $x=1$. We use $280$ eigenmodes to express the wavefield in $\Omega_{1}$, and $420$ eigenmodes in the other two regions $\Omega_{2}$ and $\Omega_{3}$. Again, reference solutions are obtained by setting $\sigma=70$ and $d=1$ in the PML, as shown in Figure 7 (a) and (b).

To validate the absorption efficiency of our PML, we compute $E_{\rm rel}$ defined in (4.47) for different choices of $\sigma$ and $d$, where now $S=\{(x,y)|x=0,1,y=-1,0,2.5\}$ is chosen to include all corner points.

Figure 8(a) shows the convergence curve of $E_{\rm rel}$ for $\sigma=70$ and for different values of $d$, ranging from $0.001$ to $1$; for a fixed $d=1$, Figure 8(b) shows the convergence curve of $E_{\rm rel}$ for different values of $\sigma$, ranging from $0.1$ to $70$. As in Example 1, we observe from Figure 8 that $E_{\rm ref}$ again decays exponentially with the PML parameters $\sigma$ and $d$ at the beginning when numerical discretization error is not dominant.

Example 3. In the last example, unlike the previous two examples, we make several indentations and put one penetrable medium on the PEC substrate, as shown in Figure 9(c). Two incident waves are considered here: (1) a plane incident wave with incident angle $\theta=\frac{\pi}{6}$; (2) a cylindrical incident wave excited at the source point $(7.2,1.2)$. With totally 3640 eigenmodes to express wavefield in all $x$-uniform regions, we compute two numerical solutions for the two incident waves respectively in $[-2.5,10.5]\times[-2.5,2.5]$, where we take $\sigma=70$ and $d=1$ in the PML. Numerical results are shown in Figure 9(a) and (b).

5 Conclusion

In this paper, we analyzed a sound-soft rough surface scattering problem in two dimensions, where the globally perturbed surface is assumed to consist of two horizontal half lines pointing oppositely and one single vertical line segment connecting their endpoints. For an incident plane wave, we enforced that the scattered wave, post-subtracting reflected plane waves by the two half lines of the scattering surface in their residential regions respectively, satisfies an integral form of SRC condition (1.7) at infinity. With this new radiation condition, we proved uniqueness and existence of weak solutions by a coupling scheme between finite element and integral equation methods. Consequently, this indicates that our new radiation condition is sharper than the ASR condition and UPRC, and generalizes the radiation condition for scattering problems in a locally perturbed half-plane.

Numerically, a NMM method was proposed based on our radiation condition. A perfectly matched layer was setup to absorb the Sommerfeld-type outgoing waves. Since the medium composes of two horizontally uniform regions, we expanded, in either uniform region, the scattered wave in terms of eigenmodes and matched the mode expansions on the common interface between the two uniform regions. This leads to an algebraic linear system and thus yields numerical solutions to our problem. Numerical experiments were carried out to validate the new radiation condition and to show the performance of our numerical method. As the NMM method is only limited to vertical and horizontal interfaces, a PML-based BIE method will be developed in an ongoing work when the scattering surface contains more general curved interfaces. Besides, we shall extend the current work to two-layer media with transmission interface conditions. In that case, the substrates below the interface become penetrable and a new radiation condition modeling the downward propagating waves should be established.

6 Appendix: existence of the improper integral

We need to show that the following improper integral

[TABLE]

exists for a fixed incident angle $\theta\in(\varepsilon,\pi-\varepsilon)$ with $\varepsilon>0$, where $\Phi(x;z)$ is the background Green’s function and

[TABLE]

Note that $\tau$ and $\nu$ denote respectively the tangential and normal directions at the ray ${\mathcal{L}}\subset\mathbb{R}^{2}_{+}$. The integral (6.48) is understood as the limit of

[TABLE]

as $l\rightarrow\infty$. By [9, 10, 41, 42, 22], the two-dimensional background Green’s function $\Phi$, which satisfies the UASP, can be written as

[TABLE]

where $G$ is the free space’s Green’s function in ${\mathbb{R}}^{2}$ and $z_{h}^{*}\in\mathbb{R}_{-}^{2}$ stands for the image of $z$ with respect to the reflection by the line $\{z=(z_{1},z_{2})\in{\mathbb{R}}^{2}:z_{2}=-h\}$. The function $v(x;z)$ is a solution to the Helmholtz equation in $\Omega_{\Gamma}$. Since $G(x,z)-G(x,z_{h}^{*})$ decays with the order $|x_{1}|^{-3/2}$ on $\Gamma$, it follows from ([10]) that the function $v$ belongs to the weighted Sobolev space H_{\varrho}^{1}({\color[rgb]{0.000,0.000,0.000}\Omega_{a}}) (see (1.5) for the definition) for any $a\geq 0$ and $\varrho\in[0,1)$. In view of the asymptotic behavior of the Hankel function for large arguments, one can readily prove that

[TABLE]

as $|x|\rightarrow\infty$ on $\mathcal{L}$, leading to

[TABLE]

Hence, it remains to consider the improper integral (6.48) for $v$ in place of $\Phi$.

Recalling the Upward Propagating Radiation Condition (UPRC), we may represent $v$ as the integral

[TABLE]

The improper integral in the above expression of $v$ can be understood as the duality between $H_{\varrho}^{1/2}({\mathbb{R}})$ and its dual space $H_{-\varrho}^{-1/2}({\mathbb{R}})$ for any $\varrho\in[0,1)$; we refer to [10] for the equivalence of the UPRC and UASR in weighted Sobolev spaces. Obviously,

[TABLE]

where the first term in the integral can be written as

[TABLE]

For $y_{2}=0$ and $x=s\tau(\theta)\in{\mathcal{L}}$, making use of $x\cdot\nu(\theta)=0$ we obtain

[TABLE]

For notational convenience, we write the distance $|x-y|$ in the previous relation as

[TABLE]

This implies that

[TABLE]

Hence, the right hand side of (6) can be bounded by

[TABLE]

where the constant $M>0$ is uniform in all $y_{1}\in{\mathbb{R}}$. This indicates the finiteness of

[TABLE]

By Fubini’s theorem, we can switch the order of integrations to obtain

[TABLE]

where we have defined the following integrals (cf. (6))

[TABLE]

As for $I_{2}^{l}$ and $I_{3}^{l}$, we can easily get the following estimates

[TABLE]

and

[TABLE]

respectively. Below we shall prove the boundedness of $I_{1}^{l}(y_{1})$ for all $y_{1}\in{\mathbb{R}}$ and $l>s_{0}$. For sufficiently large $s$, we can find a positive constant $M$ such that

[TABLE]

Hence, by (6.56),

[TABLE]

The right hand side of (6) is obviously bounded for $|y_{1}|\leq 1$, because of the first relation in (6.52). To derive an upper bound uniformly for $|y_{1}|>1$, we introduce a new variable $t(s)=d(s,y_{1})+s$. Then one can easily check that

[TABLE]

Thus, using change of variables we find

[TABLE]

Since $t-y_{1}\cos\theta\geq|y_{1}|(1-\cos\theta)=2|y_{1}|\sin^{2}(\theta/2)$ and $d(t,y_{1})\geq\max(t\sin\theta,|y_{1}|\sin\theta)$, we have

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

which is bounded as well for sufficiently large $y_{1}$, due to the fact $|\arctan x|\leq\pi/2$ and $t\geq t(s_{0})\geq s_{0}$. Consequently, by (6), (6), (6.63) and (6), we can conclude that $|I_{1}^{l}(y_{1})|$ is uniformly bounded for all $y_{1}\in{\mathbb{R}}$. This together with boundedness of $I_{2}^{l}$ and $I_{3}^{l}$ implies that $I_{j}^{l}\in L_{-\varrho}({\mathbb{R}})$ with $\varrho>1/2$ and thus lies in the dual space $H^{-1/2}_{-\varrho}({\mathbb{R}})$ of $v(y_{1},0;z)$ for $\varrho\in(1/2,1)$. By the dominated convergence theorem, we can now claim that the integral in (6.54) exists. This proves the boundedness of the improper integral (6.48).

Acknowledgement

WL thanks Prof. Ya Yan Lu for useful discussion. The work of G. Hu is supported by the NSFC grant (No. 11671028) and NSAF grant (No. U1530401).