# An adaptive finite element PML method for the elastic wave scattering problem in periodic structures

###### Abstract.

An adaptive finite element method is presented for the elastic scattering of a time-harmonic plane wave by a periodic surface. First, the unbounded physical domain is truncated into a bounded computational domain by introducing the perfectly matched layer (PML) technique. The well-posedness and exponential convergence of the solution are established for the truncated PML problem by developing an equivalent transparent boundary condition. Second, an a posteriori error estimate is deduced for the discrete problem and is used to determine the finite elements for refinements and to determine the PML parameters. Numerical experiments are included to demonstrate the competitive behavior of the proposed adaptive method.

###### Key words and phrases:

Elastic wave equation, adaptive finite element, perfectly matched layer, a posteriori error estimate###### 2010 Mathematics Subject Classification:

65N30, 78A45, 35Q60## 1. Introduction

The scattering theory in periodic diffractive structures, which are known as diffraction gratings, has many significant applications in optical industry [8, 7]. The time-harmonic problems have been studied extensively in diffraction gratings by many researchers for acoustic, electromagnetic, and elastic waves [2, 1, 4, 5, 15, 22, 23, 24, 29, 32]. The underlying equations of these waves are the Helmholtz equation, the Maxwell equations, and the Navier equation, respectively. This paper is concerned with the numerical solution of the elastic wave scattering problem in such a periodic structure. The problem has two fundamental challenges. The first one is to truncate the unbounded physical domain into a bounded computational domain. The second one is the singularity of the solution due to nonsmooth grating surfaces. Hence, the goal of this work is two fold to overcome these two issues. First, we adopt the perfectly matched layer (PML) technique to handle the domain truncation. Second, we use an a posteriori error analysis and design a finite element method with adaptive mesh refinements to deal with the singularity of the solution.

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 energies are so small that the simple 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, 27, 28, 31]. The PML technique is much less studied for the elastic wave scattering problems [25], especically for the rigorous convergence analysis. We refer to [13, 18] for recent study on convergence analysis of the elastic obstacle scattering problem. 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]. Based on the a posteriori error analysis, the adaptive finite element PML method provides an effective numerical strategy which can be used to solve a variety of wave propagation problems which are imposed in unbounded domains.

In this paper, we explore the possibility of applying such an adaptive finite element PML method to solve the diffraction grating problem of elastic waves. Specifically, we consider the incidence of a time-harmonic elastic plane wave on a one-dimensional grating surface, which is assumed to be elastically rigid. The open space, which is above the surface, is assumed to be filled with a homogeneous and isotropic elastic medium. Using the quasi-periodicity of the solution and the transparent boundary condition, we formulate the scattering problem equivalently into a boundary value problem in a bounded domain. The conservation of energy is proved for the model problem and is used to verify our numerical results when the exact solutions are not available. Following the complex coordinate stretching, we study the truncated PML problem which is an approximation to the original scattering problem. We develop the transparent boundary condition for the truncated PML problem and show that it has a unique weak solution which converges exponentially to the solution of the original scattering problem. Moreover, an a posteriori error estimate is deduced for the discrete PML problem. It consists of the finite element error and the PML modeling error. The estimate is used to design the adaptive finite element algorithm to choose elements for refinements and to determine the PML parameters. Numerical experiments show that the proposed method can effectively overcome the aforementioned two challenges.

This paper presents a nontrivial application of the adaptive finite element PML method for the grating problem from the Helmholtz (acoustic) and Maxwell (electromagnetic) equations to the Navier (elastic) equation. The elastic wave equation is complicated due to the coexistence of compressional and shear waves that have different wavenumbers and propagate at different speeds. In view of this physical feature, we introduce two scalar potential functions to split the wave field into its compressional and shear parts via the Helmholtz decomposition. As a consequence, the analysis is much more sophisticated than that for the Helmholtz equation or the Maxwell equations. We believe that this work not only enriches the range of applications for the PML technique but also is a valued contribution to the family of numerical methods for solving elastic wave scattering problems.

The paper is organized as follows. In section 2, we introduce the model problem of the elastic wave scattering by a periodic surface and formulate it into a boundary value problem by using a transparent boundary condition. The conservation of the total energy is proved for the propagating wave modes. In section 3, we introduce the PML formulation and prove the well-posedness and convergence of the truncated PML problem. Section 4 is devoted to the finite element approximation and the a posteriori error estimate. In section 5, we discuss the numerical implementation of our adaptive algorithm and present some numerical experiments to illustrate the performance of the proposed method. The paper is concluded with some general remarks and directions for future research in section 6.

## 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 reflected propagating waves of the scattering problem.

### 2.1. Navier equation

Consider the elastic scattering of a time-harmonic plane wave by a periodic surface which is assumed to be Lipschitz continuous and elastically rigid. In this work, we consider the two-dimensional problem by assuming that the surface is invariant in the direction. The three-dimensional problem will be studied as a separate work. Figure 1 shows the problem geometry in one period. Let . Denote by the artificial boundary above the scattering surface, where is the period and is a constant. Let be the bounded domain which is enclosed from below and above by and , respectively. Finally, denote by the exterior domain to .

The open space, which is above the grating surface, is assumed to be filled with a homogeneous and isotropic elastic medium with a unit mass density. The propagation of a time-harmonic elastic wave is governed by the Navier equation

(2.1) |

where is the angular frequency, and are the Lamé constants satisfying and , and is the displacement vector of the total field which satisfies

