elastic wave scattering by biperiodic structures

# Convergence of the PML solution for elastic wave scattering by biperiodic structures

Xue Jiang School of Science, Beijing University of Posts and Telecommunications, Beijing 100876, China. Peijun Li Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA. Junliang Lv School of Mathematics, Jilin University, Changchun 130012, China.  and  Weiying Zheng LSEC, NCMIS, Academy of Mathematics and System Sciences, CAS; University of Chinese Academy of Sciences, Beijing, 100190, China.
###### Abstract.

This paper is concerned with the analysis of elastic wave scattering of a time-harmonic plane wave by a biperiodic rigid surface, where the wave propagation is governed by the three-dimensional Navier equation. An exact transparent boundary condition is developed to reduce the scattering problem equivalently into a boundary value problem in a bounded domain. The perfectly matched layer (PML) technique is adopted to truncate the unbounded physical domain into a bounded computational domain. The well-posedness and exponential convergence of the solution are established for the truncated PML problem by developing a PML equivalent transparent boundary condition. The proofs rely on a careful study of the error between the two transparent boundary operators. The work significantly extend the results from the one-dimensional periodic structures to the two-dimensional biperiodic structures. Numerical experiments are included to demonstrate the competitive behavior of the proposed method.

###### Key words and phrases:
Elastic wave equation, biperiodic structure, perfectly matched layer, transparent boundary condition
###### 2010 Mathematics Subject Classification:
65N30, 78A45, 35Q60
The research of XJ was supported in part by China NSF grant 11401040 and by the Fundamental Research Funds for the Central Universities 24820152015RC17. The research of PL was supported in part by the NSF grant DMS-1151308. The research of JL was partially supported by the China NSF grants 11126040 and 11301214. The author of WZ was supported in part by China NSF 91430215 and by the National Magnetic Confinement Fusion Science Program (2015GB110003).

## 1. Introduction

Scattering theory in periodic structures has many important applications in diffractive optics [8, 7], where the periodic structures are often named as gratings. The scattering problems have been studied extensively in periodic structures by many researchers for all the commonly encountered waves including the acoustic, electromagnetic, and elastic waves [2, 1, 4, 5, 15, 22, 23, 24, 29, 32]. The governing equations of these waves are known as the Helmholtz equation, the Maxwell equations, and the Navier equation, respectively. In this paper, we consider the scattering of a time-harmonic elastic plane wave by a biperiodic rigid surface, which is also called a two-dimensional or crossed grating. The elastic wave scattering problems have received ever-increasing attention in both engineering and mathematical communities for their important applications in geophysics and seismology. The elastic wave motion is governed by the three-dimensional Navier equation. A fundamental challenge of this problem is to truncate the unbounded physical domain into a bounded computational domain. An appropriate boundary condition is needed on the boundary of the truncated domain to avoid artificial wave reflection. We adopt the perfectly matched layer (PML) technique to handle this issue.

The research on the PML technique has undergone a tremendous development since Bérenger proposed a PML for solving the time-dependent Maxwell equations [11]. The basis idea of the PML technique is to surround the domain of interest by a layer of finite thickness fictitious material which absorbs all the waves coming from inside the computational domain. When the waves reach the outer boundary of the PML region, their values are so small that the homogeneous Dirichlet boundary conditions can be imposed. Various constructions of PML absorbing layers have been proposed and investigated for the acoustic and electromagnetic wave scattering problems [10, 12, 19, 20, 21, 26, 28, 31]. In particular, combined with the PML technique, an effective adaptive finite element method was proposed in [6, 16] to solve the two-dimensional diffraction grating problem where the one-dimensional grating structure was considered. Due to the competitive numerical performance, the method was quickly adopted to solve many other scattering problems including the obstacle scattering problems [17, 14] and the three-dimensional diffraction grating problem [9]. However, the PML technique is much less studied for the elastic wave scattering problems [25], especially for the rigorous convergence analysis. We refer to [13, 18] for recent study on convergence analysis of the elastic obstacle scattering problems.

