Convergence of an adaptive finite element DtN method for the elastic wave scattering by periodic structures

TL;DR

This paper develops an adaptive finite element method with a truncated Dirichlet-to-Neumann map for elastic wave scattering by periodic structures, providing error estimates and demonstrating effectiveness through numerical experiments.

Contribution

It introduces a new a posteriori error estimate and an adaptive algorithm for elastic wave scattering problems using a truncated DtN map.

Findings

The method effectively reduces scattering problems to bounded domains.

The a posteriori error estimate guides adaptive mesh refinement.

Numerical results confirm the method's accuracy and efficiency.

Abstract

Consider the scattering of a time-harmonic elastic plane wave by a periodic rigid surface. The elastic wave propagation is governed by the two-dimensional Navier equation. Based on a Dirichlet-to-Neumann (DtN) map, a transparent boundary condition (TBC) is introduced to reduce the scattering problem into a boundary value problem in a bounded domain. By using the finite element method, the discrete problem is considered, where the TBC is replaced by the truncated DtN map. A new duality argument is developed to derive the a posteriori error estimate, which contains both the finite element approximation error and the DtN truncation error. An a posteriori error estimate based adaptive finite element algorithm is developed to solve the elastic surface scattering problem. Numerical experiments are presented to demonstrate the effectiveness of the proposed method.

Tables

Table 1. Table 1. The adaptive finite element DtN method.

$\text{if }{\eta }_{\hat{K}}>\tau \underset{K\in {ℳ}_h}{\max }{\eta }_K,\text{ refine the element }\hat{K}\in {ℳ}_h,$

Equations

$$\mu \Delta {\mathbf{u}}^{\mathrm{inc}}+(\lambda +\mu )\nabla \nabla \cdot {\mathbf{u}}^{\mathrm{inc}}+{\omega }^2{\mathbf{u}}^{\mathrm{inc}}=0\text{in}{\Omega }_f^{+},$$

$$\mu \Delta \mathbf{u}+(\lambda +\mu )\nabla \nabla \cdot \mathbf{u}+{\omega }^2\mathbf{u}=0\text{in}{\Omega }_f^{+}.$$

$$\mathbf{u}=-{\mathbf{u}}^{\mathrm{inc}}\text{on}S.$$

$$\mathbf{u}=\nabla {\varphi }_1+\mathbf{curl}{\varphi }_2,$$

$$\Delta {\varphi }_j+{\kappa }_j^2{\varphi }_j=0\text{in}{\Omega }_f^{+}.$$

$${\partial }_{\nu }{\varphi }_1-{\partial }_{\tau }{\varphi }_2={\mathbf{u}}^{\mathrm{inc}}\cdot \nu ,{\partial }_{\nu }{\varphi }_2+{\partial }_{\tau }{\varphi }_1=-{\mathbf{u}}^{\mathrm{inc}}\cdot \tau \text{on}S.$$

$$H_{\mathrm{qp}}^1(\Omega )=\{u\in H^1(\Omega ):u(\Lambda ,y)=u(0,y)e^{\mathrm{i}\alpha \Lambda }\}.$$

$$u(x,b)=\sum _{n\in \mathbb{Z}}{u}^{(n)}(b)e^{\mathrm{i}{\alpha }_{n}x},u^{(n)}(b)=\frac{1}{\Lambda }{\int }_0^{\Lambda }u(x,b)e^{-\mathrm{i}{\alpha }_{n}x}\mathrm{d}x,{\alpha }_n=\alpha +n\left(\frac{2\pi }{\Lambda }\right).$$

$$H^s(\Gamma )=\left\{u\in L^2(\Gamma ):\Vert u{\Vert }_{H^s(\Gamma )}<\infty \right\},$$

$$\Vert u{\Vert }_{H^s(\Gamma )}={\left(\Lambda \sum _{n\in \mathbb{Z}}{\left(1+{\alpha }_n^2\right)}^s|u^{(n)}(b){|}^2\right)}^{1/2}.$$

$${\varphi }_j(x,y)=\sum _{n\in \mathbb{Z}}{\varphi }_j^{(n)}(b)e^{\mathrm{i}\left({\alpha }_{n}x+{\beta }_j^{(n)}(y-b)\right)},y>b,$$

