An interior penalty discontinuous Galerkin method for the time-domain acoustic-elastic wave interaction problem

TL;DR

This paper develops and analyzes a symmetric interior penalty discontinuous Galerkin method for simulating time-domain acoustic-elastic wave interactions, providing error estimates and numerical validation.

Contribution

It introduces a novel DG method for fluid-solid interaction problems with proven error bounds and numerical validation.

Findings

The method achieves high accuracy in simulations.

Error estimates are rigorously established.

Numerical results confirm the method's effectiveness.

Abstract

In this paper, we consider numerical solutions of a time domain acoustic-elastic wave interaction problem which occurs between a bounded penetrable elastic body and a compressible inviscid fluid. It is also called the fluid-solid interaction problem. Firstly, an artificial boundary is introduced to transform the original transmission problem into a problem in a bounded region, and then its well-posed results are presented. We then employ a symmetric interior penalty discontinuous Galerkin method for the solution of the reduced interaction problem consisting of second-order wave equations. A priori error estimate in the energy norm is presented, and numerical results confirm the accuracy and validity of the proposed method.

Tables

Table 1. Table 1 : Numerical errors when T = 1 𝑇 1 T=1 for Example 1.

h	Error in energy norm	Order	$L^2$ error	Order
0.415	1.3071e+00	-	6.0158e-01	-
0.207	6.9562e-01	0.91	2.1094e-01	1.51
0.104	3.0057e-01	1.21	5.8146e-02	1.86
0.052	1.3364e-01	1.17	1.4833e-02	1.97
0.026	6.3287e-02	1.08	3.7408e-03	1.99

Table 2. Table 2 : Numerical errors when T = 1 𝑇 1 T=1 for Example 2.

h	Error in energy norm	Order	$L_2$ error	Order
0.726	3.9474e-01	-	4.0360e-02	-
0.207	2.5512e-01	0.63	1.4849e-02	1.45
0.104	1.6564e-01	0.65	5.5443e-03	1.41
0.052	1.0484e-01	0.66	2.1281e-03	1.38
0.026	6.6348e-02	0.66	8.2907e-04	1.36

Equations

$${\rho }_1\frac{\partial \mathbf{v}}{\partial t}+\nabla p=0,\frac{\partial p}{\partial t}+c^2{\rho }_1\text{div}\mathbf{v}=0,(\mathbf{x},t)\in {\Omega }^{+}\times J,$$

$$\mathbf{v}=-\nabla \phi ,\text{and}p={\rho }_1\frac{\partial \phi }{\partial t}$$

$$\frac{1}{c^2}\frac{{\partial }^2\phi }{\partial t^2}-\Delta \phi =0,(\mathbf{x},t)\in {\Omega }^{+}\times J.$$

$$\begin{array}{rl}{\rho }_2\frac{{\partial }^2\mathbf{u}}{\partial t^2}-\text{div}(\sigma (\mathbf{u}))& =0,(\mathbf{x},t)\in \Omega \times J\\ \sigma (\mathbf{u})& =\lambda \text{tr}(\epsilon (\mathbf{u}))\mathbf{I}+2\mu \epsilon (\mathbf{u})\\ \epsilon (\mathbf{u})& =\frac{1}{2}(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T).\end{array}$$

$${\rho }_2\frac{{\partial }^2\mathbf{u}}{\partial t^2}-\text{div}(\sigma (\mathbf{u}))$$

$$\frac{1}{c^2}\frac{{\partial }^2\phi }{\partial t^2}-\Delta \phi$$

$$\sigma (\mathbf{u})\mathbf{n}$$

$$\frac{\partial \phi }{\partial \mathbf{n}}$$

$$r^{\frac{1}{2}}\left(\frac{\partial \phi }{\partial \mathbf{n}}+\frac{\partial \phi }{\partial t}\right)$$

$$\mathbf{u}{|}_{t=0}={\frac{\partial \mathbf{u}}{\partial t}|}_{t=0}=0,$$

$$\Vert u{\Vert }_{0,\Omega }={\left({\int }_{\Omega }|u(\mathbf{x},t){|}^{2}d\mathbf{x}\right)}^{1/2}.$$

$$H^s(\Omega )=\left\{D^{\alpha }u\in L^2(\Omega )\text{ for all }|\alpha |\le s\right\}$$

$$\Vert u{\Vert }_{s,\Omega }={\left({\int }_{\Omega }\sum _{|\alpha |<s}|D^{\alpha }u{|}^{2}dx\right)}^{1/2},$$

$$⟨u,v{⟩}_{\Gamma }={\int }_{\Gamma }uv\gamma .$$

$$\Vert \nabla \mathbf{u}{\Vert }_{F(\Omega )}={\left(\sum _{j=1}^2{\int }_{\Omega }|\nabla u_j{|}^{2}d\mathbf{x}\right)}^{1/2}.$$

$$\Vert \nabla \mathbf{u}{\Vert }_{F(\Omega )}^2+\Vert \text{div }\mathbf{u}{\Vert }_{0,\Omega }^2\le C\Vert \mathbf{u}{\Vert }_{1,\Omega }^2,$$