Recently, we have proposed an adaptive finite element method combining with the PML technique to solve the elastic scattering problem in one-dimensional periodic structures [27]. Using the quasi-periodicity of the solution, we develop a transparent boundary condition and formulate the scattering problem equivalently into a boundary value problem in a bounded domain. Following the complex coordinate stretching, we study the truncated PML problem and show that it has a unique weak solution which converges exponentially to the solution of the original scattering problem.

The purpose of this paper is to extend our previous work on the one-dimensional periodic structures in [27] to the two-dimensional biperiodic structures. We point out that the extension is nontrivial because the more complicated three-dimensional Navier equation needs to be considered. The analysis is mathematically more sophisticated and the numerics is computationally more intense. This work presents an important application of the PML method for the scattering problem of the elastic wave equation. The elastic wave equation is complicated due to the coexistence of compressional and shear waves that have different wavenumbers. To take into account this feature, we introduce two potential functions, one scalar and one vector, to split the wave field into its compressional and shear parts via the Helmholtz decomposition. As a consequence, the scalar potential function satisfies the Helmholtz equation while the vector potential function satisfies the Maxwell equation. Using these two potential functions, we develop an exact transparent boundary condition to reduce the scattering problem from an open domain into a boundary value problem in a bounded domain. The energy conservation is proved for the propagating wave modes of the model problem and is used for verification of our numerical results. The well-posedness and exponential convergence of the solution are established for the truncated PML problem by developing a PML equivalent transparent boundary condition. The proofs rely on a careful study of the error between the two transparent boundary operators. Two numerical examples are also included to demonstrate the competitive behavior of the proposed method.

The paper is organized as follows. In section 2, we introduce the model problem of the elastic wave scattering by a biperiodic surface and formulate it into a boundary value problem by using a transparent boundary condition. In section 3, we introduce the PML formulation and prove the well-posedness and convergence of the truncated PML problem. In section 4, we discuss the numerical implementation of our numerical algorithm and present some numerical experiments to illustrate the performance of the proposed method. The paper is concluded with some general remarks in section 5.

## 2. Problem formulation

In this section, we introduce the model problem and present an exact transparent boundary condition to reduce the problem into a boundary value problem in a bounded domain. The energy distribution will be studied for the diffracted propagating waves of the scattering problem.

### 2.1. Navier equation

Let and . Consider the elastic scattering of a time-harmonic plane wave by a biperiodic surface , where is a Lipschitz continuous and biperiodic function with period in . Denote by the open space above . Let be a constant satisfying . Denote and . Let be the open space above .

The propagation of a time-harmonic elastic wave is governed by the Navier equation:

 μΔu+(λ+μ)∇∇⋅u+ω2u=0in Ωf, (2.1)

where is the displacement vector of the total elastic wave field, is the angular frequency, and are the Lamé constants satisfying and . Assuming that the surface is elastically rigid, we have

 u=0on Γf. (2.2)

Define

 κ1=ω(λ+2μ)1/2andκ2=ωμ1/2,

which are known as the compressional wavenumber and the shear wavenumber, respectively.

Let the scattering surface be illuminated from above by a time-harmonic compressional plane wave:

 uinc(x)=qeiκ1x⋅q,

where is the propagation direction vector, and are called the latitudinal and longitudinal incident angles satisfying . It can be verified that the incident wave also satisfies the Navier equation:

 μΔuinc+(λ+μ)∇∇⋅uinc+ω2uinc=0in Ωf. (2.3)
###### Remark 2.1.

The scattering surface may be also illuminated by a time-harmonic shear plane wave:

 uinc=peiκ2x⋅q,

where is the polarization vector satisfying . More generally, the scattering surface can be illuminated by any linear combination of the time-harmonic compressional and shear plane waves. For clarity, we take the time-harmonic compressional plane wave as an example since the results and analysis are the same for other forms of the incident wave.