$${\beta }_j^{(n)}=\{\begin{array}{ll}{\left({\kappa }_j^2-{\alpha }_n^2\right)}^{1/2},& |{\alpha }_n|<{\kappa }_j,\\ \mathrm{i}{\left({\alpha }_n^2-{\kappa }_j^2\right)}^{1/2},& |{\alpha }_n|>{\kappa }_j.\end{array}$$

$${\partial }_y{\varphi }_j(x,b)=\sum _{n\in \mathbb{Z}}\mathrm{i}{\beta }_j^{(n)}{\varphi }_j^{(n)}(b)e^{\mathrm{i}{\alpha }_{n}x}.$$

$$\mathbf{u}(x,y)=\sum _{n\in \mathbb{Z}}(u_1^{(n)}(y),u_2^{(n)}(y){)}^{\top }{e}^{\mathrm{i}{\alpha }_{n}x},y>b,$$

$$\mathcal{B}\mathbf{u}=\mu {\partial }_y\mathbf{u}+(\lambda +\mu )(0,1{)}^{\top }\nabla \cdot \mathbf{u}\text{on}\Gamma .$$

$$\mathcal{B}\mathbf{u}=\mathcal{T}\mathbf{u}:=\sum _{n\in \mathbb{Z}}{M}^{(n)}(u_1^{(n)}(b),u_2^{(n)}(b){)}^{\top }{e}^{\mathrm{i}{\alpha }_{n}x}\text{on}\Gamma ,$$

$$M^{(n)}=\frac{\mathrm{i}}{{\chi }_n}\left[\begin{array}{cc}{\omega }^2{\beta }_1^{(n)}& \mu {\alpha }_n{\chi }_n-{\omega }^2{\alpha }_n\\ {\omega }^2{\alpha }_n-\mu {\alpha }_n{\chi }_n& {\omega }^2{\beta }_2^{(n)}\end{array}\right].$$

$$a(\mathbf{u},\mathbf{v})=0,\forall \mathbf{v}\in {\mathbf{H}}_{S,\mathrm{qp}}^1(\Omega ),$$

$$a(\mathbf{u},\mathbf{v})=\mu {\int }_{\Omega }\nabla \mathbf{u}:\nabla \overline{\mathbf{v}}\mathrm{d}\mathbf{x}+(\lambda +\mu ){\int }_{\Omega }\left(\nabla \cdot \mathbf{u}\right)\left(\nabla \cdot \overline{\mathbf{v}}\right)\mathrm{d}\mathbf{x}$$

$$-{\omega }^2{\int }_{\Omega }\mathbf{u}\cdot \overline{\mathbf{v}}\mathrm{d}\mathbf{x}-{\int }_{\Gamma }\mathcal{T}\mathbf{u}\cdot \overline{\mathbf{v}}\mathrm{d}s.$$

$$\Vert \mathbf{u}{\Vert }_{{\mathbf{H}}^1(\Omega )}\lesssim \Vert {\mathbf{u}}^{\mathrm{inc}}{\Vert }_{{\mathbf{H}}^{1/2}(S)}\lesssim \Vert {\mathbf{u}}^{\mathrm{inc}}{\Vert }_{{\mathbf{H}}^1(\Omega )}.$$

$$\underset{0\ne \mathbf{v}\in {\mathbf{H}}_{\mathrm{qp}}^1(\Omega )}{\sup }\frac{|a(\mathbf{u},\mathbf{v})|}{\Vert \mathbf{v}{\Vert }_{{\mathbf{H}}^1(\Omega )}}\ge \gamma \Vert \mathbf{u}{\Vert }_{{\mathbf{H}}^1(\Omega )},\forall \mathbf{u}\in {\mathbf{H}}_{\mathrm{qp}}^1(\Omega ).$$

$${\mathbf{V}}_h:=\left\{\mathbf{v}\in C(\overline{\Omega }{)}^2:\mathbf{v}{|}_K\in P_m(K{)}^2\text{ for any }K\in {\mathcal{M}}_h,\mathbf{v}(0,y)=e^{-\mathrm{i}\alpha \Lambda }\mathbf{v}(\Lambda ,y)\right\},$$

$$a({\mathbf{u}}^h,{\mathbf{v}}^h)=0,\forall {\mathbf{v}}^h\in {\mathbf{V}}_{h,S},$$

$${\mathcal{T}}_N\mathbf{u}=\sum _{|n|\le N}{M}^{(n)}(u_1^{(n)}(b),u_2^{(n)}(b){)}^{\top }{e}^{\mathrm{i}{\alpha }_{n}x},$$

$$a_N({\mathbf{u}}_N^h,{\mathbf{v}}^h)=0,\forall {\mathbf{v}}^h\in {\mathbf{V}}_{h,S},$$

$$a_N(\mathbf{u},\mathbf{v})=\mu {\int }_{\Omega }\nabla \mathbf{u}:\nabla \overline{\mathbf{v}}\mathrm{d}\mathbf{x}+(\lambda +\mu ){\int }_{\Omega }(\nabla \cdot \mathbf{u})(\nabla \cdot \overline{\mathbf{v}})\mathrm{d}\mathbf{x}$$

$$-{\omega }^2{\int }_{\Omega }\mathbf{u}\cdot \overline{\mathbf{v}}\mathrm{d}\mathbf{x}-{\int }_{\Gamma }{\mathcal{T}}_N\mathbf{u}\cdot \overline{\mathbf{v}}\mathrm{d}s.$$

$$J_e=\mu \nabla {\mathbf{u}}_N^h{|}_{K_1}\cdot {\mathbf{\nu }}_1+(\lambda +\mu )(\nabla \cdot {\mathbf{u}}_N^h{|}_{K_1}){\mathbf{\nu }}_1+\mu \nabla {\mathbf{u}}_N^h{|}_{K_2}\cdot {\mathbf{\nu }}_2+(\lambda +\mu )(\nabla \cdot {\mathbf{u}}_N^h{|}_{K_2}){\mathbf{\nu }}_2,$$

$$J_e=2({\mathcal{T}}_N{\mathbf{u}}_N^h-\mathcal{B}{\mathbf{u}}_N^h).$$

$$J_e$$

Full text

1. Introduction

The scattering theory in periodic structures, which are known as gratings in optics, has many significant applications in micro-optics including the design and fabrication of optical elements such as corrective lenses, anti-reflective interfaces, beam splitters, and sensors [44, 7]. Driven by the optical industry applications, the time-harmonic scattering problems have been extensively studied for acoustic and electromagnetic waves in periodic structures. We refer to [8, 18] and the references cited therein for the mathematical results on well-posedness of the solutions for the diffraction grating problems. Computationally, various numerical methods have been developed, such as boundary integral equation method [49, 41], finite element method [5, 6], boundary perturbation method [15]. Recently, the scattering problems for elastic waves have received much attention due to the important applications in seismology and geophysics [1, 2, 42]. This paper concerns the scattering of a time-harmonic elastic plane wave by a periodic surface. Compared with acoustic and electromagnetic wave equations, the elastic wave equation is less studied due to the complexity of the coexistence of compressional and shear waves with different wavenumbers. In addition, there are two challenges for the scattering problem: the solution may have singularity due to a possible nonsmooth surface; the problem is imposed in an open domain. In this paper, we intend to address both issues.

The first issue can be overcome by using the a posteriori error estimate based adaptive finite element method. A posteriori error estimates are computable quantities from numerical solutions and measure the solution errors of discrete problems without requiring any a priori information of real solutions [4, 39]. They are crucial in designing numerical algorithms for mesh modification such as refinement and coarsening [25, 47]. The aim is to equidistribute the computational effort and optimize the computation. The a posteriori error estimate based adaptive finite element method has the ability of error control and asymptotically optimal approximation property [16, 40]. It has become an important numerical tool for solving differential equations, especially for those where the solutions have singularity or multiscale phenomena.

The second issue concerns the domain truncation. The surface scattering problem is imposed in an open domain, which needs to be truncated into a bounded computational domain. An appropriate boundary condition is required on the boundary of the truncated domain so that no artificial wave reflection occurs. Such a boundary condition is called a non-reflecting boundary condition or a transparent boundary condition (TBC) [11, 27, 28, 29, 30]. Despite a huge amount of work done so far in this aspect, it still remains to be one of the important and active research topics in the computational wave propagation. Since Bérenger proposed a perfectly matched layer (PML) technique to solve Maxwell’s equations [12], the research on PML has undergone a tremendous development due to its effectiveness and simplicity [13, 14, 10, 22, 23, 24, 31, 33, 34, 46]. Various constructions of PML have been proposed for solving a wide range of wave propagation problems. The idea of PML technique is to surround the domain of interest by a layer of finite thickness of fictitious medium that may attenuate the waves coming from inside of the computational domain. When the waves reach the outer boundary of the PML region, their amplitudes are so small that the homogeneous Dirichlet boundary condition can be imposed.

Combined with the PML technique, an adaptive finite element method was proposed in [20] to solve the two-dimensional diffraction grating problem. It was shown that the a posteriori error estimate consists of the finite element discretization error and the PML truncation error which decays exponentially with respect to the PML parameters. Due to the competitive numerical performance, the methods was quickly extended to solve the two- and three-dimensional obstacle scattering problems [19, 17] and the three-dimensional diffraction grating problem [9]. Based on the a posteriori error analysis, the adaptive finite element PML method provides an effective numerical strategy to solve a variety of acoustic, electromagnetic, and elastic wave propagation problems which are imposed in unbounded domains [21, 36].

The Dirichlet-to-Neumann (DtN) method is another approach to handle the domain truncation. The idea is to construct an explicit solution, which is usually given as an infinite Fourier series, in the exterior of the domain of interest. By taking the normal derivative of the solution, the Neumann data can be expressed in terms of the Dirichlet data. This relationship gives the DtN map and can be used as a boundary condition, which is known as the TBC. Since the TBC is exact, the artificial boundary can be put as close as possible to the scattering structures, which can reduce the size of the computational domain.

Recently, as a viable alternative to the PML, the adaptive finite element DtN method has been proposed to solve the scattering problems imposed in open domains, such as the obstacle scattering problems [37, 38], the diffraction grating problems [48]. In this approach, the TBC is applied on the artificial boundary which is chosen to enclose the domain of interest. These TBCs are based on nonlocal DtN maps and are given as infinite Fourier series. Practically, the infinite series needs to be truncated into the sum of finite number of terms by choosing an appropriate truncation parameter $N$. It is known that the convergence of the truncated DtN map could be arbitrarily slow to the original DtN map in the operator norm [35]. To overcome this issue, the duality argument has to be developed to obtain the a posteriori error estimate between the solution of the scattering problem and the finite element solution. Comparably, the a posteriori error estimates consists of the finite element discretization error and the DtN truncation error, which decays exponentially with respect to the truncation parameter $N$.

In this paper, we present an adaptive finite element DtN method for the elastic wave scattering problem in periodic structures. The goal is threefold: (1) prove the exponential convergence of the truncated DtN operator; (2) give a complete a posteriori error estimate; (3) develop an effective adaptive finite element algorithm. This paper significantly extends the work on the acoustic scattering problem [48], where the Helmholtz equation was considered. Apparently, the techniques differ greatly from the existing work because of the complicated transparent boundary condition associated with the elastic wave equation. A related work can be found in [43] for an adaptive finite element DtN method for solving the obstacle scattering problem of elastic waves.

Specifically, we consider the scattering of an elastic plane wave by a one-dimensional rigid periodic surface, where the wave motion is governed by the two-dimensional Navier equation. The open space above the surface is assumed to be filled with a homogeneous and isotropic elastic medium. The Helmholtz decomposition is utilized to reduce the elastic wave equation equivalently into a coupled boundary value problem of the Helmholtz equation. By combining the quasi-periodic boundary condition and a DtN operator, an exact TBC is introduced to reduce the original scattering problem into a boundary value problem of the elastic wave equation in a bounded domain. The discrete problem is studied by using the finite element method with the truncated DtN operator. Based on the Helmholtz decomposition, a new duality argument is developed to obtain an a posteriori error estimate between the solution of the original scattering problem and the discrete problem. The a posteriori error estimate contains the finite element approximation error and the DtN operator truncation error, which is shown to decay exponentially with respect to the truncation parameter. The estimate is used to design the adaptive finite element algorithm to choose elements for refinements and to determine the truncation parameter $N$. Due to the exponential convergence of the truncated DtN operator, the choice of the truncation parameter $N$ is not sensitive to the given tolerance. Numerical experiments are presented to demonstrate the effectiveness of the proposed method.

The outline of the paper is as follows. In Section 2, the model equation is introduced for the scattering problem. In Section 3, the boundary value problem is formulated by using the TBC and the corresponding weak formulation is studied. In Section 4, the discrete problem is considered by using the finite element method with the truncated DtN operator. Section 5 is devoted to the a posterior error estimate. In Section 6, we discuss the numerical implementation of the adaptive algorithm and present two examples to illustrate the performance of the proposed method. The paper is concluded with some general remarks and directions for future work in Section 7.

2. Problem formulation

Consider the scattering of a time-harmonic plane wave by an elastically rigid surface, which is assumed to be invariant in the $z$-axis and periodic in the $x$-axis with period $\Lambda$. Due to the periodic structure, the problem can be restricted into a single periodic cell where $x\in(0,\Lambda)$. Let $\boldsymbol{x}=(x,y)\in\mathbb{R}^{2}$. Denote the surface by $S=\{\boldsymbol{x}\in\mathbb{R}^{2}:y=f(x),\,x\in(0,\Lambda)\},$ where $f$ is a Lipschitz continuous function. Let $\nu$ and $\tau$ be the unit normal and tangent vectors on $S$, respectively. Above $S$, the open space is assumed to be filled with a homogeneous and isotropic elastic medium with unit mass density. Denote $\Omega_{f}^{+}=\{\boldsymbol{x}\in\mathbb{R}^{2}:y>f(x),\,x\in(0,\Lambda)\}.$ Let $\Gamma=\{\boldsymbol{x}\in\mathbb{R}^{2}:y=b,\,x\in(0,\Lambda)\}$ and $\Gamma^{\prime}=\{\boldsymbol{x}\in\mathbb{R}^{2}:y=b^{\prime},\,x\in(0,\Lambda)\}$, where $b$ and $b^{\prime}$ are constants satisfying $b>b^{\prime}>\max_{x\in(0,\Lambda)}f(x)$. Denote $\Omega=\{\boldsymbol{x}\in\mathbb{R}^{2}:f(x)<y<b,\,x\in(0,\Lambda)\}$. The problem geometry is shown in Figure 1.

The incident wave $\boldsymbol{u}^{\rm inc}$ satisfies the two-dimensional elastic wave equation

[TABLE]

where $\omega>0$ is the angular frequency and $\mu,\lambda$ are the Lamé parameters satisfying $\mu>0,\lambda+\mu>0$. Specifically, the incident wave can be the compressional plane wave $\boldsymbol{u}^{\rm inc}(\boldsymbol{x})=\boldsymbol{d}e^{{\rm i}\kappa_{1}\boldsymbol{x}\cdot\boldsymbol{d}}$ or the shear plane wave $\boldsymbol{u}^{\rm inc}(\boldsymbol{x})=\boldsymbol{d}^{\perp}e^{{\rm i}\kappa_{2}\boldsymbol{x}\cdot\boldsymbol{d}}$, where $\boldsymbol{d}=(\sin\theta,-\cos\theta)^{\top},\boldsymbol{d}^{\perp}=(\cos\theta,\sin\theta)^{\top},\theta=(-\pi/2,\pi/2)$ is the incident angle, $\kappa_{1}=\omega/(\lambda+2\mu)^{1/2}$ and $\kappa_{2}=\omega/\mu^{1/2}$ are known as the compressional and shear wavenumbers, respectively. For clarity, we shall take the compressional plane wave as the incident field. The results will be similar if the incident field is the shear plane wave.

Due to the interaction between the incident wave and the surface, the scattered wave is generated and satisfies

[TABLE]

Since the surface $S$ is elastically rigid, the displacement of the total field vanishes and the scattered field satisfies

[TABLE]

For any solution $\boldsymbol{u}$ of (2.1), it has the Helmholtz decomposition

[TABLE]

where $\phi_{j},j=1,2$ are scalar potential functions and ${\bf curl}\phi_{2}=(\partial_{y}\phi_{2},-\partial_{x}\phi_{2})^{\top}$. Substituting (2.3) into (2.1), we may verify that $\phi_{j}$ satisfies the Helmholtz equation

[TABLE]

Taking the dot product of (2.2) with $\nu$ and $\tau$, respectively, yields that

[TABLE]

Let $\alpha=\kappa_{\rm p}\sin\theta$. It is clear to note that $\boldsymbol{u}^{\rm inc}$ is a quasi-periodic function with respect to $x$, i.e., $\boldsymbol{u}^{\rm inc}(x,y)e^{-{\rm i}\alpha x}$ is a periodic function with respect to $x$. Motivated by uniqueness of the solution, we require that the solution $\boldsymbol{u}$ of (2.1)–(2.2) is also a quasi-periodic function of $x$ with period $\Lambda$.

We introduce some notations and functional spaces. Let $H^{1}(\Omega)$ be the standard Sobolev space. Denote a quasi-periodic functional space

[TABLE]

Let $H^{1}_{S,{\rm qp}}(\Omega)=\{u\in H^{1}_{\rm qp}(\Omega):u=0\text{ on }S\}$. Clearly, $H_{\rm qp}^{1}(\Omega)$ and $H^{1}_{S,{\rm qp}}(\Omega)$ are subspaces of $H^{1}(\Omega)$ with the standard $H^{1}$-norm. For any function $u\in H_{\rm qp}^{1}(\Omega)$, it admits the Fourier expansion on $\Gamma$:

[TABLE]

The trace functional space $H^{s}(\Gamma),s\in\mathbb{R}$ is defined by

[TABLE]

where the norm is given by

[TABLE]

Let $\boldsymbol{H}_{{\rm qp}}^{1}(\Omega),\boldsymbol{H}_{S,{\rm qp}}^{1}(\Omega),\boldsymbol{H}^{s}(\Gamma)$ be the Cartesian product spaces equipped with the corresponding 2-norms of $H^{1}_{\rm qp}(\Omega),H_{S,{\rm qp}}^{1}(\Omega),H^{s}(\Gamma)$, respectively. Throughout the paper, the notation $a\lesssim b$ stands for $a\leq Cb$, where $C$ is a positive constant whose value is not required but should be clear from the context.

3. The boundary value problem

The scattering problem (2.1)–(2.2) is formulated in the open domain $\Omega^{+}_{f}$, which needs to be truncated into the bounded domain $\Omega$. An appropriate boundary condition is required on $\Gamma$ to avoid artificial wave reflection.

Let $\phi_{j}$ be the solution of the Helmholtz equation (2.4) along with the bounded outgoing wave condition. It is shown in [42] that $\phi_{j}$ is a quasi-periodic function and admits the Fourier series expansion

[TABLE]

where

[TABLE]

We assume that $\kappa_{j}\neq|\alpha_{n}|$ for $n\in\mathbb{Z}$ to exclude possible resonance. Taking the normal derivative of (3.1) on $\Gamma$ yields

[TABLE]

As a quasi-periodic function, the solution $\boldsymbol{u}(x,y)=(u_{1}(x,y),u_{2}(x,y))^{\top}$ admits the Fourier expansion

[TABLE]

where $u_{j}^{(n)}$ is the Fourier coefficient of $u_{j}$. Define a boundary operator

[TABLE]

It is shown in [36] that the solution of (2.1) satisfies the transparent boundary condition

[TABLE]

where $\mathscr{T}$ is called the Dirichlet-to-Neumann (DtN) operator and $M^{(n)}$ is a $2\times 2$ matrix given by

[TABLE]

Here $\chi_{n}=\alpha_{n}^{2}+\beta_{1}^{(n)}\beta_{2}^{(n)}$.

By the transparent boundary condition (3.3), the variational problem of (2.1)–(2.2) is to find $\boldsymbol{u}\in\boldsymbol{H}_{{\rm qp}}^{1}(\Omega)$ with $\boldsymbol{u}=-\boldsymbol{u}^{\rm inc}$ on $S$ such that

[TABLE]

where the sesquilinear form $a:\boldsymbol{H}^{1}_{{\rm qp}}(\Omega)\times\boldsymbol{H}^{1}_{{\rm qp}}(\Omega)\to\mathbb{C}$ is defined as

[TABLE]

Here $A:B={\rm tr}(AB^{\top})$ is the Frobenius inner product of two square matrices $A$ and $B$.

The well-posedness of the variational problem (3.5) was discussed in [26]. It was shown that the variational problem (3.5) has a unique solution for all frequencies if the surface $S$ is Lipschitz continuous. Hence we may assume that the variational problem (3.5) admits a unique solution and the solution satisfies the estimate

[TABLE]

By the general theory of Babuska and Aziz [3], there exists $\gamma>0$ such that the following inf-sup condition holds

[TABLE]

4. The discrete problem

We consider the discrete problem of (3.5) by using the finite element approximation. Let $\mathcal{M}_{h}$ be a regular triangulation of $\Omega$, where $h$ denotes the maximum diameter of all the elements in $\mathcal{M}_{h}$. Since our focus is on the a posteriori error estimate, for simplicity, we assume that $S$ is polygonal and ignore the approximation error of the boundary $S$. Thus any edge $e\in\mathcal{M}_{h}$ is a subset of $\partial\Omega$ if it has two boundary vertices. Moreover, we require that if $(0,y)$ is a node on the left boundary, then $(\Lambda,y)$ is also a node on the right boundary and vice versa, which allows to define a finite element space whose functions are quasi-periodic respect to $x$.

Let $\boldsymbol{V}_{h}\subset\boldsymbol{H}_{\rm qp}^{1}(\Omega)$ be a conforming finite element space, i.e.,

[TABLE]

where $m$ is a positive integer and $P_{m}(K)$ denotes the set of all polynomials of degree no more than $m$. The finite element approximation to the variational problem (3.5) is to find $\boldsymbol{u}^{h}\in\boldsymbol{V}_{h}$ with $\boldsymbol{u}^{h}=-\boldsymbol{u}^{\rm inc}$ on $S$ such that

[TABLE]

where $\boldsymbol{V}_{h,S}=\{\boldsymbol{v}\in\boldsymbol{V}_{h}:\boldsymbol{v}=0\text{ on }S\}$.

In the variational problem (4.1), the boundary operator $\mathscr{T}$ is defined as an infinite series, in practice, it must be truncated to a sum of finitely many terms as follows

[TABLE]

where $N>0$ is a sufficiently large constant. Using the truncated boundary operator, we arrive at the truncated finite element approximation: Find $\boldsymbol{u}_{N}^{h}\in\boldsymbol{V}_{h}$ such that it satisfies $\boldsymbol{u}_{N}^{h}=-\boldsymbol{u}^{\rm inc}$ on $S$ and the variational problem

[TABLE]

where the sesquilinear form $a_{N}:\boldsymbol{V}_{h}\times\boldsymbol{V}_{h}\to\mathbb{C}$ is defined as

[TABLE]

It follows from [45] that the discrete inf-sup condition of the sesquilinear form $a_{N}$ can be established for sufficient large $N$ and small enough $h$. Based on the general theory in [3], it can be shown that the discrete variational problem (4.3) has a unique solution $\boldsymbol{u}_{N}^{h}\in\boldsymbol{V}_{h}$. The details are omitted for brevity.

5. The a posteriori error analysis

For any triangular element $K\in\mathcal{M}_{h}$, denoted by $h_{K}$ its diameter. Let $\mathcal{B}_{h}$ denote the set of all the edges of $K$. For any $e\in\mathcal{B}_{h}$, denoted by $h_{e}$ its length. For any interior edge $e$ which is the common side of $K_{1}$ and $K_{2}\in\mathcal{M}_{h}$, we define the jump residual across $e$ as

[TABLE]

where $\boldsymbol{\nu}_{j}$ is the unit outward normal vector on the boundary of $K_{j},j=1,2$. For any boundary edge $e\subset\Gamma$, we define the jump residual

[TABLE]

For any boundary edge on the left line segment of $\partial\Omega$, i.e., $e\in\{x=0\}\cap\partial K_{1}$ for some $K_{1}\in\mathcal{M}_{h}$, and its corresponding edge on the right line segment of $\partial\Omega$, i.e., $e^{\prime}\in\{x=\Lambda\}\cap\partial K_{2}$ for some $K_{2}\in\mathcal{M}_{h}$, the jump residual is

[TABLE]

For any triangular element $K\in\mathcal{M}_{h}$, denote by $\eta_{K}$ the local error estimator which is given by

[TABLE]

where $\mathscr{R}$ is the residual operator defined by

[TABLE]

For convenience, we introduce a weighted norm of $\boldsymbol{H}^{1}(\Omega)$ as

[TABLE]

It is easy to check that

[TABLE]

which implies that the weighted norm is equivalent to standard $\boldsymbol{H}^{1}(\Omega)$ norm.

Now we state the main result of this paper.

Theorem 5.1**.**

Let $\boldsymbol{u}$ and $\boldsymbol{u}_{N}^{h}$ be the solutions of the variational problem (3.5) and (4.3), respectively. Then for sufficient large $N$, the following a posteriori error estimate holds

[TABLE]

It is easy to note that the a posteriori error consists of two parts: the finite element discretization error and the truncation error of the DtN operator. We point out that the latter is almost exponentially decaying since $b>b^{\prime}$ and $|\beta_{2}^{(n)}|>0$. In practice, the DtN truncated error can be controlled to be small enough such that it does not contaminate the finite element discretization error.

In the rest of the paper, we shall prove the a posteriori error estimate in Theorem 5.1. First, let’s state the trace regularity for functions in $H^{1}_{\rm qp}(\Omega)$. The proof can be found in [20].

Lemma 5.2**.**

For any $u\in H^{1}_{\rm qp}(\Omega)$, the following estimates hold

[TABLE]

Denote by $\boldsymbol{\xi}=\boldsymbol{u}-\boldsymbol{u}_{N}^{h}$ the error between the solutions of (3.5) and (4.3). It can be verified that

[TABLE]

In the following, we shall discuss the four terms in the right hand side of (5.2). Lemma 5.3 gives the error estimate of the truncated DtN operator. Lemma 5.4 presents the a posteriori error estimate for the finite element approximation and the truncated DtN operator.

Lemma 5.3**.**

Let $\boldsymbol{u}\in\boldsymbol{H}^{1}_{\rm qp}(\Omega)$ be the solution of the variational problem (3.5). For any $\boldsymbol{v}\in\boldsymbol{H}^{1}_{\rm qp}(\Omega)$, the following estimate holds:

[TABLE]

where $C>0$ is a constant independent of $N$.

Proof.

Using (2.3) and (3.1) yields

[TABLE]

It follows from the straightforward calculations that we obtain

[TABLE]

where

[TABLE]

It is clear to note from (3.2) that $\beta_{j}^{(n)}$ is purely imaginary for sufficiently large $|n|$. By the mean value theorem, for sufficiently large $|n|$, there exists $\tau\in({\rm i}\beta_{1}^{(n)},{\rm i}\beta_{2}^{(n)})$ such that

[TABLE]

A simple calculation yields

[TABLE]

and

[TABLE]

which give

[TABLE]

Similarly, we may show that

[TABLE]

Combining the above estimates lead to

[TABLE]

By (3.3) and (4.2), we have from Lemma 5.2 that

[TABLE]

Using (3.6), we get

[TABLE]

which completes the proof. ∎

In the following lemmas, the first two terms in (5.2) are estimated.

Lemma 5.4**.**

Let $\boldsymbol{v}$ be any function in $\boldsymbol{H}^{1}_{S,{\rm qp}}(\Omega)$, the following estimate holds

[TABLE]

Proof.

For any function $\boldsymbol{v}\in H^{1}_{S,{\rm qp}}(\Omega)$, we have

[TABLE]

For any function $\boldsymbol{v}\in\boldsymbol{H}_{S,{\rm qp}}^{1}(\Omega)$ and $\boldsymbol{v}^{h}\in\boldsymbol{V}_{h,S}$, it follows from the integration by parts that

[TABLE]

We take $\boldsymbol{v}^{h}=\Pi_{h}\boldsymbol{v}\in\boldsymbol{V}_{h,S}$, where $\Pi_{h}$ is the Scott–Zhang interpolation operator and has the following interpolation estimates

[TABLE]

Here $\tilde{K}$ and $\tilde{K}_{e}$ are the unions of all the triangular elements in $\mathcal{M}_{h}$, which have nonempty intersection with the element $K$ and the side $e$, respectively. By the Hölder equality, we get from (5.5) that

[TABLE]

which completes the proof. ∎

Lemma 5.5**.**

Let $\hat{M}^{(n)}=-\frac{1}{2}(M^{(n)}+(M^{(n)})^{*})$, where $M^{(n)}$ is defined in (3.4). Then $\hat{M}^{(n)}$ is positive definite for sufficiently large $|n|$.

Proof.

It follows from (3.2) that $\beta_{j}^{(n)}$ is purely imaginary for sufficiently large $|n|$. By (3.4), we have

[TABLE]

Since $\chi_{n}=\alpha_{n}^{2}-(\alpha_{n}^{2}-\kappa_{1}^{2})^{1/2}(\alpha_{n}^{2}-\kappa_{2}^{2})^{1/2}>0$, we get

[TABLE]

A simple calculation yields that

[TABLE]

Since $\kappa_{2}>\kappa_{1}$ and $\alpha^{2}_{n}$ has an order of $n^{2}$ for sufficiently large $|n|$, we obtain

[TABLE]

which gives that $\det\hat{M}^{(n)}>0$ and completes the proof. ∎

Lemma 5.6**.**

Let $\Omega^{\prime}=\{\boldsymbol{x}\in\mathbb{R}^{2}:b^{\prime}<y<b,\,0<x<\Lambda\}$. Then for any $\delta>0$, there exists a positive constant $C(\delta)$ independent of $N$ such that

[TABLE]

Proof.

Using (4.2), we get from a simple calculation that

[TABLE]

By Lemma 5.5, $\hat{M}^{(n)}$ is positive definite for sufficiently large $|n|$. Hence, for fixed $\omega,\lambda,\mu$, there exists $N^{*}$ such that $-\left(\hat{M}^{(n)}\boldsymbol{\xi}^{(n)}\right)\cdot\overline{\boldsymbol{\xi}^{(n)}}\leq 0$ for $n>N^{*}$. Correspondingly, we split $\Re\int_{\Gamma}\mathscr{T}_{N}\boldsymbol{\xi}\cdot\overline{\boldsymbol{\xi}}{\rm d}s$ into two parts:

[TABLE]

where $\sum\limits_{N>|n|>\min(N^{*},N)}\left(\hat{M}^{(n)}\boldsymbol{\xi}_{n}\right)\cdot\overline{\boldsymbol{\xi}_{n}}=0$ if $N>N^{*}$. Since the second part in the right hand side of (5.6) is non-positive, we only need to estimate the first part in the right hand side of (5.6), which has finitely many terms. Hence there exists a constant $C$ depending only on $\omega,\mu,\lambda$ such that $|\left(\hat{M}^{(n)}\boldsymbol{\xi}^{(n)}\right)\cdot\overline{\boldsymbol{\xi}^{(n)}}|\leq C|\boldsymbol{\xi}^{(n)}|^{2}$ for all $|n|\leq\min(N^{*},N)$.

For any $\delta>0$, it follows from Yong’s inequality that

[TABLE]

which gives

[TABLE]

Let $\phi(x,y)=\sum\limits_{n\in\mathbb{Z}}\phi_{n}(y)e^{{\rm i}\alpha_{n}x}$. A simple calculation yields that

[TABLE]

Using the above estimates, we have for any $\phi\in H^{1}(\Omega^{\prime})$ that

[TABLE]

Combining the above estimates, we obtain

[TABLE]

which completes the proof. ∎

To estimate $\int_{\Omega}|\boldsymbol{\xi}|^{2}{\rm d}\boldsymbol{x}$ in (5.2) , we introduce the dual problem

[TABLE]

It can be verified that $\boldsymbol{p}$ is the weak solution of the boundary value problem

[TABLE]

where $\mathscr{T}^{*}$ is the adjoint operator to the DtN operator $\mathscr{T}$.

It requires to explicitly solve the boundary value problem (5.8). We consider the Helmholtz decomposition and let

[TABLE]

where $\zeta_{j},j=1,2$ has the Fourier series expansion

[TABLE]

Consider the following coupled first order ordinary different equations

[TABLE]

It follows from straightforward calculations that the solution is

[TABLE]

It is easy to verify the following estimate

[TABLE]

Let $\boldsymbol{p}$ be the solution of the dual problem (5.8). Then it satisfies the following boundary value problem

[TABLE]

Let function $q_{j},j=1,2$ have the Fourier expansion in $\Omega^{\prime}$:

[TABLE]

The Fourier coefficients $q_{j}^{(n)}$ are required to satisfy the two point boundary value problem

[TABLE]

where $c_{1}=(\lambda+2\mu)^{-1}$ and $c_{2}=\mu^{-1}$, $\zeta_{j}^{(n)}$ are the Fourier coefficients of the potential functions $\zeta_{j}$ for the Helmholtz decomposition of $\boldsymbol{\xi}$ in (5.9).

Lemma 5.7**.**

Let $\boldsymbol{p}=\nabla q_{1}+{\bf curl}q_{2}$. Then $\boldsymbol{p}$ satisfies (5.10).

Proof.

If (5.11) holds, then it is easy to check that

[TABLE]

Noting $\boldsymbol{p}=\nabla q_{1}+{\bf curl}q_{2}$, we obtain

[TABLE]

Next is to verify that the boundary condition on $y=b$. Assume that $\boldsymbol{p}$ admits the Fourier expansion $\boldsymbol{p}=\sum\limits_{n\in\mathbb{Z}}(p_{1}^{(n)}(y),p_{2}^{(n)}(y))^{\top}e^{{\rm i}\alpha_{n}x}$. It follows from the Helmholtz decomposition that

[TABLE]

which gives

[TABLE]

A straightforward calculation yields that

[TABLE]

Evaluating the above equations at $y=b$, we get

[TABLE]

Noting $\zeta_{j}^{(n)}(b)=0$, we have from (5.11) that $q^{(n)^{\prime\prime}}_{j}(b)=-(\kappa_{j}^{2}-\alpha_{n}^{2})q^{(n)}_{j}(b)$. Hence

[TABLE]

On the other hand, we have

[TABLE]

which shows $\mathscr{B}\boldsymbol{p}=\mathscr{T}^{*}\boldsymbol{p}$ and completes the proof. ∎

It follows from the classic theory of second order differential equations that the solution of the system

[TABLE]

is

[TABLE]

Lemma 5.8**.**

Let $\boldsymbol{p}=(p_{1},p_{2})^{\top}$ be the solution of the dual problem problem (5.7). For sufficiently large $|n|$, the following estimate hold

[TABLE]

where $p_{j}^{(n)}$ is the Fourier coefficient of $p_{j},j=1,2$.

Proof.

Evaluating (5) at $y=b$ yields

[TABLE]

Taking the derivative of $q^{(n)}_{j}$ with respect to $y$ in (5) and then evaluating at $y=b^{\prime}$, we have

[TABLE]

which is equivalent to

[TABLE]

where

[TABLE]

It follows from Lemma 5.7 and the Helmholtz decomposition $\boldsymbol{p}=\nabla q_{1}+{\bf curl}q_{2}$ that

[TABLE]

which gives

[TABLE]

Substituting the boundary condition

[TABLE]

into the Helmholtz decomposition $\boldsymbol{p}=\nabla q_{1}+{\bf curl}q_{2}$, i.e.,

[TABLE]

we obtain

[TABLE]

By (5),

[TABLE]

where

[TABLE]

Combining the above equations leads to

[TABLE]

where $P^{(n)}$ is defined in (5.3).

Recall that

[TABLE]

Since $s-b\geq 2b^{\prime}-b-s$ and $|\alpha_{n}|\sim|n|,|\beta_{j}^{(n)}|\sim|n|$ for sufficiently large $|n|$, we have from (5.4) and the mean-value theorem that

[TABLE]

Combining the above estimates yields

[TABLE]

Following the similar steps of the estimate for $\eta_{j}^{(n)}$, we can show that

[TABLE]

which gives

[TABLE]

Since for sufficiently large $|n|$, we have

[TABLE]

Hence

[TABLE]

which proves

[TABLE]

The proof is completed. ∎

Taking $\boldsymbol{v}=\boldsymbol{\xi}$ in (5.7), we have

[TABLE]

By Lemma 5.8, we obtain

[TABLE]

Following the similar proof in [37, eq. (30)], we may show that

[TABLE]

It follows from the Cauchy–Schwarz inequality that

[TABLE]

Noting that the function $t^{4}e^{-2t}$ is bounded on $(0,+\infty)$, we have

[TABLE]

Substituting (5.16) into $I_{2}$, we get

[TABLE]

A simple calculation yields

[TABLE]

It is easy to note that

[TABLE]

Then

[TABLE]

Therefore,

[TABLE]

Plugging (5.17) to (5.15), we obtain

[TABLE]

Now, we prove Theorem 5.1.

Proof.

By Lemma 5.3, Lemma 5.4, and Lemma 5.6, we have

[TABLE]

where $C_{1},C_{2},C(\delta)$ are positive constants. From (5.1), by choosing $\delta$ such that $\frac{\delta}{\min(\mu,\omega^{2})}<\frac{1}{2}$, we get

[TABLE]

It follows from (5.14) and (5.18) that

[TABLE]

Taking sufficiently large $N$ such that $\frac{2\left(C_{2}+C(\delta)\right)}{N}\frac{1}{\min(\mu,\omega^{2})}<1$ and substituting (5.20) into (5), we obtain

[TABLE]

The proof is completed by noting the equivalence of the norms $\vvvert\cdot\vvvert_{\boldsymbol{H}^{1}(\Omega)}$ and $\|\cdot\|_{\boldsymbol{H}^{1}(\Omega)}$. ∎

6. Numerical experiments

In this section, we introduce the algorithmic implementation of the adaptive finite element DtN method and present two numerical examples to demonstrate the effectiveness of the proposed method.

6.1. Adaptive algorithm

Our implementation is based on the FreeFem [32]. The first-order linear element is used to solve the problem. It is shown in Theorem 5.1 that the a posteriori error consists of two parts: the finite element discretization error $\epsilon_{h}$ and the DtN operator truncation error $\epsilon_{N}$, where

[TABLE]

In the implementation, we choose the parameters $b,b^{\prime}$ and $N$ based on (6.1) to make sure that the DtN operator truncation error is smaller than the finite element discretization error. In the following numerical experiments, $b^{\prime}$ is chosen such that $b^{\prime}=\max_{x\in(0,\Lambda)}f(x)$ and $N$ is the smallest positive integer that makes $\epsilon_{N}\leq 10^{-8}$. The adaptive finite element algorithm is shown in Table 1.

6.2. Numerical experiments

We report two examples to illustrate the numerical performance of the proposed method. The first example concerns the scattering by a flat surface and has an exact solution; the second example is constructed such that the solution has corner singularity.

Example 1. We consider the simplest periodic structure, a straight line, where the exact solution is available. Let $S=\{y=0\}$ and take the artificial boundary $\Gamma=\{y=0.25\}$. The space above the flat surface is filled with a homogenenous and isotropic elastic medium, which is characterized by the Lamé constants $\lambda=2$, $\mu=1$. The rigid surface is impinged by the compressional plane wave $\boldsymbol{u}^{\rm inc}=\boldsymbol{d}e^{{\rm i}\kappa_{1}\boldsymbol{x}\cdot\boldsymbol{d}}$, where the incident angle is $\theta=\pi/3$. The compressional and shear wavenumbers are $\kappa_{1}=\omega/2$ and $\kappa_{2}=\omega$, respectively, where $\omega$ is the angular frequency. It can be verified that the exact solution is

[TABLE]

where $\alpha=\kappa_{1}\sin\theta,\beta=\kappa_{1}\cos\theta,\gamma=(\kappa_{2}^{2}-\alpha^{2})^{1/2}$. The period $\Lambda=0.5$. Figure 2 shows the curves of $\log e_{h}$ versus $\log{\rm DoF}_{h}$ with different angular frequencies, where $e_{h}=\|\boldsymbol{u}-\boldsymbol{u}_{N}^{h}\|_{\boldsymbol{H}^{1}(\Omega)}$ is the a priori error and ${\rm DoF}_{h}$ stands for the degree of freedom or the number of nodal points. It indicates that the meshes and the associated numerical complexity are quasi-optimal, i.e., $e_{h}=\mathcal{O}({\rm DoF}_{h}^{-1/2})$ holds asymptotically.

Example 2. This example concerns the scattering of the compressional plane wave by a piecewise linear surface, which has multiple sharp angles. The incident wave $\boldsymbol{u}^{\rm inc}$ and the parameters are chosen the same as Example 1, i.e., $b=0.25,\Lambda=0.5,\theta=\pi/3,\lambda=1,\mu=2$. Clearly, the solution has singularity around the corners of the surface. Since there is no exact solution for this example, we plot in Figure 3 the curves of $\log\epsilon_{h}$ versus $\log{\rm DoF}_{h}$ at different angular frequencies, where $\epsilon_{h}$ is the a posteriori error. Again, it indicates that the meshes and the associated numerical complexity are quasi-optimal, i.e., $\epsilon_{h}=\mathcal{O}({\rm DoF}_{h}^{-1/2})$. Figure 4 plots the contour of the magnitude of the numerical solution and its corresponding mesh at the angular frequency $\omega=2$. It is clear to note that the algorithm does capture the solution feature and adaptively refines the mesh around the corners where solution displays singularity.

7. Conclusion

In this paper, we have presented an adaptive finite element DtN method for the elastic scattering problem in periodic structures. Based on the Helmholtz decomposition, a new duality argument is developed to obtain the a posteriori error estimate. It contains both the finite element discretization error and the DtN operator truncation error, which is shown to decay exponentially with respect to the truncation parameter. Numerical results show that the proposed method is effective and accurate. This work provides a viable alternative to the adaptive finite element PML method for solving the elastic surface scattering problem. It also enriches the range of choices available for solving wave propagation problems imposed in unbounded domains. One possible future work is to extend our analysis to the adaptive finite element DtN method for solving the three-dimensional elastic surface scattering problem, where a more complicated TBC needs to be considered.