An adaptive finite element DtN method for the elastic wave scattering by biperiodic structures

Consider the scattering of a time-harmonic elastic plane wave by a bi-periodic rigid surface. The displacement of elastic wave motion is modeled by the three-dimensional Navier equation in an open domain above the surface. Based on the Dirichlet-to-Neumann (DtN) operator, which is given as an infinite series, an exact transparent boundary condition is introduced and the scattering problem is formulated equivalently into a boundary value problem in a bounded domain. An a posteriori error estimate based adaptive finite element DtN method is proposed to solve the discrete variational problem where the DtN operator is truncated into a finite number of terms. The a posteriori error estimate takes account of the finite element approximation error and the truncation error of the DtN operator which is shown to decay exponentially with respect to the truncation parameter. Numerical experiments are presented to illustrate the effectiveness of the proposed method.


Introduction
This paper concerns the scattering of a time-harmonic elastic plane wave by a bi-periodic surface in three dimensions. Due to the wide and significant applications in seismology and geophysics, the elastic wave scattering problems have received ever increasing attention in both mathematical and engineering communities [1,2,28]. Compared with the acoustic and electromagnetic wave scattering problems, the elastic wave scattering problems are less studied due to the fact that the elastic wave consists of coupled compressional and shear wave components with different wavenumbers, which makes the analysis of the problems more complicated. In addition, there are two challenges for the elastic surface scattering problem: the solution may have singularity due to a possible non-smooth surface; the problem is imposed in an open domain. In this paper, we intend to address both of these two issues by proposing an a posteriori error estimate based adaptive finite element method with the transparent boundary condition.
The a posteriori error estimates are computable quantities from numerical solutions. They can be used to measure the solution errors of discrete problems without requiring any a priori information of exact solutions [4,31]. Since the a posteriori error estimate based adaptive finite element method has the ability to control the error and to asymptotically optimize the approximation, it is crucial for mesh modification such as refinement and coarsening [17,34]. The method has become an important numerical tool for solving boundary value problems of partial differential equations, especially for those where the solutions have singularity or multiscale phenomena.
The key of overcoming the second issue is to reformulate the open domain problem into a boundary value problem in a bounded domain without generating artificial wave reflection. One possible approach is to make use of the perfectly matched layer (PML) techniques. The basic idea of the PML is to surround the domain of interest by a layer of finite thickness of fictitious medium that may attenuate the waves propagating 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. Due to the effectiveness and simplicity, since Bérenger proposed the technique to solve the time domain Maxwell equations [7], it has undergone a tremendous development of designing various PML methods to solving a wide range of open domain scattering problems [6, 8-10, 15, 16, 19, 20, 27]. Combined with adaptive finite element methods, the PML method has been investigated to solve the two-and three-dimensional obstacle scattering problems [11,12,14] and the two-and three-dimensional diffraction grating problems [5,13,22]. The a posteriori error estimates based adaptive finite element PML methods take account of the finite element discretization errors and the PML truncation errors which decay exponentially with respect to the PML parameters.
Alternatively, another effective approach to truncate the open domain is to construct the Dirichletto-Neumann (DtN) map and introduce the transparent boundary condition to enclose the domain of interest. Since the DtN operator is exact, the transparent boundary condition can be imposed on the boundary which is chosen as close as possible to the scattering structure. Compared to the PML method, the DtN method can reduce the size of the computational domain. As a viable alternative to the PML method, the adaptive finite element DtN methods have also been developed recently to solve many two-and three-dimensional scattering problems, such as the acoustic scattering problems [23,25,35], the three-dimensional electromagnetic scattering problem [26], and the two-dimensional elastic wave scattering problems [29,30].
This paper concerns the numerical solution of the elastic wave scattering by biperiodic structures in three dimensions. It is a non-trivial extension of the elastic wave scattering by periodic structures in two dimensions [29]. There are two challenges for the three-dimensional problem. First, the Helmholtz decomposition of the elastic wave equation gives two two-dimensional Helmholtz equations in the two-dimensional case; however, for the Helmholtz decomposition in the three-dimensional case, we have to consider a three-dimensional Helmholtz equation and a three-dimensional Maxwell equation, which makes the analysis much more complicated. Second, from the computational point of view, it is much more time-consuming to solve the three-dimensional problem than to solve the two-dimensional problem.
Specifically, we consider the scattering of a time-harmonic plane elastic wave by a biperiodic rigid surface. The elastic wave propagation is modeled by the three-dimensional Navier equation in the open domain above the scattering surface. By the Helmholtz decomposition, a DtN operator is constructed in terms of Fourier series expansions for the compressional and shear wave components, then an exact transparent boundary condition is introduced to reduce the open domain problem into an equivalent boundary value problem in a bounded domain. The nonlocal DtN operator needs to be truncated into a sum of finitely many terms in actual computation. However, it is known that the convergence of the truncated DtN operator could be arbitrary slow to the original DtN operator in the operator norm [21]. By carefully examining the properties of the exact solution, we observe that the truncated DtN operator converges exponentially to the original DtN operator when acting on the solution of the elastic wave equation, which enables the analysis of exponential convergence for this work. Combined with the truncated DtN operator and finite element method, the discrete problem is studied. We develop a new duality argument to deduce the a posteriori error estimate. The a posteriori error estimate takes account of the finite element approximation error and the DtN operator truncation error which is shown to decay exponentially with respect to the truncation parameter N . Moreover, an a posteriori error estimate based adaptive finite element algorithm is presented to solve the discrete problem, where the estimate is used to design the 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 turns out not to be sensitive to the given tolerance of accuracy. Numerical examples are presented to demonstrate the effectiveness of the proposed method.
The paper is organized as follows. In Section 2, the model equation is introduced for the scattering problem. Section 3 concerns the variational problem. By the Helmholtz decomposition, the DtN operator is constructed and the transparent boundary condition is introduced to reformulate the scattering problem into a boundary value problem in a bounded domain, and the corresponding weak formulation is presented. In Section 4, the discrete problem is studied by using the finite element method with the truncated DtN operator. Section 5 is devoted to the a posteriori error analysis for the discrete problem and the exponential convergence is proved for the truncated DtN operator. In Section 6, an adaptive finite element algorithm is described and numerical experiments are carried out to illustrate the competitive behavior of the proposed method. The paper is concluded with some general remarks and directions for future work in Section 7.