(2.2) |

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

or a time-harmonic shear plane wave

where is the incident angle and

(2.3) |

are the compressional and shear wavenumbers, respectively. It can be verified that the incident wave also satisfies the Navier equation:

(2.4) |

###### Remark 2.1.

Our method works for either the compressional plane incident wave, or the shear plane incident wave, or any linear combination of these two plane incident waves. For clarity, we will take the compressional plane incident wave as an example to present the results in our subsequent analysis.

Motivated by uniqueness, we are interested in a quasi-periodic solution of , i.e., is periodic in with period where . 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 and be a vector and scalar function, respectively. Define the Jacobian matrix of as

and two curl operators

Define a quasi-periodic functional space

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

We define a trace functional space with the norm given by

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

where the bar denotes the complex conjuate.

### 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.,

(2.5) |

Noting (2.5) and subtracting (2.4) from (2.1), we obtain the Navier equation for the diffracted field :

(2.6) |

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

(2.7) |

where and are scalar potential functions. Substituting (2.7) into (2.6) gives

which is fulfilled if satisfy the Helmholtz equation

(2.8) |

where is the wavenumber defined in (2.3).

Since is a quasi-periodic function, we have from (2.7) that is also a quasi-periodic function in the direction with period and it has the Fourier series expansion

(2.9) |

Plugging (2.9) into (2.8) yields

(2.10) |

where

(2.11) |

Note that . We assume that for all to exclude possible resonance. Noting (2.11) and using the bounded outgoing radiation condition, we obtain the solution of (2.10):

which gives Rayleigh’s expansion for :

(2.12) |

Combining (2.12) and the Helmholtz decomposition (2.7) yields

(2.13) |

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

(2.14) |

From (2.14) and (2.13), we obtain a linear system of algebraic equations for :

Solving the above equations via Cramer’s rule gives

(2.15a) | ||||

(2.15b) |

where

(2.16) |

Plugging (2.15) into (2.13), we obtain Rayleigh’s expansion for the diffracted field in :

(2.17) |

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

(2.18) |

By (2.18), and (2.2), we deduce the transparent boundary condition

where the matrix

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

where .

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

(2.19) |

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

(2.20) |

where the sesquilinear form is defined by

(2.21) |

Here is the Frobenius inner product of square matrices and .

### 2.3. Energy distribution

We study the energy distribution for the propagating reflected wave modes of the displacement. The result will be used to verify the accuracy of our numerical method when the analytic solution is not available.

Denote by and the unit normal and tangential vectors on , where and . Let and . We point out that and are the collections of all the propagating modes for the compressional and shear waves, respectively. It is clear to note that for and for .

Consider the Helmholtz decomposition for the total field:

(2.23) |

Substituting (2.23) into (2.1), we may verify that also satisfies the Helmholtz equation

Using the boundary condition (2.2), we have

Correspondingly, we introduce the Helmholtz decomposition for the incident field:

which gives explicitly that

Hence we have

Using the Rayleigh expansions (2.12), we get

(2.24) | ||||

(2.25) |

where

The grating efficiency is defined by

(2.26) |

where and are the efficiency of the -th order reflected modes for the compressional wave and the shear wave, respectively. We have the following conservation of energy.

###### Theorem 2.2.

The total energy is conserved, i.e.,

###### Proof.

Consider the following coupled problem:

(2.27) |

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

(2.28) |

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

which yields after taking the imaginary part of (2.3) that

(2.29) |

It follows from (2.24) and (2.25) that we have

and

Substituting the above four functions into (2.29) and using the orthogonality of Fourier series, we get

which completes the proof. ∎

## 3. The PML problem

In this section, we shall 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. As is shown in Figure 2, the domain is covered by a chunck of PML layer of thickness in . Let be the PML function which is continuous and satisfies

We introduce the PML by complex coordinate stretching:

(3.1) |

Let . Introduce the new field

(3.2) |

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

Here the PML differential operator

Define the PML region

Clearly, we have from (3.2) and (2.2) that the outgoing wave in decays exponentially as . Therefore, the homogeneous Dirichlet boundary condition can be imposed on

to truncate the PML problem. Define the computational domain for the PML problem . We arrive at the following truncated PML problem: Find a quasi-periodic solution such that

(3.3) |

where

Define . The weak formulation of the PML problem (3.3) reads as follows: Find such that on and

(3.4) |

Here for any domain , the sesquilinear form is defined by

We will reformulate the variational problem (3.4) in the domain into an equivalent variational formulation in the domain , and discuss the existence and uniqueness of the weak solution to the equivalent weak formulation. To do so, we need to introduce the transparent boundary condition for the truncated PML problem.

### 3.2. Transparent boundary condition of the PML problem

Let . It is clear to note that satisfies the Navier equation in the complex coordinate

(3.5) |

where with .

We introduce the Helmholtz decomposition to the solution of (3.5):

(3.6) |

where and satisfies the Helmholtz equation

(3.7) |

Due to the quasi-periodicity of the solution, we have the Fourier series expansion

(3.8) |

Substituting (3.8) into (3.7) yields

(3.9) |

The general solutions of (3.9) is

Denote by

(3.10) |

The coefficients and can be uniquely determined by solving the following linear equations

(3.11) |

where we have used the Helmholtz decomposition (3.6) and the homogeneous Dirichlet boundary condition

due to the PML absorbing layer. Solving the linear equations (3.11), we obtain