Motivated by uniqueness, we are interested in a quasi-periodic solution of , i.e., is biperiodic in and with periods and , respectively. Here with . In addition, the following radiation condition is imposed: the total displacement consists of bounded outgoing waves plus the incident wave in .

We introduce some notation and Sobolev spaces. Let be a vector function. Define the Jacobian matrix of :

 ∇u=⎡⎢ ⎢⎣∂x1u1∂x2u1∂x3u1∂x1u2∂x2u2∂x3u2∂x1u3∂x2u3∂x3u3⎤⎥ ⎥⎦.

Define a quasi-biperiodic functional space

 H1qp(Ω)={u∈H1(Ω): u(x1+n1Λ1,x2+n2Λ2,x3) =u(x1,x2,x3)ei(n1α1Λ1+n2α2Λ2),n=(n1,n2)⊤∈Z2},

which is a subspace of with the norm . For any quasi-biperiodic function defined on , it admits the Fourier series expansion:

 u(r,h)=∑n∈Z2u(n)(h)eiα(n)⋅r,

where , and

 u(n)(h)=1Λ1Λ2∫Λ10∫Λ20u(r,h)e−iα(n)⋅rdr.

We define a trace functional space with the norm given by

 ∥u∥Hs(Γh)=(Λ1Λ2∑n∈Z2(1+|α(n)|2)s|u(n)(h)|2)1/2.

Let and be the Cartesian product spaces equipped with the corresponding 2-norms of and , respectively. It is known that is the dual space of with respect to the inner product

 ⟨u,v⟩Γh=∫Γhu⋅¯vdr,

where the bar denotes the complex conjugate.

### 2.2. Boundary value problem

We wish to reduce the problem equivalently into a boundary value problem in by introducing an exact transparent boundary condition on .

The total field consists of the incident field and the diffracted field , i.e.,

 u=uinc+v. (2.4)

Subtracting (2.3) from (2.1) and noting (2.4), we obtain the Navier equation for the diffracted field :

 μΔv+(λ+μ)∇∇⋅v+ω2v=0in Ωh. (2.5)

For any solution of (2.5), we introduce the Helmholtz decomposition to split it into the compressional and shear parts:

 v=∇ϕ+∇×ψ,∇⋅ψ=0, (2.6)

where is a scalar potential function and is a vector potential function. Substituting (2.6) into (2.5) gives

 (λ+2μ)∇(Δϕ+ω2ϕ)+μ∇×(Δψ+ω2ψ)=0,

which is fulfilled if and satisfy the Helmholtz equation:

 Δϕ+κ21ϕ=0,Δψ+κ22ψ=0. (2.7)

It follows from and (2.7) that the vector potential function satisfies the Maxwell equation:

 ∇×(∇×ψ)−κ22ψ=0.

Since is a quasi-biperiodic function, we have from (2.6) that and are also quasi-biperiodic functions. They have the Fourier series expansions:

 ϕ(x)=∑n∈Z2ϕ(n)(x3)eiα(n)⋅r,ψ(n)(x)=∑n∈Z2ψ(n)(x3)eiα(n)⋅r.

Plugging the above Fourier series into (2.7) yields

 d2ϕ(n)(x3)dx23+(β(n)1)2ϕ(n)(x3)=0,d2ψ(n)(x3)dx23+(β(n)2)2ψ(n)(x3)=0,

where

 β(n)j=⎧⎨⎩(κ2j−|α(n)|2)1/2,|α(n)|<κj,i(|α(n)|2−κ2j)1/2,|α(n)|>κj. (2.8)

Note that . We assume that for all to exclude all possible resonances. Noting (2.8) and using the bounded outgoing radiation condition, we obtain

 ϕ(n)(x3)=ϕ(n)(h)eiβ(n)1(x3−h),ψ(n)(x3)=ψ(n)(h)eiβ(n)2(x3−h).