$$\|w(x,t)\|_{L^{q}(J;H^{1}(\Omega))}=\left\{\begin{array}[]{cc}{\left(\int_{J}\|w\|^{q}_{1,\Omega}\right)}^{1/q}&1\leq q<\infty,\\ \text{ess}\sup\limits_{0\leq t\leq T}\|w\|_{1,\Omega}\,dt&q=\infty,\end{array}\right.$$

$$\Vert w(x,t){\Vert }_{H^1(J;H^s(\Omega ))}={\left({\int }_J\left(\Vert w{\Vert }_{s,\Omega }^2+\Vert w_t{\Vert }_{s,\Omega }^{2}dt\right)\right)}^{1/2}.$$

$$\frac{\partial \phi }{\partial \mathbf{n}}=-\frac{\partial \phi }{\partial t},\text{ on }{\Gamma }_R\times J,$$

$${\rho }_2\frac{{\partial }^2\mathbf{u}}{\partial t^2}-\text{div}(\sigma (\mathbf{u}))$$

$$\frac{1}{c^2}\frac{{\partial }^2\phi }{\partial t^2}-\Delta \phi$$

$$\sigma (\mathbf{u})\mathbf{n}$$

$$\frac{\partial \phi }{\partial \mathbf{n}}$$

$$\mathbf{u}{|}_{t=0}={\frac{\partial \mathbf{u}}{\partial t}|}_{t=0}$$

$$\mathbf{u}\in L^2\left(J;{\mathbf{H}}^1(\Omega )\right)\cap H^1\left(J;{\mathbf{L}}^2(\Omega )\right),$$

$$\phi \in L^2\left(J;H^1({\Omega }_R)\right)\cap H^1\left(J;L^2({\Omega }_R)\right),$$

$$\underset{t\in [0,T]}{\max }\left(\Vert {\mathbf{u}}_t{\Vert }_{0,\Omega }+\Vert \nabla \mathbf{u}{\Vert }_{F(\Omega )}+\Vert \nabla \cdot \mathbf{u}{\Vert }_{0,\Omega }\right)$$

$$\le$$

$$+\underset{t\in [0,T]}{\max }{\Vert \nabla {\phi }_t^i\Vert }_{-\frac{1}{2},\Gamma }^2+\Vert {\phi }_{ttt}^i{\Vert }_{L^1(J;H^{-1/2}(\Gamma ))}^2+\Vert \nabla {\phi }_{tt}^i{\Vert }_{L^1(J;H^{-1/2}(\Gamma ))}^2),$$

$$\underset{t\in [0,T]}{\max }\left(\Vert {\phi }_t{\Vert }_{0,{\Omega }_R}+\Vert \nabla \phi {\Vert }_{0,{\Omega }_R}\right)$$

$$\le$$

$$+\underset{t\in [0,T]}{\max }{\Vert \nabla {\phi }_t^i\Vert }_{-\frac{1}{2},\Gamma }^2+\Vert {\phi }_{ttt}^i{\Vert }_{L^1(J;H^{-1/2}(\Gamma ))}^2+\Vert \nabla {\phi }_{tt}^i{\Vert }_{L^1(J;H^{-1/2}(\Gamma ))}^2).$$

$$(\mathbf{u},\phi )\in L^2\left(J;{\mathbf{H}}^1(\Omega )\right)\times L^2\left(J;H^1({\Omega }_R)\right),$$

$${\mathbf{u}}_t\in L^2\left(J;{\mathbf{L}}^2(\Omega )\right),{\phi }_t\in L^2\left(J;L^2({\Omega }_R)\right),{\mathbf{u}}_{tt}\in L^2\left(J;{\mathbf{H}}^{-1}(\Omega )\right),{\phi }_{tt}\in L^2\left(J;H^{-1}({\Omega }_R)\right),$$

$$\frac{{\rho }_2}{{\rho }_1}⟨{\mathbf{u}}_{tt},\mathbf{v}⟩+\frac{1}{c^2}⟨{\phi }_{tt},\varphi ⟩+a(\mathbf{u},\phi ;\mathbf{v},\varphi )=L(\mathbf{v},\varphi ).$$

$$\mathbf{u}{|}_{t=0}={\mathbf{u}}_t{|}_{t=0}=0,\phi {|}_{t=0}={\phi }_t{|}_{t=0}=0.$$

Full text

Keywords:* Fluid-solid interaction problem; Time-domain wave propagation; Second wave equation; Interior penalty discontinuous Galerkin methods; Error analysis*

1 Introduction

In a time-domain fluid-solid interaction (FSI) problem, an incident acoustic wave is scattered by a bounded elastic obstacle immersed in a homogeneous, compressible and inviscid fluid. The problem of determining the scattered wave field plays prominent roles in many scientific and engineering areas, such as detecting and identifying submerged objects, geophysical exploration, medical imaging, seismology, and oceanography (such as **[30]**). Because of the difficulties of dealing with the time dependence, the FSI problem is usually studied based on a time-harmonic setting. Various efficient and accurate numerical methods have been developed to solve the time-harmonic FSI problem. Some of them provide proper guidance also for time-domain problems. Popular methods include the boundary integral equation (BIE) method, see **[33, 46]** and coupling methods such as the so-called coupled FEM-BEM method, see **[47]. An artificial boundary or an absorbing layer ([43]**) can be introduced to reduce the original unbounded problem to a bounded problem which could be solved using field equation solvers such as the finite element method.

The studies of this time-domain scattering problems now gain more and more attention since the time-domain model gives more information about the wave, more general material, and nonlinearity, see **[7, 39]**. There are relatively fewer mathematical analysis and numerical studies for the time-domain wave scattering problems. As to the numerical studies towards the time-domain wave scattering problems, the main challenge is how to handle the problem defined on the unbounded domain. Many approaches attempted to solve the time-domain problems numerically are developed, such as coupling of boundary element and finite element with different time quadratures, see **[17, 21]**, time-domain boundary elemment methods  **[26, 25]**, and to name a few. The authors of **[27]** gave the exact non-reflective boundary condition for the three dimensional time-domain wave scattering problem in 1995. For basic isotropy wave equation with constant coefficients, a planar PML method in one space direction designed for some particular domain is considered in **[29]**. In **[8]**, the mathematical analysis of a time-domain DtN operator and the convergence analysis of the PML method for acoustic wave scattering were given. In **[39]**, the time-domain exact nonreflective boundary conditions in both two and three dimensions were computed and analyzed. For the time-domain wave FSI problem, the rigorous mathematical study is still an open challenge, and related work include **[19, 31, 22, 6]**. Concerning its numerical solutions, the authors of  **[32]** applied time-domain boundary integral equations to time-domain FSI problem and analyzed the resulting nonlocal initial-boundary problem, motivated by the time-harmonic FSI problems. Coupling methods are also utilized in  **[17, 18, 20, 23, 24]** for solving the time-domain FSI problems.

