An adaptive finite element PML method for the acoustic-elastic interaction in three dimensions

TL;DR

This paper develops an adaptive finite element PML method with an a posteriori error estimate for 3D acoustic-elastic interaction problems, ensuring well-posedness, exponential convergence, and demonstrating its effectiveness through numerical experiments.

Contribution

It introduces a novel adaptive finite element PML approach with error estimation for 3D acoustic-elastic scattering, improving accuracy and efficiency.

Findings

The method achieves exponential convergence.

Numerical results validate the method's effectiveness.

The approach is competitive for complex scattering problems.

Abstract

Consider the scattering of a time-harmonic acoustic incident wave by a bounded, penetrable, and isotropic elastic solid, which is immersed in a homogeneous compressible air or fluid. The paper concerns the numerical solution for such an acoustic-elastic interaction problem in three dimensions. An exact transparent boundary condition (TBC) is developed to reduce the 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 using a PML equivalent TBC. An a posteriori error estimate based adaptive finite element method is developed to solve the scattering problem. Numerical experiments are included to demonstrate the competitive behavior of the proposed method.

Tables

Table 1. Table 1. The adaptive FEM algorithm.

1	Given a tolerance $ϵ>0$ and mesh refinement threshold $\tau \in (0,1)$;
2	Choose $d_j$ and $\sigma$ such that ;
3	Construct an initial tetrahedral partition ${ℳ}_h$ over $D$ and compute error estimators;
4	While ${ϵ}_h>ϵ$ do
5	choose ${\widehat{ℳ}}_h\subset {ℳ}_h$ according to the strategy ${\eta }_{{\widehat{ℳ}}_h}>\tau {\eta }_{{ℳ}_h}$;
6	refine all the elements in ${\widehat{ℳ}}_h$ and obtain a new mesh denoted still by ${ℳ}_h$;
7	solve the discrete problem (4.1) on the new mesh ${ℳ}_h$;
8	compute the corresponding error estimators;
9	End while.

Equations

$$\Delta p^{\mathrm{inc}}+{\kappa }^2{p}^{\mathrm{inc}}=0\text{in}{\Omega }_e,$$

$$\Delta p+{\kappa }^{2}p=0\text{in}{\Omega }_e.$$

$$p=p^{\mathrm{inc}}+p^{\mathrm{sc}}\text{in}{\Omega }_e,$$

$$\underset{r\to \infty }{\lim }r({\partial }_r{p}^{\mathrm{sc}}-\mathrm{i}\kappa p^{\mathrm{sc}})=0,r=|\mathbf{x}|.$$

$$\nabla \cdot \mathbf{\sigma }(\mathbf{u})+{\omega }^2\mathbf{u}=0\text{in}{\Omega }_s,$$

$$\mathbf{\sigma }(\mathbf{u})=2\mu \mathbf{ϵ}(\mathbf{u})+\lambda \mathrm{tr}(\mathbf{ϵ}(\mathbf{u}))I,\mathbf{ϵ}(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\top }).$$

$$\nabla \mathbf{u}=\left[\begin{array}{ccc}{\partial }_{x_1}{u}_1& {\partial }_{x_2}{u}_1& {\partial }_{x_3}{u}_1\\ {\partial }_{x_1}{u}_2& {\partial }_{x_2}{u}_2& {\partial }_{x_3}{u}_2\\ {\partial }_{x_1}{u}_3& {\partial }_{x_2}{u}_3& {\partial }_{x_3}{u}_3\end{array}\right].$$

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

$${\partial }_{{\mathbf{n}}_1}p={\rho }_a{\omega }^2{\mathbf{n}}_1\cdot \mathbf{u}\text{on}{\Gamma }_s,$$

$$-p{\mathbf{n}}_1=\mathbf{\sigma }(\mathbf{u})\cdot {\mathbf{n}}_1\text{on}{\Gamma }_s,$$