Hence we deduce Rayleigh’s expansions of and for :

 ϕ(x)=∑n∈Z2ϕ(n)(h)ei(α(n)⋅r+β(n)1(x3−h)),ψ(x)=∑n∈Z2ψ(n)(h)ei(α(n)⋅r+β(n)2(x3−h)). (2.9)

Combining (2.9) and the Helmholtz decomposition (2.6) yields

 v(x)=i∑n∈Z2 ⎡⎢ ⎢ ⎢⎣α(n)1α(n)2β(n)1⎤⎥ ⎥ ⎥⎦ϕ(n)(h)ei(α(n)⋅r+β(n)1(x3−h)) (2.10)

On the other hand, as a quasi-biperiodic function, the diffracted field has the Fourier series expansion:

 v(r,h)=∑n∈Z2v(n)(h)eiα(n)⋅r. (2.11)

It follows from (2.10)–(2.11) and that we obtain a linear system of algebraic equations for and :

 (2.12)

Solving the above linear system directly via Cramer’s rule gives

 ϕ(n)(h) =−iχ(n)(α(n)1v(n)1(h)+α(n)2v(n)2(h)+β(n)2v(n)3(h)) ψ(n)1(h) =−iχ(n)(α(n)1α(n)2(β(n)1−β(n)2)v(n)1(h)/κ22 +[(α(n)1)2β(n)2+(α(n)2)2β(n)1+β(n)1(β(n)2)2]v(n)2(h)/κ22−α(n)2v(n)3(h)) ψ(n)2(h) =−iχ(n)(−[(α(n)1)2β(n)1+(α(n)2)2β(n)2+β(n)1(β(n)2)2]v(n)1(h)/κ22 −α(n)1α(n)2(β(n)1−β(n)2)v(n)2(h)/κ22+α(n)1v(n)3(h)) ψ(n)3(h) =−iκ22(α(n)2v(n)1(h)−α(n)1v(n)2(h)),

where

 χ(n)=|α(n)|2+β(n)1β(n)2. (2.13)

Given a vector field , we define a differential operator on :

 Dv=μ∂x3v+(λ+μ)(∇⋅v)e3, (2.14)

where . Substituting the Helmholtz decomposition (2.6) into (2.14) and using (2.7), we get

 Dv=μ∂x3(∇ϕ+∇×ψ)−(λ+μ)κ21ϕe3.

It follows from (2.10) that

 (Dv)(n)=−μ⎡⎢ ⎢ ⎢⎣α(n)1β(n)10−(β(n)2)2α(n)2β(n)2α(n)2β(n)1(β(n)2)20−α(n)1β(n)2(β(n)2)2−α(n)2β(n)2α(n)1β(n)20⎤⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ϕ(n)(h)ψ(n)1(h)ψ(n)2(h)ψ(n)3(h)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (2.15)

By (2.12) and (2.15), we deduce the transparent boundary conditions for the diffracted field:

 Dv=Tv:=∑n∈Z2M(n)v(n)(h)eiα(n)⋅ron Γ,

where the matrix

 M(n) =iμχ(n) ⎡⎢ ⎢ ⎢⎣(α(n)1)2(β(n)1−β(n)2)+β(n)2χ(n)α(n)1α(n)2(β(n)1−β(n)2)α(n)1β(n)2(β(n)1−β(n)2)α(n)1α(n)2(β(n)1−β(n)2)(α(n)2)2(β(n)1−β(n)2)+β(n)2χ(n)α(n)2β(n)2(β(n)1−β(n)2)−α(n)1β(n)2(β(n)1−β(n)2)−α(n)2β(n)2(β(n)1−β(n)2)κ22β(n)2⎤⎥ ⎥ ⎥⎦.

Equivalently, we have the transparent boundary condition for the total field on :

 Du=Tu+g,

where

 g=Duinc−Tuinc=−2iω2β(0)1κ1χ(0)(α1,α2,−β(0)2)⊤ei(α1x1+α2x2−β(0)1b).