Instead of numerical methods mentioned above, in this paper, we focus on using discontinuous Galerkin (DG) methods, which has natural advantages for dealing with the time-domain FSI problem, see **[16]**. The original DG methods were proposed for the numerical solution of hyperbolic neutron transport equations, as well as for discontinuities in some elliptic and parabolic problems. For second order wave equations, various DG methods have been proposed, which first reformulate the original wave problem to a first-order hyperbolic system, such as the local discontinuous Galerkin (LDG) methods in **[13, 42, 9]**, and we refer to **[10, 11, 34, 40]** for a review of some other DG methods for the first order wave equations. The first DG method for the original second-order formulation of acoustic wave equation was proposed in **[36]**, which is based on a nonsymmetric formulation. Here, we propose and analyze a symmetric interior penalty DG (SIPDG) method to solve the time-domain FSI problem. It is the same method used in **[28]** for the spatial discretization of the second-order scalar wave equation. Compared to the nonsymmetric formulation in **[36]**, the symmetric discretization of the second-order form wave equation offers extra benefits such as a positive definite stiffness matrix and hence is free of any (unnecessary) damping. One can refer to **[3]** for details of the DG methods for second order equations. Finally we note that another DG formulation for wave equations in second order for is the energy based method proposed in **[4]** and extended to the coupled acoustic-elastic problem in **[5]**.

The remainder of the paper is organized as follows. We first describe the original time-domain FSI problem in Section 2. Then in Section 3, the unbounded problem is reduced to a bounded initial-boundary value problem. In Section 4, we establish a priori error estimates for the IPDG solution of the reduced problem. Numerical experiments are presented in Section 5 to confirm the theoretical results, and finally, a conclusion is made in Section 6. For the sake of completeness, we provide the mathematical analysis towards well-posedness of the reduced problem introduced in Section 3 in the appendix.

2 Model

Here we study the same model as in **[31]**, where the authors gave mathematical analysis from the aspect of integral equation method. The statement of the model is as follows. Suppose $\Omega\subset\mathbb{R}^{2}$ is a bounded domain of a elastic body with boundary $\Gamma=\partial\Omega$, which is enclosed by the unbounded homogenous compressible inviscid fluid domain $\Omega^{+}=\mathbb{R}^{2}\backslash\overline{\Omega}$ (see Figure 1), and a finite time interval $J=(0,T)$. Given an incident wave $\varphi^{i}$, the scattered wave $\varphi$ is generated by the induced solid $\Omega$ in the fluid domain $\Omega^{+}$. In the fluid domain $\Omega^{+}$, the governing equations are the linearized Euler equation and the linearized equation of continuity

[TABLE]

where $\mathbf{v}$ is the velocity field, $p$, $\rho_{1}$ and $c$ are the pressure, the density and the speed of sound in the fluid respectively.

For an irrational flow, the velocity potential $\varphi$ s be the velocity potential satisfies:

[TABLE]

Then, the wave equation for $\varphi$ takes the form:

[TABLE]

The elastic wave in $\Omega$ is described by the displacement field $\mathbf{u}={(u_{1},u_{2})}^{T}$, which satisfies the dynamic linear elastic equation:

[TABLE]

Here $\rho_{2}$ is the constant density of the elastic body, which is assumed to be homogeneous and isotropic with the Lamé constants $\lambda$ and $\mu$ such that $3\lambda+2\mu>0$ and $\mu>0$.

The velocity potential $\varphi$ in equation (2.2) and the elastic displacement field $\mathbf{u}$ in equation (2.3) are coupled by the transmission conditions on $\Gamma\times J$, together with the homogeneous initial conditions, we can get an initial-boundary value problem:

[TABLE]

Here $\mathbf{n}$ is the unit outward vector on $\Gamma$ from $\Omega$ towards $\Omega^{+}$, $\varphi^{i}$ is the given incident field and equation (2.8) is the radiation condition for $\varphi$.

Before the discussion, we first introduce the definitions of some relevant Sobolev spaces and norms. Let $L^{2}(\Omega)$ be the function space consisting of all square integrable functions over $\Omega$ equipped with the norm

[TABLE]

For $s>0$, the standard Sobolev space is denoted by

[TABLE]

with the norm

[TABLE]

and $H^{s}(\Gamma)$ the trace functional space for $\Gamma=\partial\Omega$ under the $L^{2}(\Gamma)$ inner produce

[TABLE]

It is clear to note that $H^{-s}(\Omega)$ and $H^{-s}(\Gamma)$ are the dual space of $H^{s}(\Omega)$ and $H^{s}(\Gamma)$. We also denote $\mathbf{L}^{2}(\Omega)={(L^{2}(\Omega))}^{2}$ and $\mathbf{H}^{1}(\Omega)={(H^{1}(\Omega))}^{2}$. For any $\mathbf{u}={(u_{1},u_{2})}^{T}\in\mathbf{H}^{1}(\Omega)$, the Frobenius norm is defined as:

[TABLE]

And a simple calculation gives:

[TABLE]

where $C$ is a positive constant.

Furthermore, for $1\leq q\leq\infty$ we will make use of the Bochner space $L^{q}(J;H^{1}(\Omega))$, endowed with the norm:

[TABLE]

and the Bochner space $H^{1}(J;H^{s}(\Omega))$, endowed with norm:

[TABLE]

3 The reduced problem in bounded domain

To reduce this exterior problem to a problem in a bounded domain, we impose the first order approximate boundary condition on the artificial boundary $\Gamma_{R}$ (see Figure 2):

[TABLE]

where $\mathbf{n}$ is the outward unit vector from $\Omega_{R}$ on $\Gamma_{R}$. Then the reduced problem on the bounded domain reads:

[TABLE]

Then, we present the well-posedness and stability of the reduced problem in the follows.

3.1 Well-posedness and stability

First of all, we need to show that the reduced interaction problem in the bounded domain is well-posed and stable.

Theorem 3.1**.**

Suppose the incident wave $\varphi^{i}\in H^{1}(J;L^{2}(\Omega_{R}))$, $\varphi^{i}|_{t=0}=0$ and $spt(\varphi^{i})\subset\Omega_{R}\times J$. Then (3.1)-(3.6) has a unique solution satisfying:

[TABLE]

[TABLE]

and we have the following stability estimate:

[TABLE]

and

[TABLE]

Although different boundary conditions are considered to develop the reduced problem, the proof of above theorem 3.1 is analogous to the proof of Theorem 3.2 in **[6]**. We omit to present the proof of Theorem 3.1 in this section, and report the detailed proof of Theorem 3.1 in the appendix for the sake of completeness.

3.2 Variational formulation

The standard variational formulation of problem (3.1)-(3.6) is as follows: For given incident wave $\varphi^{i}$, find

[TABLE]

with

[TABLE]

such that, $\forall~{}\mathbf{v}\in\mathbf{H^{1}}(\Omega)\forall~{}\phi\in H^{1}(\Omega_{R})$:

[TABLE]

with the initial conditions

[TABLE]

Here $\langle\cdot,\cdot\rangle$ denotes the duality pairing between spaces $H^{-1}$ and $H^{0}$ with the associated domain, and $a(\mathbf{u},\varphi;\mathbf{v},\phi)$ and $L(\mathbf{v},\phi)$ are defined as:

[TABLE]

with $(\cdot,\cdot)$ the standard inner product in $L^{2}$-space.

4 The IPDG method

4.1 Spaces, jumps and averages

Assume that $\Omega$ with Lipschitz boundary $\Gamma$ is regularly divided into disjoint elements ${E}$ by mesh ${\mathbf{\mathcal{E}_{h}}}$ such that $\overline{\Omega}=\bigcup_{E\in\mathbf{\mathcal{E}_{h}}}\overline{E}$, where $E$ is a triangle or quadrilateral in 2D, or a tetrahedron or hexahedron in 3D. Similarly, The regular meshes $\mathcal{K}_{h}$ partitions the fluid domain $\Omega_{R}$ into disjoint elements $K$ such that $\overline{\Omega}_{R}=\bigcup_{K\in\mathcal{K}_{h}}\overline{K}$. The diameter of element $E/K$ is denoted by $h_{E/K}$, and $h$ is the mesh size given by $h=\max\limits_{E\in\mathbf{\mathcal{E}}_{h},K\in\mathcal{K}_{h}}(h_{E},h_{K})$.

The boundary of the elastic body $\Omega$, i.e., $\Gamma$ is approximated by the boundary edges of the subdivision: The set of boundary edges of $\mathbf{\mathcal{E}}_{h}$ is $\Gamma_{h}:=\partial E\bigcap\Gamma$, since $\Omega\bigcap\Omega_{R}=\Gamma$, the edges in $\Gamma_{h}$ are also the boundary edges of $\mathcal{K}_{h}$, so the sets of boundary edges of $\mathcal{K}_{h}$ are $\Gamma_{h}\bigcup\Gamma_{h}^{R}$, $\Gamma_{h}^{R}:=\partial K\bigcap\Gamma_{R}$. We denote by $\Gamma_{I}^{1}$ and $\Gamma_{I}^{2}$ the set of interior edges of the subdivision $\mathbf{\mathcal{E}}_{h}$ and $\mathcal{K}_{h}$ respectively.

For any real number $m$, the broken Sobolev spaces $H^{m}(\mathcal{E}_{h})$ and $H^{m}(\mathcal{K}_{h})$ are defined as :

[TABLE]

and

[TABLE]

with the broken Sobolev norms:

[TABLE]

In particular, we also define the broken gradient seminorm:

[TABLE]

If $\varphi\in H^{m}(\mathcal{K}_{h})$, the trace of $\varphi$ along any side of each element $K\in\mathcal{K}_{h}$ is well defined. Let $e$ be the edge between the elements $K_{1}$ and $K_{2}$, then the jump and average of a scalar function $\varphi$ on $e$ are given by:

[TABLE]

respectively, where $\mathbf{n}_{K_{1}}$ is the outward unit normal from $K_{1}$ to $K_{2},$ and likewise for $\mathbf{n}_{K_{2}}$. On the boundary edge $e\in\Gamma_{h}$ or $e\in\Gamma_{h}^{R}$, we extend the definition: $\{\varphi\}=[\varphi]=\varphi|_{K^{e}}$, where $K_{e}$ is the element that $K_{e}\bigcap\Gamma_{h}=e$ or $K_{e}\bigcap\Gamma_{h}^{R}=e$.

Similarly, for a vector valued function $\mathbf{v}\in\mathbf{H}^{s}(\mathbf{\mathcal{E}}_{h})$, the jump and average of a vector function $\mathbf{v}$ on $e$ are given by:

[TABLE]

respectively, where $\mathbf{n}_{E_{1}}$ is outward unit normal from $E_{1}$ to $E_{2},$ and likewise for $\mathbf{n}_{E_{2}}$. Similarly, on the boundary edge $e$: $\{\mathbf{v}\}=[\mathbf{v}]=\mathbf{v}|_{E^{e}}$, where $E_{e}$ is the element that $E_{e}\bigcap\Gamma_{h}=e$.

4.2 Spatial discretization

For given partitions $\mathbf{\mathcal{E}}_{h}$ of $\Omega$, $\mathcal{K}_{h}$ of $\Omega_{R}$ and an approximation order $k\geq 1$, we can approximate the solution $(\mathbf{u},\varphi)$ of (3.1)-(3.6) in the finite element subspaces