Problem formulation
Consider the scattering of a time-harmonic plane elastic wave by a rigid biperiodic surface. Due to the biperiodic structure, the scattering problem can be restricted into a single biperiodic cell, as shown in Figure 1. Let S " tx " px 1 , x 2 , x 3 q J P R 3 : px 1 , x 2 q P p0, Λ 1 qˆp0, Λ 2 q, x 3 " f px 1 , x 2 qu be the scattering surface, where f is a Lipschitz continuous biperiodic function with periods Λ 1 and Λ 2 in the x 1 and x 2 directions, respectively. Denote the open space above S by which is assumed to be filled with an isotropic and homogeneous elastic medium. The medium may be characterized by the Lamé parameters λ, µ and the mass density ρ which is assumed to be unit for simplicity. Furthermore, we assume that the Lamé constants satisfy µ ą 0, λ`µ ą 0. Define where h is a constant satisfying h ą max px 1 ,x 2 qPp0,Λ 1 qˆp0,Λ 2 q f px 1 , x 2 q. Denote by Ω the bounded domain enclosed by S and Γ h , i.e., Ω " tx P R 3 : px 1 , x 2 q P p0, Λ 1 qˆp0, Λ 2 q, f px 1 , x 2 q ă x 3 ă hu.
Let a compressional plane wave u inc pxq " qe iκ 1 q¨x be sent from the above to impinge the surface, where q " psin θ 1 cos θ 2 , sin θ 1 sin θ 2 ,´cos θ 1 q J , θ 1 P r0, π{2q and θ 2 P r0, 2πs are the incident angles, and κ 1 " ω{ ? λ`2µ is the compressional wavenumber with ω being the angular frequency. We mention that the results are the same for the incidence of a shear plane wave u inc pxq " pe iκ 2 q¨x , where p is a unit vector satisfying p¨q " 0 and κ 2 " ω{ ? µ is the shear wavenumber, or a linear combination of the shear and compressional plane waves. Denote the displacement of the scattered wave by u, which satisfies the elastic wave equation Since the surface is assumed to be elastically rigid, we have u "´u inc on S.
Define a trace function space H s pΓ h q, s P R`by where the norm is given by It is clear that the dual space of H s pΓ h q is H´spΓ h q with respect to the scalar product in L 2 pΓ h q given by Let H 1 qp pΩq, H 1 S,qp pΩq and H s pΓ h q be the Cartesian product spaces equipped with the corresponding 2-norms of H 1 qp pΩq, H 1 S,qp pΩq and H s pΓ h q, respectively. Throughout the paper, the notation a À b stands for a ď Cb, where C is a positive constant whose value is not required but should be clear from the context.

The boundary value problem
In this section, we introduce the DtN operator to reduce the problem (2.1)-(2.2) into a boundary value problem in the bounded domain Ω and present the well-posedness of its variational formulation.
Consider the Helmholtz decomposition u " ∇φ`∇ˆψ, ∇¨ψ " 0 in Ω, (3.1) where φ is a scalar potential function and ψ " pψ 1 , ψ 2 , ψ 3 q is a vector potential function. Substituting (3.1) into the elastic wave equation (2.1), we may verify that φ and ψ satisfy the following Helmholtz equation and the Maxwell equation, respectively: It is also easy to verify from the Helmholtz decomposition (3.1) and the boundary condition (2.2) that φ and ψ satisfy the following coupled boundary conditions on S: where ν is the unit normal vector on S. The potential functions φ and ψ are required to be quasi-periodic in x 1 and x 2 directions with periods Λ 1 and Λ 2 . Hence they have the Fourier series expansions Plugging (3.4) into (3.2) and using the bounded outgoing wave condition, we have from a straightforward calculation that φ and ψ admit the following expansions: It follows from (3.4)-(3.5) that φ n px 3 q " φ n phqe iβ 1n px 3´h q , ψ jn px 3 q " ψ jn phqe iβ 2n px 3´h q , j " 1, 2, 3. (3.7) We observe from (3.5)-(3.6) that β jn is a pure imaginary number and thus φ n and ψ n are known as surface wave modes when |α n | ą κ j . Substituting (3.5) into (3.1), we obtain the representation of the scattered field u in terms of the Fourier coefficients of the potential functions φ and ψ: fi fl φ n phqe iβ 1n px 3´h q » -α 2n ψ 3n phq´β 2n ψ 2n phq β 2n ψ 1n phq´α 1n ψ 3n phq α 1n ψ 2n phq´α 2n ψ 1n phq fi fl e iβ 2n px 3´h q , .
e iαn¨r . (3.8) Noting ∇¨ψ " 0, we may represent conversely the coefficients of the potential functions of φ and ψ by the coefficients of the scattered field u " pu 1 , u 2 , u 3 q J : φ n phq "´i χ n pα 1n u 1n phq`α 2n u 2n phq`β 2n u 3n phqq , (3.9) ψ 1n phq "´i χ nˆ1 κ 2 2 α 1n α 2n pβ 1n´β2n q u 1n phq´α 2n u 3n phq where χ n " |α n | 2`β 1n β 2n . It is easy to verify that χ n ‰ 0 for n P Z 2 . Define a differential operator where e 3 " p0, 0, 1q J . Substituting (3.8)-(3.12) into the differential operator D, we may deduce the DtN operator T u " where the matrix M n is defined as Here β pnq 12 " β 1n´β2n . The details can be found in [24] for the derivation. Based on the DtN operator (3.14), the scattering problem (2.1)-(2.2) can be equivalently reduced to the following boundary value problem: The variational problem of (3.16) is to find u P H 1 qp pΩq with u "´u inc on S such that where the sesquilinear form a : Here A : B " trpAB J q is the Frobenius inner product of two square matrices A and B.
The well-posedness of the variational problem (3.17) was discussed in [18]. It was shown that the variational problem has a unique weak solution for all but a discrete set of frequencies.
Here we simply assume that the variational problem (3.17) admits a unique solution which satisfies the estimates }u} H 1 pΩq À }u inc } H 1{2 pSq À }u inc } H 1 pΩq .
By the general theory of Babuška and Aziz [3], there exists a γ ą 0 such that the following inf-sup condition holds:

The finite element approximation
Let M h be a regular tetrahedral mesh of the domain Ω, where h denotes the maximum diameter of all the elements in M h . To handle the quasi-periodic solution, we assume that the mesh is periodic in both x 1 and x 2 directions, i.e., the surface meshes on the planes x 1 " 0 and x 2 " 0 coincide with the surface meshes on the planes x 1 " Λ 1 and x 2 " Λ 2 , respectively. We also assume for simplicity that S is polygonal to keep from using the isoparametric finite element space and deriving the approximation error of the boundary S in order to avoid being distracted from the main focus of the a posteriori error analysis.
Let V h Ă H 1 qp pΩq be a conforming finite element space, i.e.
where C qp pΩq is the set of all continuous functions satisfying the quasi-periodic boundary condition, m is a positive integer, and P m denotes the set of all polynomials with degree no more than m.
Since the non-local DtN operator (3.14) is given as an infinite series, in practice, it needs to be truncated into a sum of finitely many terms where N ą 0 is a sufficiently large integer. Using (4.1), we arrive at the truncated finite element approximation: find u h N P V h such that u h N "´u inc on S and satisfies the variational problem where V h,S " tv P V h : v " 0 on Su and the sesquilinear form a N : By [33], 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.2) has a unique solution u h N P V h . The details are omitted for brevity since our focus is on the a posteriori error estimate.