The scattering problem can be reduced to the following boundary value problem:

 ⎧⎪⎨⎪⎩μΔu+(λ+μ)∇∇⋅u+ω2u=0in Ω,Du=Tu+gon Γh,u=0on Γf. (2.16)

The weak formulation of (2.16) reads as follows: Find such that

 a(u,v)=⟨g,v⟩Γh,∀v∈H1qp(Ω)3, (2.17)

where the sesquilinear form is defined by

 a(u,v)=μ∫Ω∇u:∇¯vdx+(λ+μ)∫Ω(∇⋅u)(∇⋅¯v)dx −ω2∫Ωu⋅¯vdx−⟨Tu,v⟩Γh. (2.18)

Here is the Frobenius inner product of square matrices and .

In this paper, we assume that the variational problem (2.17) admits a unique solution. It follows from the general theory in [3] that there exists a constant such that the following inf-sup condition holds:

 sup0≠v∈H1qp(Ω)3|a(u,v)|∥v∥H1(Ω)3≥γ1∥u∥H1(Ω)3,∀u∈H1qp(Ω)3. (2.19)

### 2.3. Energy distribution

We study the energy distribution for the scattering problem. The result will be used to verify the accuracy of our numerical method for examples where the analytical solutions are not available. In general, the energy is distributed away from the scattering surface through propagating wave modes.

Consider the Helmholtz decomposition for the total field:

 u=∇ϕt+∇×ψt,∇⋅ψt=0. (2.20)

Substituting (2.20) into (2.1), we may verify that the scalar potential function and the vector potential function satisfy

 Δϕt+κ21ϕt=0,∇×(∇×ψt)−κ22ψt=0%in Ωf.

We also introduce the Helmholtz decomposition for the incident field

 uinc=∇ϕinc+∇×ψinc,∇⋅ψinc=0,

which gives explicitly that

 ϕinc=−1κ21∇⋅uinc=−iκ1ei(α⋅r−βx3),ψinc=1κ22∇×uinc=0.

Hence we have

 ϕt=ϕinc+ϕ,ψt=ψ.

Using the Rayleigh expansions (2.9), we get

 ϕt(x) =a0ei(α⋅r−βx3)+∑n∈Z2a(n)1ei(α(n)⋅r+β(n)1x3) (2.21) ψt(x) =∑n∈Z2b(n)ei(α(n)⋅r+β(n)2x3), (2.22)

where

 a0=−iκ1,a(n)1=ϕ(n)(h)e−iβ(n)1h,b(n)=ψ(n)(h)e−iβ(n)2h.

The grating efficiency is defined by

 e(n)1=β(n)1|a(n)1|2β|a0|2,e(n)2=β(n)2|b(n)|2β|a0|2, (2.23)

where and are the efficiency of the -th order reflected modes for the compressional wave and the shear wave, respectively. In practice, the grating efficiency (2.23) can be computed from (2.12) once the scattering problem is solved and the diffracted field is available on .

###### Theorem 2.2.

The total energy is conserved, i.e.,

 ∑n∈U1e(n)1+∑n∈U2e(n)2=1,

where .

###### Proof.

It follows from the boundary condition (2.2) and the Helmholtz decomposition (2.20) that

 ∇ϕt+∇×ψt=0on Γf,

which gives

 ν⋅∇ϕt+ν⋅(∇×ψt)=0,ν×∇ϕt+ν×(∇×ψt)=0.

Here is the unit normal vector on .