[TABLE]

and

[TABLE]

of the broken Sobolev space $H^{m}(\mathcal{E}_{h})$ and $H^{m}(\mathcal{F}_{h})$ for $m>3/2$. Here $\mathbb{P}_{k}$ denotes the space of polynomials of degree at most $k$.

Then we consider the following semidiscrete DG approximation of (3.1)-(3.6): find $\mathbf{u}^{h}\in\mathbf{\mathcal{D}}_{k}(\mathbf{\mathcal{E}}_{h})$ and $\varphi^{h}\in\mathcal{D}_{k}(\mathcal{K}_{h})$ such that

[TABLE]

holds for any $\mathbf{v}\in\mathbf{\mathcal{D}}_{k}(\mathbf{\mathcal{E}}_{h})$ and $\phi\in\mathcal{D}_{k}(\mathcal{K}_{h})$. Here, the discrete bilinear form $a_{h}(\mathbf{u}^{h},\varphi^{h};\mathbf{v},\phi)$ and linear form $L_{h}(\mathbf{v},\phi)$ are given by the IPDG discretization as:

[TABLE]

The interior penalty stabilization function $\frac{\alpha}{|e|^{\beta}}$ penalizes the jumps over the edges of $\mathbf{\mathcal{E}}_{h}$ and $\mathcal{K}_{h}$, where $\alpha$ is a positive parameter independent of the local mesh sizes. Here $|e|$ simply means the length of $e$ and we have

[TABLE]

with $\beta>0$, and when $\beta>1$ the method is superpenalized.

We define the space $\mathbf{V}_{1}(h)=\mathbf{H}^{1}(\Omega)+\mathbf{\mathcal{D}}_{k}(\mathbf{\mathcal{E}}_{h})$ and $V_{2}(h)=H^{1}(\Omega_{R})+\mathcal{D}_{k}(\mathcal{K}_{h})$. On $\mathbf{V}_{1}(h)$ and $V_{2}(h)$, we define the DG energy norm as

[TABLE]

and

[TABLE]

The consistency of the scheme is straightforward since the jumps at element boundaries vanishes when $\mathbf{u},\mathbf{v}\in\mathbf{H}^{m}(\mathbf{\mathcal{E}}_{h})$ and $\varphi,\phi\in H^{m}(\mathcal{K}_{h})$ for $m>2/3$. Next we will discuss the property of the bilinear form $a_{h}(\cdot,\cdot;\cdot,\cdot)$. For more details about the IPDG method, one can refer to **[28]**

Lemma 4.1** (Coercivity).**

There exists a positive constant $C_{\text{coer}}$ independent of $h$ such that for all $\mathbf{v}\in\mathbf{V}_{1}(h)$ and $\phi\in V_{2}(h)$

[TABLE]

Proof:* By Cauchy-Schwarz inequality,*

[TABLE]

Consider the average of the fluxes for an interior edge $e$ shared by $E_{1}^{e}$ and $E_{2}^{e}$ and apply the trace theorem, we have:

[TABLE]

Assume $\beta\geq 1$ and $h\leq 1$, we have:

[TABLE]

We denote $n_{0}$ the maximum number of neighbors an element can have.

[TABLE]

Using Young’s inequality, we have for $\delta>0$

[TABLE]

Following a similar steps, we get:

[TABLE]

Thus we obtain a lower bound for $a_{h}(\mathbf{v},\phi;\mathbf{v},\phi)$:

[TABLE]

For any $0<t\leq T$, using the Young’s inequality, similar as (38) we have

[TABLE]

Then, we have

[TABLE]

Since $\sigma(\mathbf{v})=(\lambda\,\text{div}(\mathbf{v}))\mathrm{I}+2\mu\varepsilon(\mathbf{v})$ and triangle inequality,

[TABLE]

Therefore, we have

[TABLE]

From the definition of the strain tensor, $\|\varepsilon(\mathbf{v})\|_{F(E)}\leq\|\nabla\cdot\mathbf{v}\|_{0,E}$, we have

[TABLE]

Assume $\alpha$ is large enough (e.g. $\alpha\geq\frac{C_{+}^{2}n_{0}+C}{\delta}$), take $C_{\text{coer}}=\lambda+2\mu-2\delta$, we obtain the coercivity result for $a_{h}(\mathbf{v},\phi;\mathbf{v},\phi)$. ∎

Lemma 4.2** (Continuity).**

There exists a positive constant $C_{\text{cont}}$ independent of $h$ such that

[TABLE]

Proof:* Using Cauchy-Schwartz inequality, we obtain:*

[TABLE]

Similarly, we have

[TABLE]

According to Cauchy-Schwartz inequality and the trace theorem:

[TABLE]

[TABLE]

and

[TABLE]

where all constants $C>0$ are independent of $\mathbf{u},\varphi,\mathbf{v},$ and $\phi$. Set $C_{\text{cont}}=\max\left\{C,\frac{\lambda+8\mu}{\rho_{1}}\right\}$, then combining (4.6)-(4.8) together and applying the arithmetic geometric average inequality, we obtain (4.5). ∎

4.3 A priori error estimate

In this section, we will give a priori error estimates in the energy norm for the IPDG method. First of all, we will derive an error equation. For $(\mathbf{u},\varphi)\in\mathbf{H}^{1+m}\times H^{1+m}$ with $m>1/2$ and $(\mathbf{v},\phi)\in\mathbf{V}_{1}(h)\times V_{2}(h)$, we define

[TABLE]

Here $\Pi_{h}$ denotes the $L^{2}$-projection onto the associated finite element space ($\mathbf{V}_{1}(h)$* or $V_{2}(h)$). The assumption $(\mathbf{u},\varphi)\in\mathbf{H}^{1+m}\times H^{1+m}$ ensures that $r_{h}(\mathbf{u},\varphi;\mathbf{v},\phi)$ is well defined. From the definition (4.11) it is obvious that $r_{h}(\mathbf{u},\varphi;\mathbf{v},\phi)=0$ when $(\mathbf{v},\phi)\in\mathbf{H}^{1+m}\times H^{1+m}$.*