The a posteriori error analysis
For any tetrahedral element K P M h , denote by h K its diameter. Define the operator residual in K as Let F h be the set of all the faces on M h . Given any interior face F P F h , which is the common face of tetrahedral element K 1 and K 2 , we define the jump residual across F as where ν j , j " 1, 2 is the unit outward normal vector on the face of K j . For any boundary face F P F h X Γ h , define the jump residual as Denote the four lateral boundary surfaces by For any boundary face F P Γ 10 and the corresponding face F 1 P Γ 11 , if F P K 1 and F 1 P K 2 , then the jump residual is defined as Similarly, for any face F P Γ 20 and its corresponding face F 1 P Γ 21 , the jump residual is defined as For any tetrahedral element K P M h , denote by η K the local error estimator as follows: For convenience, we introduce a weighted H 1 pΩq norm Since µ and λ`µ are positive, it is easy to check that which implies that the weighted H 1 pΩq norm (5.1) is equivalent to the standard H 1 pΩq norm.
The following theorem is the main result of the paper. It presents the a posteriori error estimate between the solutions of the original scattering problem (3.17) and the truncated finite element approximation (4.2).
Theorem 5.1. Let u and u h N be the solutions of the variational problems (3.17) and (4.2), respectively. Then for anyĥ such that max rPR 2 f prq ăĥ ă h and for sufficiently large N , the following a posteriori error estimate holds: where |n| min " minp|n 1 |, |n 2 |q and |n| max " maxp|n 1 |, |n 2 |q.
The a posteriori error (5.2) contains two parts: the finite element discretization error and the truncation error of the DtN operator. Sinceĥ ă h, the latter is almost exponentially decaying. Hence the DtN truncated error can be controlled to be small enough so that it does not contaminate the finite element discretization error.
To prove Theorem 5.1, let us begin with the following trace result in H 1 qp pΩq. The proof can be found in [24,Lemma 3.3].
Lemma 5.2. Let a " min rPR 2 f prq. Then for any u P H 1 qp pΩq the following estimate holds: Denote by ξ " u´u h N the error between the solutions of (3.17) and (4.2), then a simple calculation yields~ξ~2 Due to the equivalence of the weighted norm~¨~H1 pΩq to the standard norm }¨} H 1 pΩq , it suffices to estimate the four terms on the right hand side of (5.3) one by one. The estimates of the first two terms are given in Lemmas 5.3 and 5.4.
Since the proof of Lemma 5.4 is essentially the same as that for [29, Lemma 5.4], we omit it for brevity. The following two lemmas are to estimate the last term in (5.3).
Lemma 5.6. Let Ω 1 " x P R 3 : px 1 , x 2 q P p0, Λ 1 qˆp0, Λ 2 q,ĥ ă x 3 ă h ( . Then for any δ ą 0, there exists a positive constant Cpδq independent of N such that Proof. It follows from the definition of the DtN operator (3.14) that we have On the other hand, there exists a constant C depending only on ω, µ, λ such thaťˇpM n ξ n q¨ξ nˇď C|ξ n | 2 @ |n| max ď minpN˚, N q. For any δ ą 0, it follows from Young's inequality that which gives that Let φpxq " ř nPZ 2 φ n px 3 qe iαn¨r . It is easy to get Hence, we have for any φ P H 1 pΩ 1 q that which completes the proof.
To estimate the third term of (5.3), we introduce the dual problem apv, pq " ż Ω v¨ξ dx @ v P H 1 S,qp pΩq.
Proof. Denote the coefficient matrix of system (5.19) by A n , which is the matrix A n defined in (5.5).
By the method of the variation of parameters, the solution of (5.19) is pZ 1n px 3 q, Z 2n px 3 q, Z 3n px 3 q, ζ n px 3 qq J " Φ n px 3 qC n px 3 q.