$$\{\begin{array}{ll}\Delta p+{\kappa }^{2}p=0,p=p^{\mathrm{inc}}+p^{\mathrm{sc}}& \text{in}{\Omega }_e,\\ \nabla \cdot \mathbf{\sigma }(\mathbf{u})+{\omega }^2\mathbf{u}=0& \text{in}{\Omega }_s,\\ {\partial }_{{\mathbf{n}}_1}p={\rho }_a{\omega }^2{\mathbf{n}}_1\cdot \mathbf{u},-p{\mathbf{n}}_1=\mathbf{\sigma }(\mathbf{u})\cdot {\mathbf{n}}_1& \text{on}{\Gamma }_s,\\ {\partial }_r{p}^{\mathrm{sc}}-\mathrm{i}\kappa p^{\mathrm{sc}}=o(r^{-1})& \text{as}r\to \infty .\end{array}$$

$$\mathcal{T}v={\partial }_{{\mathbf{n}}_2}u\text{on}\partial B,$$

$$\{\begin{array}{ll}\Delta u+{\kappa }^{2}u=0& \text{in}{\mathbb{R}}^3\setminus \bar{B},\\ u=v& \text{on}\partial B,\\ {\partial }_{r}u-\mathrm{i}\kappa u=o(r^{-1})& \text{as}r\to \infty .\end{array}$$

$$\{\begin{array}{ll}\Delta p+{\kappa }^{2}p=0& \text{in}{\Omega }_a,\\ \nabla \cdot \mathbf{\sigma }(\mathbf{u})+{\omega }^2\mathbf{u}=0& \text{in}{\Omega }_s,\\ {\partial }_{{\mathbf{n}}_1}p={\rho }_a{\omega }^2{\mathbf{n}}_1\cdot \mathbf{u},-p{\mathbf{n}}_1=\mathbf{\sigma }(\mathbf{u})\cdot {\mathbf{n}}_1& \text{on}{\Gamma }_s,\\ {\partial }_{{\mathbf{n}}_2}p=\mathcal{T}p+f& \text{on}\partial B,\end{array}$$

$$\mathbf{X}:=H^1({\Omega }_a)\times H^1({\Omega }_s{)}^3=\{\mathbf{\Phi }=(p,\mathbf{u}):p\in H^1({\Omega }_a),\mathbf{u}\in H^1({\Omega }_s{)}^3\},$$

$$(\mathbf{\Phi },\mathbf{\Psi }{)}_{\mathbf{X}}:={\int }_{{\Omega }_a}\left(\nabla p\cdot \nabla \bar{q}+p\bar{q}\right)\mathrm{d}\mathbf{x}+{\int }_{{\Omega }_s}\left(\nabla \mathbf{u}:\nabla \bar{\mathbf{v}}+\mathbf{u}\cdot \bar{\mathbf{v}}\right)\mathrm{d}\mathbf{x}$$

$$a(p,\mathbf{u};q,\mathbf{v})=$$

$$+{\int }_{{\Omega }_s}\left(\mathbf{\sigma }(\mathbf{u}):\nabla \bar{\mathbf{v}}-{\omega }^2\mathbf{u}\cdot \bar{\mathbf{v}}\right)\mathrm{d}\mathbf{x}+{\int }_{{\Gamma }_s}(p{\mathbf{n}}_1)\cdot \bar{\mathbf{v}}\mathrm{d}s.$$

$$a(p,\mathbf{u};q,\mathbf{v})={\int }_{\partial B}f\bar{q}\mathrm{d}s,\forall \mathbf{\Psi }=(q,\mathbf{v})\in \mathbf{X}.$$

$$\underset{0\ne (q,\mathbf{v})\in \mathbf{X}}{\sup }\frac{|a(p,\mathbf{u};q,\mathbf{v})|}{\Vert (q,\mathbf{v}){\Vert }_{\mathbf{X}}}\ge {\gamma }_0\Vert (p,\mathbf{u}){\Vert }_{\mathbf{X}},\forall (p,\mathbf{u})\in \mathbf{X}.$$

$${\sigma }_j(t)=0\text{for}|t|<L_j\text{and}{\sigma }_j(t)={\sigma }_0{\left(\frac{|t|-L_j}{d_j}\right)}^m\text{otherwise}.$$

$${\tilde{x}}_j={\int }_0^{x_j}{\alpha }_j(\tau )\mathrm{d}\tau ,1\le j\le 3.$$

$$\tilde{p}(\mathbf{x})=\{\begin{array}{ll}{p}^{\mathrm{inc}}(\mathbf{x})+(p(\tilde{\mathbf{x}})-p^{\mathrm{inc}}(\tilde{\mathbf{x}})),& \mathbf{x}\in {\Omega }_{\mathrm{PML}},\\ p(\tilde{\mathbf{x}}),& \mathbf{x}\in {\Omega }_a.\end{array}$$

$$\mathcal{L}(\tilde{p}-p^{\mathrm{inc}})=0\text{in}\Omega ,$$

$$\mathcal{L}p=\nabla \cdot (A\nabla p)+{\kappa }^{2}bp,$$

$$A=\mathrm{diag}\left(\frac{{\alpha }_2{\alpha }_3}{{\alpha }_1},\frac{{\alpha }_1{\alpha }_3}{{\alpha }_2},\frac{{\alpha }_1{\alpha }_2}{{\alpha }_3}\right),b={\alpha }_1{\alpha }_2{\alpha }_3.$$

$$\{\begin{array}{ll}\mathcal{L}\hat{p}=g& \text{in}\Omega ,\\ \nabla \cdot \mathbf{\sigma }(\hat{\mathbf{u}})+{\omega }^2\hat{\mathbf{u}}=0& \text{in}{\Omega }_s,\\ {\partial }_{{\mathbf{n}}_1}\hat{p}={\rho }_a{\omega }^2{\mathbf{n}}_1\cdot \hat{\mathbf{u}},-\hat{p}{\mathbf{n}}_1=\mathbf{\sigma }(\hat{\mathbf{u}})\cdot {\mathbf{n}}_1& \text{on}{\Gamma }_s,\\ \hat{p}=p^{\mathrm{inc}}& \text{on}\Gamma ,\end{array}$$

$$g=\{\begin{array}{ll}\mathcal{L}{p}^{\mathrm{inc}}& \text{in}{\Omega }_{\mathrm{PML}},\\ 0& \text{in}{\Omega }_a.\end{array}$$

$$\mathbf{Y}:=H^1(\Omega )\times H^1({\Omega }_s{)}^3=\{\mathbf{\Phi }=(p,\mathbf{u}):p\in H^1(\Omega ),\mathbf{u}\in H^1({\Omega }_s{)}^3\},$$

$$(\mathbf{\Phi },\mathbf{\Psi }{)}_{\mathbf{Y}}:={\int }_{\Omega }\left(\nabla p\cdot \nabla \bar{q}+p\bar{q}\right)\mathrm{d}\mathbf{x}+{\int }_{{\Omega }_s}\left(\nabla \mathbf{u}:\nabla \bar{\mathbf{v}}+\mathbf{u}\cdot \bar{\mathbf{v}}\right)\mathrm{d}\mathbf{x}$$

$$b(\hat{p},\hat{\mathbf{u}};q,\mathbf{v})=-{\int }_{\Omega }g\bar{q}\mathrm{d}\mathbf{x},\forall (q,\mathbf{v})\in {\mathbf{Y}}_0,$$

Full text

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.

1. Introduction

Consider the incidence of a time-harmonic acoustic wave onto a bounded, penetrable, and isotropic elastic solid, which is immersed in a homogeneous and compressible air or fluid. Due to the interaction between the incident wave and the solid obstacle, an elastic wave is excited inside the solid region, while the acoustic incident wave is scattered in the air/fluid region. This scattering phenomenon leads to an air/fluid-solid interaction problem. The surface of the elastic solid divides the whole three-dimensional space into a bounded interior domain and an open exterior domain where the elastic wave and the acoustic wave occupies, respectively. The two waves are coupled together on the surface via the interface conditions: continuity of the normal component of velocity and the continuity of traction. The acoustic-elastic interaction problems have received ever-increasing attention due to their significant applications in geophysics and seismology [22, 23]. These problems have been examined mathematically by using either variational method [18, 19] or boundary integral equation method [28, 24]. Many computational approaches have also been developed to numerically solve these problems such as boundary element method [17, 31] and coupling of finite and boundary element methods [16].

Since the work by Bérenger [4], the perfectly matched layer (PML) technique has been extensively studied and widely used to simulate various wave propagation problems, which include acoustic waves [5, 12, 21, 27, 32], elastic waves [6, 11, 13, 20, 26], and electromagnetic waves [3, 15]. The PML 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. It has been proven to be an effective approach to truncated open domains in the wave computation. Combined with the PML technique, the adaptive finite element method (FEM) has recently been developed to solve the diffraction grating problems [2, 8, 25] and the obstacle scattering problems [7, 9, 10]. Despite the large number of work done so far, they were concerned with a single wave propagation problem, i.e., either an acoustic wave, or an elastic wave, or an electromagnetic wave. It is very rare to study rigorously the PML problem for the interaction of multiple waves.

This paper aims to investigate the adaptive finite element PML method for solving the acoustic-elastic interaction problem. An exact transparent boundary condition (TBC) is developed to reduce the problem equivalently into a boundary value problem in a bounded domain. The PML technique is adopted to truncated the unbounded physical domain into a bounded computational domain. The variational approach is taken to incorporate naturally the interface conditions which couple the two waves. The well-posedness and exponential convergence of the solution are established for the truncated PML problem by using a PML equivalent TBC. The proofs rely on the error estimate between the two transparent boundary operators. To effciently resolve the solution with possible singularities, the a posteriori error estimate based adaptive FEM is developed to solve the truncated PML problem. The error estimate consists of the PML error and the finite element discretization error, and provides a theoretical basis for the mesh refinement. Numerical experiments are reported to show the competitive behavior of the proposed method.

The paper is organized as follows. In section 2, we introduce the model equations for the acoustic-elastic interaction problem. In section 3, we present the PML formulation and prove the well-posedness and convergence of the solution for the truncated PML problem. In section 4, we discuss the numerical implementation and show some numerical experiments. The paper is concluded with some general remarks in section 5.

2. Problem formulation

In this section, we introduce the model equations for acoustic and elastic waves, and present an interface problem for the acoustic-elastic interaction. In addition, an exact transparent boundary condition is introduced to reformulate the scattering problem into an boundary value problem in an bounded domain.

2.1. Problem geometry

Consider an acoustic plane wave incident on a bounded elastic solid which is immersed in a homogeneous compressible air/fluid in three dimensions. The problem geometry is shown in Figure 1. Due to the wave interaction, an elastic wave is induced inside the solid region, while the scattered acoustic wave is generated in the open air/fluid region. The wave propagation described above leads to an air/fluid-solid interaction problem. The surface of the solid divides the whole three-dimensional space into the interior domain and the exterior domain, where the elastic wave and the acoustic wave occupies, respectively. Let the solid $\Omega_{s}\subset\mathbb{R}^{3}$ be a bounded domain with a Lipschitz boundary $\Gamma_{s}$. The exterior domain $\Omega_{e}=\mathbb{R}^{3}\setminus\bar{\Omega}_{s}$ is assumed to be connected and filled with a homogeneous, compressible, and inviscid air/fluid with a constant density $\rho_{a}>0$. Denote by $B=\{\boldsymbol{x}=(x_{1},x_{2},x_{3})^{\top}\in\mathbb{R}^{3}:|x_{j}|<L_{j},j=1,2,3\}$ the rectangular box with the boundary $\partial B$, where $L_{j}$ are sufficiently large such that $\bar{\Omega}_{s}\subset B$. Define $\Omega_{a}=B\setminus\bar{\Omega}_{s}$. Let $\boldsymbol{n}_{1}$ be the unit normal vector on $\Gamma_{s}$ directed from $\Omega_{s}$ into $\Omega_{e}$, and let $\boldsymbol{n}_{2}$ be the unit outward normal vector on $\partial B$.

2.2. Wave equations

Let the elastic solid be impinged by a time-harmonic sound wave $p^{\rm inc}$, which satisfies the three-dimensional Helmholtz equation:

[TABLE]

where $\kappa=\omega/c$ is the wavenumber, $\omega>0$ is the angular frequency, and $c$ is the speed of sound in the air/fluid. The total acoustic wave field $p$ also satisfies the Helmholtz equation:

[TABLE]

The total field $p$ consists of the incident field $p^{\rm inc}$ and the scattered field $p^{\rm sc}$:

[TABLE]

where scattered field $p^{\rm sc}$ is required to satisfy the Sommerfeld radiation condition:

[TABLE]

The time-harmonic elastic wave satisfies the three-dimensional Navier equation:

[TABLE]

where $\boldsymbol{u}=(u_{1},u_{2},u_{3})^{\top}$ is the displacement of the elastic wave, and the stress tensor $\boldsymbol{\sigma}(\boldsymbol{u})$ is given by the generalized Hook law:

[TABLE]

Here $\mu(\boldsymbol{x})\in L^{\infty}(\Omega_{s}),\lambda(\boldsymbol{x})\in L^{\infty}(\Omega_{s})$ are the Lamé parameters satisfying $\mu>0,\lambda>0$, and $\nabla\boldsymbol{u}$ is the displacement gradient tensor given by

[TABLE]

Substituting (2.3) into (2.2) yields

[TABLE]

2.3. Interface conditions

To couple the acoustic wave equation and the elastic wave equation, the kinematic interface condition is imposed to ensure the continuity of the normal component of the velocity:

[TABLE]

In addition, the dynamic interface condition is required to ensure the continuity of traction:

[TABLE]

where $\boldsymbol{\sigma}(\boldsymbol{u})\cdot\boldsymbol{n}_{1}$ denotes the matrix-vector multiplication.

2.4. Acoustic-elastic interaction problem

The acoustic-elastic interaction problem can be formulated into the following coupled boundary value problem: Given $p^{\rm inc}$, to find $(p,\boldsymbol{u})$ such that

[TABLE]

We refer to [28] for the discussion on the well-posedness of the boundary value problem (2.7). From now on, we assume that the acoustic-elastic interaction problem has a unique solution.

2.5. Transparent boundary condition

Given $v\in H^{1/2}(\partial B)$, we define the Dirichlet-to-Neumann (DtN) operator $\mathscr{T}:H^{1/2}(\partial B)\to H^{-1/2}(\partial B)$ as follows:

[TABLE]

where $u$ is the solution of the exterior Dirichlet problem of the Helmholtz equation:

[TABLE]

It is well-known that the exterior problem (2.8) has a unique solution $u\in H^{1}_{loc}(\mathbb{R}^{3}\setminus\bar{B})$ (cf., e.g., [14]). Thus the DtN operator $\mathscr{T}:H^{1/2}(\partial B)\to H^{-1/2}(\partial B)$ is well-defined and is a bounded linear operator.

Using the DtN operator $\mathscr{T}$, we reformulate the boundary value problem (2.7) from the open domain into the bounded domain: Given $p^{\rm inc}$, to find $(p,\boldsymbol{u})$ such that

[TABLE]

where $f=\partial_{\boldsymbol{n}_{2}}p^{\rm inc}-\mathscr{T}p^{\rm inc}$.

To study the well-posedness of (2.9), we define

[TABLE]

which is endowed with the inner product:

[TABLE]

for any $\boldsymbol{\Phi}=(p,\boldsymbol{u})$ and $\boldsymbol{\Psi}=(q,\boldsymbol{v})$, where $A:B={\rm tr}(AB^{\top})$ is the Frobenius inner product of square matrices $A$ and $B$. Clearly, $\|\cdot\|_{\boldsymbol{X}}=\sqrt{(\cdot,\cdot)_{\boldsymbol{X}}}$ is a norm on $\boldsymbol{X}$.

Let $a:\boldsymbol{X}\times\boldsymbol{X}\to\mathbb{C}$ be the sesquilinear form:

[TABLE]

The acoustic-elastic interaction problem (2.9) is equivalent to the following weak formulation: Find $\boldsymbol{\Phi}=(p,\boldsymbol{u})\in\boldsymbol{X}$ such that

[TABLE]

Since we assume that the variational problem (2.11) has a unique weak solution $(p,\boldsymbol{u})\in\boldsymbol{X}$, the general theory in Babuška and Aziz [1, Chap. 5] implies that there exists a constant $\gamma_{0}$ such that the following inf-sup condition is satisfied

[TABLE]

3. The PML problem

In this section, we introduce the PML formulation for the acoustic-elastic interaction problem and establish its well-posedness. 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 $\Omega_{a}$ is surrounded by a PML layer of thickness $d_{j}$ which is denoted as $\Omega_{\rm PML}$. Define $\Omega:=\Omega_{a}\cup\partial B\cup\Omega_{\rm PML}$. Let $\alpha_{j}(t)=1+{\rm i}\sigma_{j}(t)$ be the PML function which is continuous and satisfies

[TABLE]

Here $\sigma_{0}>0$ is a constant and $m$ is an integer. Following [15], we introduce the PML by the complex coordinate stretching:

[TABLE]

Let $\tilde{\boldsymbol{x}}=(\tilde{x}_{1},\tilde{x}_{2},\tilde{x}_{3})$. Introduce the new function:

[TABLE]

It is clear to note that $\tilde{p}({\boldsymbol{x}})=p(\boldsymbol{x})$ in $\Omega_{a}$ since $\tilde{\boldsymbol{x}}=\boldsymbol{x}$ in $\Omega_{a}$. It can be verified from (2.1) and (3.1) that $\tilde{p}$ satisfies

[TABLE]

where the PML differential operator is defined by

[TABLE]

where

[TABLE]

It can be verified from (2.1) and (3.2) that the outgoing wave $\tilde{p}(\boldsymbol{x})-p^{\rm inc}(\boldsymbol{x})$ in $\Omega_{\rm PML}$ decays exponentially. Therefore, the homogeneous Dirichlet boundary condition can be imposed on $\Gamma:=\partial\Omega_{\rm PML}\setminus\partial B$ to truncate the PML problem. We arrive at the following truncated PML problem: Find $(\hat{p},\hat{\boldsymbol{u}})$ such that

[TABLE]

where

[TABLE]

Define

[TABLE]

which is endowed with the inner product

[TABLE]

for any $\boldsymbol{\Phi}=(p,\boldsymbol{u})$ and $\boldsymbol{\Psi}=(q,\boldsymbol{v})$. Obviously, $\|\cdot\|_{\boldsymbol{Y}}=\sqrt{(\cdot,\cdot)_{\boldsymbol{Y}}}$ is a norm on $\boldsymbol{Y}$.

The weak formulation of the truncated PML problem (3.3) reads as follows: Find $(\hat{p},\hat{\boldsymbol{u}})\in\boldsymbol{Y}$ such that $\hat{p}=p^{\rm inc}$ on $\Gamma$ and

[TABLE]

where $\boldsymbol{Y}_{0}=\{\boldsymbol{\Phi}=(p,\boldsymbol{u})\in\boldsymbol{Y}:p=0~{}\text{on}~{}\Gamma\}$, and the sesquilinear form $b:\boldsymbol{Y}\times\boldsymbol{Y}\to\mathbb{C}$ is defined by

[TABLE]

We will reformulate the variational problem (3.4) imposed in the domain $\Omega\cup\bar{\Omega}_{s}$ into an equivalent variational formulation in the domain $B=\Omega_{a}\cup\bar{\Omega}_{s}$, 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

We start by introducing the approximate DtN operator $\mathscr{T}^{\rm PML}:H^{1/2}(\partial B)\to H^{-1/2}(\partial B)$ associated with the PML problem.

Given $\psi\in H^{1/2}(\partial B)$, let $\mathscr{T}^{\rm PML}\psi=\partial_{\boldsymbol{n}_{2}}\phi$ on $\partial B$, where $\phi\in H^{1}(\Omega_{\rm PML})$ is the solution of the following boundary value problem in the PML layer:

[TABLE]

The PML problem (3.3) can be reduced to the following boundary value problem: Find $(p^{\rm PML},\boldsymbol{u}^{\rm PML})$ such that

[TABLE]

where $f^{\rm PML}=\partial_{\boldsymbol{n}_{2}}p^{\rm inc}-{\mathscr{T}}^{\rm PML}p^{\rm inc}$.

The weak formulation of (3.5) is to find $(p^{\rm PML},\boldsymbol{u}^{\rm PML})\in\boldsymbol{X}$ such that

[TABLE]

where the sesquilinear form $a^{\rm PML}:\boldsymbol{X}\times\boldsymbol{X}\to\mathbb{C}$ is defined by

[TABLE]

The following lemma establishes the relationship between the variational problem (3.6) and the weak formulation (3.4). The proof is straightforward based on our constructions of the transparent boundary conditions for the PML problem. The details of the proof is omitted for simplicity.

Lemma 3.1**.**

Any solution $\hat{p}$ of the variational problem (3.4) restricted to $\Omega_{a}$ is a solution of the variational (3.6); conversely, any solution $p^{\rm PML}$ of the variational problem (3.6) can be uniquely extended to the whole domain to be a solution $\hat{p}$ of the variational problem (3.4) in $\Omega$.

3.3. Convergence of the PML solution

Now we turn to estimating the error between $(p^{\rm PML},\boldsymbol{u}^{\rm PML})$ and $(p,\boldsymbol{u})$. The key is to estimate the error of the boundary operators $\mathscr{T}^{\rm PML}$ and $\mathscr{T}$.

Lemma 3.2**.**

For any $p,q\in H^{1}(\Omega_{a})$, there exists a constant $C>0$ such that

[TABLE]

where $L=\max_{1\leq j\leq 3}L_{j},\alpha_{0}=\max_{\boldsymbol{x}\in\Gamma}(|\alpha_{1}(x_{1})|,|\alpha_{2}(x_{2})|,|\alpha_{3}(x_{3})|)$,

[TABLE]

and $\sigma>0$ is a sufficiently large constant such that $\gamma_{1}\sigma\geq 1$.

Proof.

The proof can follow similar arguments as that in [5, Theorem 3.8]. For the sake of simplicity, we do not elaborate on the details here. ∎

Theorem 3.3**.**

Let $\gamma_{0}$ be the constant in the inf-sup condition (2.12). If

[TABLE]

then the PML variational problem (3.6) has a unique weak solution $(p^{\rm PML},\boldsymbol{u}^{\rm PML})$, which satisfies the error estimate

[TABLE]

where $(p,\boldsymbol{u})$ is the unique weak solution of the variational problem (2.11).

Proof.

It suffices to show the coercivity of the sesquilinear form $a^{\rm PML}$ defined in (3.2) in order to prove the unique solvability of the weak problem (3.6). Using Lemma 3.2, and the assumption $\gamma_{2}<\gamma_{0}$, we get for any $(p,\boldsymbol{u}),(q,\boldsymbol{v})$ in $\boldsymbol{X}$ that

[TABLE]

It remains to show the error estimate (3.8). It follows from (3.6)–(3.2) that

[TABLE]

which completes the proof upon using Lemma 3.2 and the trace theorem. ∎

4. Finite element approximation

In this section we introduce the finite element approximations of the PML problem (3.4).

4.1. Error representation formula

Let $\mathcal{M}_{h}$ be a regular tetrahedral partition of the domain $D=\Omega\cup\Gamma_{s}\cup\Omega_{s}=\{\boldsymbol{x}\in\mathbb{R}^{3}:|x_{j}|<L_{j}+d_{j},1\leq j\leq 3\}$ such that $\mathcal{M}_{h}|_{\Omega}$ and $\mathcal{M}_{h}|_{\Omega_{s}}$ are also regular tetrahedral partitions of $\Omega$ and $\Omega_{s}$, respectively. Let $V_{h}\subset H^{1}(\Omega)$ and $\boldsymbol{U}_{h}\subset H^{1}(\Omega_{s})^{3}$ be the conforming linear finite element space over $\Omega$ and $\Omega_{s}$, respectively, and

[TABLE]

The finite element approximation to the PML problem (3.4) reads as follows: Find $(p_{h},\boldsymbol{u}_{h})\in V_{h}\times\boldsymbol{U}_{h}$ such that $p_{h}=I_{h}p^{\rm inc}$ on $\Gamma$ and

[TABLE]

For any $\varphi\in H^{1}(\Omega_{a})$, let $\tilde{\varphi}$ be its extension in $\Omega_{\rm PML}$ such that

[TABLE]

Introduce the sesquilinear form $c:H^{1}(\Omega_{\rm PML})\times H^{1}(\Omega_{\rm PML})\to\mathbb{C}$ as follows:

[TABLE]

The weak formulation for (4.2)–(4.3) is: Given $\varphi\in H^{1/2}(\partial B)$, find $\tilde{\varphi}\in H^{1}(\Omega_{\rm PML})$ such that $\tilde{\varphi}=0$ on $\Gamma$, $\tilde{\varphi}=\varphi$ on $\partial B$, and

[TABLE]

In this paper we will not elaborate on the well-posedness of (4.4) and simply make the following assumption: There exists a unique solution to the boundary value problem (4.4) in the PML layer.

In order to obtain a constant independent of PML parameter $\sigma$ in the inf-sup condition, we define

[TABLE]

By using the general theory in [1, Chap. 5], we know that there exists a constant $\hat{C}>0$ such that

[TABLE]

The constant $\hat{C}$ depends on the domain $\Omega_{\rm PML}$ and the wave number $\kappa$.

Lemma 4.1** (Estimates for the extension).**

For any $\varphi\in H^{1}(\Omega_{a})$, which is extended to be a function $\tilde{\varphi}\in H^{1}(\Omega)$ according to (4.2)–(4.3). Then there exists a constant $C>0$ independent of $\kappa$ and $\sigma$ such that

[TABLE]

where $\boldsymbol{n}_{3}$ is the unit outward normal vector on $\Gamma$.

Proof.

For any $\zeta\in H^{1}(\Omega_{\rm PML})$ such that $\zeta=\varphi$ on $\partial B$ and $\zeta=0$ on $\Gamma$. By the inf-sup condition in (4.5) and using (4.4), we know that

[TABLE]

By Cauchy–Schwarz inequality

[TABLE]

Noting

[TABLE]

using the triangle inequality and the trace inequality, we conclude that

[TABLE]

which shows the first estimate in the theorem by using the definition of $|||\cdot|||_{\Omega_{\rm PML}}$.

Next, for any $\psi\in H^{1}(\Omega_{\rm PML})$ such that $\psi=0$ on $\partial B$, using (4.2) and the integration by parts, we obtain

[TABLE]

It follows from the Cauchy–Schwarz inequality and (4.8) that

[TABLE]

which completes the proof after using the trace inequality. ∎

Lemma 4.2** (Error representation formula).**

For any $\varphi\in H^{1}(\Omega_{a})$, which is extended to be a function $\tilde{\varphi}\in H^{1}(\Omega)$ according to (4.2)–(4.3), and $\varphi_{h}\in V_{\Gamma,h}$, we have

[TABLE]

Proof.

First by (2.5), (2.11), (3.6), and (3.2), we have

[TABLE]

Using (4.1) yields

[TABLE]

Recalling that $\boldsymbol{n}_{2}$ is the unit outer normal to $\partial B$ which points outside $B$ and $\boldsymbol{n}_{3}$ is the unit outer normal vector on $\Gamma$ directed outside $\Omega_{\rm PML}$, we deduce that

[TABLE]

where we have used (4.2)–(4.3), the definition of $\mathscr{T}^{\rm PML}$, and the identity (c.f., [9, Lemma 5.1])

[TABLE]

By (3.4), (4.1), and (4.1)–(4.1),

[TABLE]

which completes the proof. ∎

4.2. A posteriori error analysis

For any $K\in\mathcal{M}_{h}$, we denote by $h_{K}$ its diameter. Let $\mathcal{B}_{h}$ denote the set of all sides that do not lie on $\Gamma$. For any $e\in\mathcal{B}_{h}$, $h_{e}$ stands for its length. For any $K\in\mathcal{M}_{h}$, we introduce the residual:

[TABLE]

For any interior side $e\in\mathcal{B}_{h}$ not lying on the interface $\Gamma_{s}$ which is the common side of $K_{1},K_{2}\in\mathcal{M}_{h}$, we define the jump residual across $e$:

[TABLE]

where we have used the notation that the unit normal vector $\boldsymbol{\nu}$ on $e$ points from $K_{2}$ to $K_{1}$. If $e$ lies on the interface $\Gamma_{s}$, then we define the jump residual as

[TABLE]

For any $K\in\mathcal{M}_{h}$, we define the local error estimator $\eta_{K}$ as

[TABLE]

Theorem 4.3**.**

There exists a constant $C>0$ depending only on $\gamma_{1}$ and the minimum angle of the mesh $\mathcal{M}_{h}$ such that the following a posterior error estimate holds

[TABLE]

Proof.

Let $\Pi_{h}:H^{1}_{\Gamma}(\Omega)\to V_{\Gamma,h}$ and $\boldsymbol{\Pi}_{h}:H^{1}(\Omega_{s})^{3}\to\boldsymbol{U}_{h}$ be Scott–Zhang [30] interpolation operators satisfying the following interpolation estimates: For any $\varphi\in H^{1}(\Omega)$ and $\boldsymbol{v}\in H^{1}(\Omega_{s})^{3}$,

[TABLE]

and

[TABLE]

where $\tilde{K}$ and $\tilde{e}$ are the union of all elements in $\mathcal{M}_{h}$ having a non-empty intersection with $K\in\mathcal{M}_{h}$ and the side $e$, respectively.

Taking $\tilde{\varphi}_{h}=\Pi_{h}\tilde{\varphi}\in V_{\Gamma,h}$ and $\boldsymbol{v}_{h}=\boldsymbol{\Pi}_{h}\boldsymbol{v}\in\boldsymbol{U}_{h}$ in the error representation formula (4.2), we get

[TABLE]

It follows from the integration by parts and (4.13)–(4.15) that

[TABLE]

By (4.16)–(4.17) and the estimate (4.6), we have

[TABLE]

By Lemma 3.2, we have

[TABLE]

It follows from (4.7) that

[TABLE]

The proof is completed by using the above estimates in (4.2) and the inf-sup condition (2.12). ∎

5. Numerical experiments

According to the discussion in section 4, we choose the PML medium property as the power function and need to specify the thickness $d_{j}$ of the layers and the medium parameter $\sigma$. It is clear to note from Theorem 4.3 that the a posteriori error estimate consists of two parts: the PML error $\epsilon_{\rm PML}$ and the finite element discretization error $\epsilon_{\rm FEM}$, where

[TABLE]

In our implementation, we first choose $d_{j}$ and $\sigma$ such that $\alpha_{0}^{3}(1+\kappa L)^{3}e^{-\gamma_{1}\kappa\sigma}\leq 10^{-8}$, which makes the PML error (5.2) negligible compared with the finite element discretization error (5.1). Once the PML region and the medium property are fixed, we use the standard finite element adaptive strategy to modify the mesh according to the a posteriori error estimate. For any $K\in\mathcal{M}_{h}$, we define the local a posteriori error estimator

[TABLE]

The adaptive FEM algorithm is summarized in Table 1.

In the following, we present two examples to demonstrate the competitive numerical performance of the proposed algorithm. The first-order linear element is used for solving the problem. Our implementation is based on parallel hierarchical grid (PHG) [29], which is a toolbox for developing parallel adaptive finite element programs on unstructured tetrahedral meshes. The linear system resulted from finite element discretization is solved by the PCG solver.

Example 1. We consider a problem with an exact solution. We set the elastic region $\Omega_{s}:=B(0,0.2)$ and the acoustic region $\Omega_{a}:=B(0,0.5)\setminus\bar{\Omega}_{s}$, where $B(0,R):=\{\boldsymbol{x}\in\mathbb{R}^{3}:|\boldsymbol{x}|<R\}$ denotes the ball with radius $R>0$ and centering at the origin. Let

[TABLE]

where $\boldsymbol{x}_{0}=(1,0,0)^{\top}$. The parameters are chosen as $\kappa=1$, $\omega=1$, $\lambda=0.5$, $\mu=0.25$, and $\rho_{a}=1$ such that

[TABLE]

First it is easy to verify that

[TABLE]

When $\mu$ and $\lambda$ are constants, the Navier equation (2.4) reduces to

[TABLE]

Using (5.3) and (5.4), we have from a straightforward calculation that

[TABLE]

which shows that $\boldsymbol{u}=\omega^{2}\nabla p$ satisfies (5.5) in $\Omega_{s}$. It can be verified that the interface conditions (2.5)–(2.6) are also satisfied by letting $\rho_{a}=1$.

Let $q=p|_{\partial B(0,0.5)}$ and consider the following acoustic-elastic interaction problem with the Dirichlet boundary condition:

[TABLE]

We may test the adaptive FEM algorithm by solving the above boundary value problem.

Figure 3 displays the errors of $p$ and $\boldsymbol{u}$ against the number of nodal points $N_{p}$ in $B(0,0.5)\setminus\bar{B}(0,0.2)$ and $N_{\boldsymbol{u}}$ in $B(0,0.2)$, respectively. It clearly shows that the adaptive FEM yields quasi-optimal convergence rates, i.e.,

[TABLE]

and

[TABLE]

where $\eta_{p,h}$ and $\eta_{\boldsymbol{u},h}$ are the a posterior error estimators for $p$ and $\boldsymbol{u}$, respectively. Figure 4 plots the adaptive mesh of $\Omega_{a}$ for solving $p_{h}$ and Figure 5 plots the mesh on a cross section of the domain $\Omega_{a}$ on the $xz$-plane. Figure 6 plots the adaptive mesh of $\Omega_{s}$ for solving $\boldsymbol{u}_{h}$ and Figure 7 plot the mesh on the cross section of the domain $\Omega_{s}$ on the $xz$-plane.

Example 2. This example concerns the scattering of the incident plane wave

[TABLE]

The Dirichlet boundary condition on the PML layer outer boundary $\Gamma$ is set by $p=p^{\rm inc}$. We choose $\kappa=2$, $\omega=2\pi$, $\lambda=1$, $\mu=2$, and $\rho_{a}=1$. Let the elastic region and the acoustic region be $\Omega_{s}=B_{1}\backslash\bar{B}_{0}$ and $\Omega_{a}=B_{2}\backslash\bar{\Omega}_{s}$, respectively. Here $B_{0}=(-0.1,0.1)\times(-0.1,0.1)\times(-0.2,0.0),B_{1}=(-0.2,0.2)\times(-0.2,0.2)\times(-0.2,0.2)$, and $B_{2}=[-0.6,0.6]\times[-0.6,0.6]\times[-0.6,0.6]$. The PML domain is $\Omega_{\rm PML}=(0,1)\times(0,1)\times(0,1)\setminus\bar{B}_{2}$, i.e., the thickness of the PML layer is 0.4 in each direction. In this example, the elastic solid is a rectangular box with a small rectuangular dent on the surface. The solutions of $p$ and $\boldsymbol{u}$ may have singularities around the corners of the dent. We choose $\sigma=16$ and $m=2$ for the medium property to ensure the PML error is negligible compared to the finite element error.

For this example, we set the numerical solution on the very fine mesh to be a reference solution since there is no analytic solution. Figure 8 shows the errors of $p$ and $\boldsymbol{u}$ against the number of nodal points $N_{p}$ and $N_{\boldsymbol{u}}$. It is clear to note that the FEM algorithm yields a quasi-optimal convergence rate. The surface plots of the amplitude of the fields are shown as follows: Figure 9 shows the real part of $p_{h}$ for the cross section in $\Omega_{a}$ on the $yz$-plane and Figure 10 shows the real part of $\boldsymbol{u}_{h}$ for the cross section in $\Omega_{s}$ on the $yz$-plane.

6. Concluding remarks

We have studied a variational formulation for the acoustic-elastic interaction problem in $\mathbb{R}^{3}$ and adopted the PML to truncate the unbounded physical domain. The scattering problem is reduced to a boundary value problem by using transparent boundary conditions. We prove that the truncated PML problem has a unique weak solution which converges exponentially to the solution of the original problem by increasing the PML parameters. We incorporate the adaptive mesh refinement with a posteriori error estimate for the finite element method to handle the problem where the solution may have singularities. Numerical results show that the proposed method is effective to solve the acoustic-elastic interaction problem.