Lemma 4.3**.**

Assume the analytical solution $(\mathbf{u},\varphi)$ of (3.1)-(3.6) satisfy

[TABLE]

and

[TABLE]

Let $(\mathbf{u}^{h}$, $\varphi^{h})$ be the semidiscrete DG approximation solution. Then, the error $(\mathbf{e}^{u},e^{\varphi})=(\mathbf{u}-\mathbf{u}^{h},\varphi-\varphi^{h})$ satisfies:

[TABLE]

for $(\mathbf{v},\phi)\in\mathbf{\mathcal{D}}_{k}(\mathbf{\mathcal{E}}_{h})\times\mathcal{D}_{k}(\mathcal{K}_{h}),$ and $r_{h}$ is given in (4.11).

Proof:* Let $(\mathbf{v},\phi)\in\mathbf{\mathcal{D}}_{k}(\mathbf{\mathcal{E}}_{h})\times\mathcal{D}_{k}(\mathcal{K}_{h})$. Since $\mathbf{u}_{tt}\in L^{1}\left(J;\mathbf{L}^{2}(\Omega)\right)$, $\varphi_{tt}\in L^{1}\left(J;L^{2}(\Omega_{R})\right)$, we have $\langle\mathbf{u}_{tt},\mathbf{v}\rangle=(\mathbf{u}_{tt},\mathbf{v}),~{}\langle\varphi_{tt},\phi\rangle=(\varphi_{tt},\phi)$ almost everywhere in $[0,T]$. Hence, using the discrete formulation in (4.1), we obtain that*

[TABLE]

By the definition of $a_{h}$, the fact that $[\mathbf{u}]=0$, $[\varphi]=0$ on all edges and the properties of the $L^{2}$-projection $\Pi_{h}$, we obtain

[TABLE]

Since $\mathbf{u}_{tt}\in L^{1}\left(J;\mathbf{L}^{2}(\Omega)\right)$, $\varphi_{tt}\in L^{1}\left(J;\mathbf{L}^{2}(\Omega_{R})\right)$ and $\varphi^{i}\in L^{1}\left(J;L^{2}(\Omega_{R})\right)$, we have that $\text{div}(\sigma(\mathbf{u}))\in L^{2}(\Omega)$ and $\Delta\varphi\in L^{2}(\Omega_{R})$ almost everywhere in $[0,T]$, which means that $\text{div}(\sigma(\mathbf{u}))$ and $\Delta\varphi$ have continuous normal components across all interior faces. Therefore, using integration by parts elementwisely and combining with the trace operators yield

[TABLE]

From the definition of $r_{h}(\mathbf{u},\varphi;\mathbf{v},\phi)$ in (4.11), we conclude that

[TABLE]

and obtain

[TABLE]

where we have used the differential equations (3.1) and (3.2). ∎

Next, we recall some approximation properties, see **[12]** for more detail.

Lemma 4.4**.**

*Let $E\in\mathbf{\mathcal{E}}_{h}$, $K\in\mathcal{K}_{h}$. Then the following hold:

(i) For $(\mathbf{v},\phi)\in\mathbf{H}^{m}(E)\times H^{m}(K)$, $m\geq 0$, we have*

[TABLE]

(ii) For $(\mathbf{v},\phi)\in\mathbf{H}^{1+m}(E)\times H^{1+m}(K)$, $m>\frac{1}{2}$, we have

[TABLE]

where $C$ is independent of local mesh size $h$.

As a consequence of the approximation properties in Lemma 4.4, we obtain the following results.

Lemma 4.5**.**

*Let $(\mathbf{v},\phi)\in H^{1+m}(E)\times H^{1+m}(K)$, $m>\frac{1}{2}$. Then the following hold:

(i) We have*

[TABLE]

*with a constant $C_{A}$ that is independent of the mesh size.

(ii) For $(\mathbf{v},\phi)\in\mathbf{V}_{1}(h)\times V_{2}(h)$, the form $r_{h}(\mathbf{u},\varphi;\mathbf{v},\phi)$ in (4.11) can be bounded by*

[TABLE]

with a constant $C_{R}$ that is independent of $h$ and depends only $\alpha$ and the constants in Lemma 4.4.

Lemma 4.6**.**

Assume the analytical solution $(\mathbf{u},\varphi)$ of (3.1)-(3.6) satisfy

[TABLE]

and

[TABLE]

for a regularity exponent $m>1/2$. Let

[TABLE]

Then we have

[TABLE]

where $C_{R}$ is the constant from the bound (ii) in Lemma 4.5.

Proof:* From the definition of $r_{h}$ in (4.11) and integration by parts, we obtain*

[TABLE]

Lemma 4.5 implies

[TABLE]

and

[TABLE]

which finishes the proof of the lemma.∎

Theorem 4.7**.**

Assume that the solution of (3.1)-(3.6) satisfy

[TABLE]

and

[TABLE]

for a regularity exponent $m>1/2$, and let $\mathbf{u}^{h}$, $\varphi^{h}$ be the semidiscrete discontinuous Galerkin approximation obtained by (4.1). Then the error $\mathbf{e}^{u}=\mathbf{u}-\mathbf{u}^{h}$ and $e^{\varphi}=\varphi-\varphi^{h}$ satisfy the estimate:

[TABLE]

and

[TABLE]

with the constant $C$ that is independent of $h$ and $T$.

Proof:* Because of Theorem 3.1, we have*

[TABLE]

Using the symmetry of $a_{h}$ and the error equation in Lemma 4.3, we obtain

[TABLE]

We fix $\tau\in[0,T]$ and integrate (4.13) over the time interval $(0,\tau)$. This yields

[TABLE]