Consider the following coupled problem:

 {Δϕt+κ21ϕt=0,∇×(∇×ψt)−κ22ψt=0in Ω,ν⋅∇ϕt+ν⋅(∇×ψt)=0,ν×∇ϕt+ν×(∇×ψt)=0on Γf. (2.24)

It is clear to note that also satisfies the problem (2.24) since the wavenumbers are real. Using Green’s theorem and quasi-periodicity of the solution, we get

 0= ∫Ω(¯ϕtΔϕt−ϕtΔ¯ϕt)dx−∫Ω(¯ψt⋅∇×(∇×ψt)−ψt⋅∇×(∇×¯ψt))dx = ∫Γf(¯ϕt∂νϕt−ϕt∂ν¯ϕt)dγ−∫Γf(¯ψt⋅(ν×∇×ψt)−ψt⋅(ν×∇×¯ψt))dγ +∫Γh(¯ϕt∂x3ϕt−ϕt∂x3¯ϕt)dr−∫Γh(¯ψt⋅(e3×∇×ψt)−ψt⋅(e3×∇×¯ψt))dr. (2.25)

It follows from the integration by parts and the boundary conditions in (2.24) that

 ∫Γf∂νϕt¯ϕtdγ=−∫Γfν⋅(∇×ψt)¯ϕtdγ=∫Γfψt⋅(ν×∇¯ϕt)dγ=−∫Γfψt⋅(ν×(∇×¯ψt))dγ,

which gives after taking the imaginary part of (2.25) that

 Im∫Γh(¯ϕt∂x3ϕt−¯ψt⋅(e3×∇×ψt))dr=0. (2.26)

Let . It is clear to note that for and for . It follows from (2.21) and (2.22) that we have

 ϕt(r,h) =a0ei(α⋅r−βh)+∑n∈U1a(n)1e(iα(n)⋅r+iΔ(n)1h)+∑n∉U1a(n)1e(iα(n)⋅r−Δ(n)1h), ψt(r,h) =∑n∈U2b(n)e(iα(n)⋅r+iΔ(n)2h)+∑n∉U2b(n)e(iα(n)⋅r−Δ(n)2h),

and

 ∂x3ϕt(r,h) =−iβa0ei(α⋅r−βh)+∑n∈U1iΔ(n)1a(n)1e(iα(n)⋅r+iΔ(n)1h)−∑n∉U1Δ(n)1a(n)1e(iα(n)⋅r−Δ(n)1h), e3×∇×ψt(r,h) =∑n∈U2⎡⎢ ⎢⎣iα(n)1b(n)3−iΔ(n)2b(n)1iα(n)2b(n)3−iΔ(n)2b(n)20⎤⎥ ⎥⎦e(iα(n)⋅r+iΔ(n)2h) +∑n∉U2⎡⎢ ⎢⎣iα(n)1b(n)3+Δ(n)2b(n)1iα(n)2b(n)3+Δ(n)2b(n)20⎤⎥ ⎥⎦e(iα(n)⋅r−Δ(n)2h),

where . Substituting the above four functions into (2.26), using the orthogonality of the Fourier series and the divergence free condition, we obtain

 β|a0|2=∑n∈U1Δ(n)1|a(n)1|2+∑n∈U2Δ(n)2|b(n)|2,

which completes the proof. ∎

## 3. The PML problem

In this section, we introduce the PML formulation for the scattering problem and establish the well-posedness of the PML problem. An error estimate will be shown for the solutions between the original scattering problem and the PML problem.

### 3.1. PML formulation

Now we turn to the introduction of an absorbing PML layer. The domain is covered by a PML layer of thickness in . Let be the PML function which is continuous and satisfies

 ρ1=1,ρ2=0for τ0otherwise.

We introduce the PML by complex coordinate stretching:

 ^x3=∫x30ρ(τ)dτ. (3.1)

Let . Introduce the new field

 ^u(x)={uinc(x)+(u(^x)−uinc(^x)),x∈Ωh,u(^x),x∈Ω. (3.2)

It is clear to note that in since in . It can be verified from (2.1) and (3.1) that satisfies

 L(^u−uinc)=0in Ωf.

Here the PML differential operator

 Lu=(w1,w2,w3)⊤,

where

 w1= (λ+2μ)∂2x1x1u1+μ(∂2x2x2u1+ρ−1(x3)∂x3(ρ−1(x3)∂x3u1)) +(λ+μ)(∂2x1x2u2