(5.29)
Then p has the Helmholtz decomposition p " ∇g`∇ˆq and satisfies the boundary value problem (5.26).
In our experiments, the parameters are chosen as λ " 1, µ " 1, θ 1 " θ 2 " π{6, ω " 2π. The computational domain Ω " p0, 1qˆp0, 1qˆp0, 0.3q. The mesh and surface plots of the amplitude of the scattered field v h are shown in Figure 2. The mesh has 228400 tetrahedrons and the total number of degrees of freedom (DoFs) on the mesh is 253200. The grating efficiencies are displayed in Figure 3, which verifies the conservation of the energy in [24, Theorem 2.1]. Figure 4 shows the curves of log }∇pu´u k q} 0,Ω versus log N k for both the a priori and the a posteriori error estimates, where N k is the total number of DoFs of the mesh. It indicates that the meshes and the associated numerical complexity are quasi-optimal, i.e., log }∇pu´u k q} 0,Ω " OpN´1 {3 k q is valid asymptotically. Example 2. This example concerns the scattering of a time-harmonic compressional plane wave u inc on a flat grating surface with two square bumps, as seen in Figure 5. The parameters are chosen as λ " 1, µ " 1, θ 1 " θ 2 " π{6, ω " 2π. The computational domain is Ω " p0, 1qˆp0, 1qˆp0, 0.6q. Since there is no exact solution for this example, we plot in Figure 6 the curves of log }∇pu´u k q} 0,Ω versus log N k for the a posteriori error estimates, where N k is the total number of DoFs of the mesh. Again, the result shows that the meshes and the associated numerical complexity are quasi-optimal for the proposed method. We also plot the grating efficiencies against the DoFs in Figure 7 to verify the conservation of the energy. Figures 8 and 9 show the meshes and the amplitude of the associated solution for the scattered field u h when the mesh has 346734 tetrahedrons.

conclusion
In this paper, we have presented an adaptive finite element DtN method for the elastic scattering problem in bi-periodic structures. Based on the Helmholtz decomposition, a new duality argument is developed to obtain the a posteriori error estimate. It takes account of 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 elastic wave propagation problems imposed in unbounded domains. Along the line of this work, a possible continuation is to extend our analysis to the adaptive finite element DtN method for solving the three-dimensional obstacle scattering problem and acoustic-elastic interactive problem. The progress will be reported elsewhere on these problems in the future.