Integration by parts of the fifth and sixth terms on the right-hand side yields

[TABLE]

and

[TABLE]

From the stability properties of $a_{h}$ in Lemma 4.1 and Lemma 4.2 and the Hölder’s inequalities, we conclude that

[TABLE]

Since this inequality holds for every $\tau\in J$, it also holds for the maximum over $J$, that is

[TABLE]

with

[TABLE]

Using the inequality that $|ab|\leq\frac{1}{2\epsilon}a^{2}+\frac{\epsilon}{2}b^{2}$, and the approximation results in Lemma 4.4, we conclude that

[TABLE]

and

[TABLE]

with the constant $C$ depends only on the constants in Lemma 4.4. Similarly,

[TABLE]

and

[TABLE]

with the constant $C$ depends only on $C_{\text{cont}}$, $C_{\text{coer}}$ and the constants $C_{A}$ in Lemma 4.5.

To bound the term $T_{5}$, we use the Lemma 4.6 to obtain

[TABLE]

with

[TABLE]

The triangle inequality, the geometric-arithmetic mean, and the approximation properties of $\Pi_{h}$ in Lemma 4.5 then yield

[TABLE]

with the constant $C$ depends only on $C_{R}$, $C_{\text{coer}}$ and the constants $C_{A}$. Combining the above estimates for $T_{1},T_{2},T_{3},T_{4}$ and $T_{5}$ then shows that

[TABLE]

with a constant that is independent of $T$ and the mesh size. This concludes the proof of Theorem 4.7.∎

Hence, with a standard projection approximation for the initial condition, Theorem 4.7 yields an optimal convergence estimate in the (DG) energy norm

[TABLE]

and

[TABLE]

with a constant $C$ that is independent of $h$.

5 Numerical results

We will present the numerical results to verify the theoretical error analysis, and the IPDG penalty parameters are chosen as follows: $\alpha=100,\beta=1$. The discretization of the FSI problem in space by the IPDG method leads to the linear second-order system of ordinary differential equation as below

[TABLE]

with the initial conditions

[TABLE]

Here $M$, $N$ and $A$ are the known coefficient matrices. We use the second order Newmark time stepping scheme **[35]** to discretize (5.1) in time domain, and $l$ denotes the time step size with $t_{n}=nl$. Then the Newmark method consists of finding $\{U^{h}_{n}\}_{n}$ to $U^{h}(t_{n})$ such that

[TABLE]

for $n=1,\cdots,N-1$. Here, $\mathbf{f}_{n}:=\mathbf{f}(t_{n})$, while $\delta\geq 0$ and $\gamma\geq 1/2$ are free parameters that still can be chosen. In our examples, we use the explicit second-order Newmark scheme, setting $\gamma=1/2$ and $\delta=0$. In our accuracy tests, we consider the following two models with the parameters $c=1$, $\rho_{1}=1$, $\rho_{2}=1$, $\lambda=1$, $\mu=1$. The computational domain $\Omega$ and $\Omega_{R}$ are approximated by uniform triangle elements and piecewise linear basis functions ($k=1$) are employed.

Example 1

Let $\Omega$ be a disk with radius $R_{0}=1$, a plane incident wave $\varphi^{i}=\cos(\mathbf{x}\cdot\mathbf{d})\cos(t)$ is given, where the direction is $\mathbf{d}=(1,0)$. Choosing the artificial boundary $\Gamma_{R}$ to be a circle with the radius $R=2$, sharing the same center of $\Gamma$,.

*Let the final time $T=1$ and use the numerical solution $U^{r}$ with mesh size $h=0.0134$ as a reference solution. We set time step $l=h/20$ to ensure the stability of time discretization. The numerical errors $U^{r}-U^{h}$ and convergence order are presented in Table 1 for various mesh sizes. The numerical results confirm the expected rates of *$k-$*th order for the energy norm and *$(k+1)-$th order for the $L^{2}$ norm.

Example 2

We consider the FSI problem on a $L$-shaped domain $\Omega=(-1,1)^{2}\setminus[0,1)^{2}$. The incident wave is a point source at $(2,0)$:

[TABLE]

The artificial boundary $\Gamma_{R}$ is a circle centering at $(0,0)$ with radius $R=3$. Here, although solution $(\mathbf{u},\varphi)$ is continuous respect to time, it has a singular point $(0,0)$, that is $\mathbf{u}\in C^{\infty}(J;\mathbf{H}^{2/3}(\Omega))$, $\varphi\in C^{\infty}(J;H^{2/3}(\Omega_{R}))$.

Similarly, we let the numerical solution $U^{r}$ on the mesh with mesh size $h=0.0227$ be the reference solution and $l=h/20$, and compute the errors $U^{r}-U^{h}$ in the energy norm and $L^{2}$ norm at $T=1$. In this case, the parameter $m$ in Theorem 4.7 is $2/3$, so the theoretical convergence rate in energy norm should be $2/3$. The results in Table 2 validate this conclusion.

6 Conclusion

In this work, we have proposed an interior penalty discontinuous Galerkin method for solving the acoustic-elastic wave interaction problem. A low-order approximate absorbing boundary condition has been used to deal with the acoustic waves in the unbounded domain. Essential analysis, including a priori error analysis, have been performed for the discontinuous Galerkin solution. We will consider the analysis of stable and conserved fully discrete schemes for the interaction problem, and the interaction problem of wave propagation in orthotropic porous elastic media in the future. Numerical schemes using more accurate coupling method of discontinuous Galerkin and boundary integral equation methods will also be envisioned in our future work.

Appendix. Proof of Theorem 3.1

This appendix is devoted to the proof of Theorem 3.1. As we mentioned in Section 3, such a proof is technically complicated. Here, we prove this theorem by using a priori estimate of an elliptic problem and the abstract inversion theorem of the Laplace transform. The following lemma (Treve **[38]**, Theorem 43.1) is the analog of the Paley-Wiener-Schiwarz theorem for the Fourier transform of the distributions with compact support in the case of the Laplace transform.

Lemma 1**.**

*Let $\mathbf{h}(s)$ denote a holomorphic function in the half-plane $Re(s)>\alpha_{0}$, valued in the Banach space $E$. The following conditions are equivalent:

(i) there is a distribution $T\in D^{\prime}_{+}(E)$ whose Laplace transform is equal to $\mathbf{h}(s)$;

(ii) there is a real $\alpha_{1}$, $\alpha_{0}\leq\alpha_{1}<\infty$, a constant $C>0$, and an integer $m\geq 0$ such that, for all complex numbers $s$ with $Re(s)>\alpha_{1}$, it holds that*

[TABLE]

where $\|\cdot\|_{E}$ is the norm of Banach space $E$, and $D^{\prime}_{+}$ is the space of distributions on the real line which vanish identically in the open negative half-line.

Take the Laplace transform of (3.1)-(3.5) and denote $\tilde{\mathbf{u}}=\mathcal{L}(\mathbf{u})$, $\tilde{\varphi}=\mathcal{L}(\varphi)$, where $\mathcal{L}$ is the Laplace transform operator, we obtain:

[TABLE]

Lemma 2**.**

Let $s=s_{1}+is_{2}$, $s_{1}>0$, $s_{2}\in\mathbb{R}$. For any $g\in L^{2}(\Omega_{R})$, let $(\tilde{\mathbf{u}},\tilde{\varphi})$ be the solution of (1)-(5). Then there exists a constant $C$ independent of $s$ such that:

[TABLE]

Proof.

Multiplying (1) and (2) by test functions $\mathbf{v}\in\mathbf{H}^{1}(\Omega)$ and $\phi\in H^{1}(\Omega_{R})$ respectively, then using integration by parts, the symmetry of the stress tensor, and the transmission and boundary conditions, we get

[TABLE]

Multiply (8) with $\frac{1}{\rho_{1}s}$ and add it to (9), we get:

[TABLE]

where

[TABLE]

Using Cauchy-Schwarz inequality and the trace theorem, we have:

[TABLE]

where $C_{1}$ is a constant that independent of $\tilde{\mathbf{u}},\tilde{\varphi},\mathbf{v},\phi$, hence the sesquilinear form $\tilde{a}(\tilde{\mathbf{u}},\tilde{\varphi};\mathbf{v},\phi)$ is bounded.

Letting $(\mathbf{v},\phi)=(\tilde{\mathbf{u}},\tilde{\varphi})$ in (11) yields:

[TABLE]

Taking real part of (13), we have

[TABLE]

where $C_{2}$ is a constant that independent of $s$. According to Lax-Milgram lemma, there exists a unique solution $(\tilde{\mathbf{u}},\tilde{\varphi})\in\mathbf{H}^{1}(\Omega)\times H^{1}(\Omega_{R}),$ satisfies the following problem

[TABLE]

Based on equation (15) and Cauchy-Schwarz inequality, there exists a constant $C_{3}$ such that

[TABLE]

Combining (14) and (16):

[TABLE]

Letting $C=\left(\frac{C_{3}}{C_{2}}\right)^{2}$, we obtain (6) and (7). ∎

Now we are ready to prove Theorem 3.1.

Proof of Theorem 3.1:

Using Lemma 2, $\tilde{\mathbf{u}},\tilde{\varphi}$ are the holomorphic functions of $s$ on the half plane $s_{1}>\gamma>0$. Hence the inverse Laplace transform of $\tilde{\mathbf{u}}$ and $\tilde{\varphi}$ exist i.e., the problem (3.1)-(3.6) has a unique solution $(\mathbf{u},\varphi)=(\mathcal{L}^{-1}(\tilde{\mathbf{u}}),\mathcal{L}^{-1}(\tilde{\varphi}))$.

For $\varphi$,we have:

[TABLE]

Similarly, for $\mathbf{u}$, we have:

[TABLE]

From estimate (7) and the trace theorem, we have:

[TABLE]

Using Parseval identity (see **[14]**) and estimate (20):

[TABLE]

which shows that

[TABLE]

Similarly, we have

[TABLE]

It follows from (23) that

[TABLE]

Next we prove the stability, which also helps establishing numerical stability of the IPDG scheme in the following section. For any $0<t<T$, consider the energy function:

[TABLE]

where

[TABLE]

and

[TABLE]

Obviously,

[TABLE]

By using the integration by parts and (3.1),(3.3),(3.6) that

[TABLE]

Similarly, follows from the integration by parts and (3.1),(3.3),(3.6), we have

[TABLE]

Since $E(0)=0$, by combining equations (25)-(27) and using the trace theorem, we have:

[TABLE]

Therefore, using Young’s inequality, we obtain:

[TABLE]

Since the norm on the right-hand side in (29) contains $\mathbf{u}_{t}$ and $\varphi_{t}$, which cannot be bounded by the left-hand side, here a new system is needed. Take the first partial derivative of (3.1)-(3.6) with respect to $t$, we have:

[TABLE]

We consider the energy function:

[TABLE]

where

[TABLE]

and

[TABLE]

It is clear that $F(0)=0$. Similarly, follows the same step to prove inequality (28), we obtain that

[TABLE]

Combining (29) and (36), using Young’s inequality, we have:

[TABLE]

Hence we choose $\epsilon>0$ small enough such that $2C\epsilon<1/2$. It follows from (37) that

[TABLE]

For any $0<t\leq T$, using the Young’s inequality leads to

[TABLE]

Here we choose $\epsilon$ small enough so that $\epsilon T<1$ (e.g. $\epsilon=\frac{1}{2T}$). Hence, we have:

[TABLE]

Similarly, we obtain

[TABLE]

Combining (37)-(39), we have

[TABLE]

after applying the Cauchy-Schwarz inequality we have the estimate (3.7) and (3.9). ∎