Acoustic vibration problem for dissipative fluids

Abstract

In this paper we analyze a finite element method for solving a quadratic eigenvalue problem derived from the acoustic vibration problem for a heterogeneous dissipative fluid. The problem is shown to be equivalent to the spectral problem for a non-compact operator and a thorough spectral characterization is given. The numerical discretization of the problem is based on Raviart-Thomas finite elements. The method is proved to be free of spurious modes and to converge with optimal order. Finally, we report numerical tests which allow us to assess the performance of the method.

1. Introduction

This paper deals with the numerical approximation of an acoustic dissipative fluid system. This kind of problem has attracted much interest, since it is frequently encountered in engineering applications (see, for instance, 52032). One typical example is to achieve optimal designs that reduce noise and vibrations in fluid-structure systems like cars, aircraft, or dams.

Although dissipation is usually neglected in standard acoustics, modeling this phenomenon is essential to study the effect of noise reduction techniques. Indeed, in most real situations, damping mechanisms that transform mechanical energy into heat do exist. Sometimes these mechanisms are based on surface damping arising from viscoelastic materials placed on the boundary of the propagation domain. In these cases, the dissipative effects are typically included in the model by means of a surface impedance in the boundary conditions (see, for instance, 467). The present paper addresses damping when it arises in the propagation media itself due to friction and heat conduction. A general approach to this topic can be found in the books by Landau and Lifshitz 22, Morse 31, and Pierce 37, all of which include extensive bibliographic references on the subject.

This paper focuses on computing the (complex) vibration frequencies and modes of an acoustic dissipative fluid system within a rigid cavity. One motivation for considering this problem is that it constitutes a stepping stone towards the more challenging goal of devising numerical approximations for coupled systems involving fluid-structure interaction between viscous fluids and solid structures. The natural model for the fluid system should be based on the Stokes equations for compressible fluids. However, since in real applications the viscosity is typically very small, the resulting problem turns out to be a singular perturbation of the problem for an inviscid fluid. This fact leads to a kind of dilemma, since appropriate finite elements for the Stokes equations introduce spurious modes in the limit case of a vanishing viscosity, whereas the finite elements that avoid such spectral pollution fail when applied to the Stokes equation.

To circumvent this drawback, we resort to an alternative model based on a curl-free displacement formulation (see 8 for the derivation of a similar model in the time domain from basic mechanical laws). Let us remark that in principle the fluid displacement does not need to be curl-free. However, since the viscosity term due to vorticity is typically very small, except perhaps near the walls of the enclosure, it may be neglected in the interior of the enclosure and eventually modeled as a wall impedance on its boundary (see 32 for a similar model).

This curl-free displacement formulation for a viscous fluid leads to a quadratic eigenvalue problem (QEP), as happens in 4. However, the resulting problem involves additional challenges related to the fact that the essential spectrum does not reduce to a single point as in 4. In fact, in this case, we can only prove that the essential spectrum is well-separated from the physically relevant eigenvalues when the viscosity is sufficiently small (as happens in practice). On the other hand, the associated solution operator is not regularizing. Because of this, we need to split it into two terms for the numerical analysis. One of them is dealt with from the techniques in 4, but the spectral approximation analysis for the other is new.

As is shown below (cf. Remark 2.1), the vibration frequencies and modes of a viscous homogeneous irrotational fluid within a rigid cavity can be obtained without actually solving a QEP. In fact, these frequencies can be algebraically computed from those of the analogous inviscid fluid, whose approximation is nowadays a well-known subject (see, for instance, 5). However, this is not the case for a heterogeneous fluid and this is the reason why we choose this as our model problem. In particular, we consider the QEP arising from the acoustic vibration problem for a dissipative fluid system that consists of two homogeneous viscous immiscible fluids contained in a rigid cavity.

QEP has many applications in the study of vibration for solid systems, acoustic fluids, electrical circuits, etc., where the damping effects are involved. A state of the art work on the QEP up to the beginning of this century can be found in 39. However, there are not many works with a rigorous mathematical framework in the context of the numerical approximation of the eigenvalue problem of a partial differential operator involving damping. The first article proving this type of results is 4, where the authors have considered a displacement formulation for a fluid in a rigid cavity with absorbing walls. The theory of non-compact operators of 16 is used to obtain error estimates with minor modifications due to non-conformity.

On the other hand, alternative formulations for the absorbing wall problem have been studied in the engineering literature. For instance, a formulation for the QEP in terms of the fluid pressure has been proposed in 23. This type of formulation leads to a rational eigenvalue problem, for which different algorithms to compute the vibration frequencies have been introduced. Alternatively, two formulations, one based on the fluid pressure and the other on the fluid displacement, have been considered in 14, where an improved Arnoldi algorithm has been proposed to solve the discrete problem. On the other hand, an application of the damping effects in electromechanical-thermoelastic systems is presented in 1. Moreover, an a posteriori estimator for a QEP contextualized in the photonic crystal applications is proposed in 17. Nevertheless, all the previously mentioned studies present different numerical technologies to solve the QEP, but without a rigorous mathematical analysis. Such a rigorous analysis is present instead in the recent paper 27, where an efficient multiscale technique based on localized orthogonal decompositions to solve discrete generic damped vibration problems has been proposed (see also 26).

In the present paper, we consider a displacement-based variational formulation of the transmission eigenproblem resulting from our physical model. This approach leads to a QEP, which is transformed into an equivalent double-size linear eigenvalue problem that fits within the functional framework for non-self-adjoint and non-compact bounded operators. At the continuous level, we follow 21 to obtain an appropriate spectral characterization. Next, we propose an $\mathrm{H}(\text{div})$-conforming mixed finite element approximation of the problem and adapt the abstract spectral approximation theory for non-compact operators developed in 1516 to prove that the spectrum is correctly approximated and to obtain error estimates.

The discrete analysis relies partly on the techniques used in the Raviart-Thomas mixed approximation of the grad-div eigenvalue problem. This spectral problem emerged in the study of coupled fluid-solid systems in 3 (see also 28 for a similar setting in elasticity). The grad-div spectral problem is posed in $\mathrm{H}(\text{div})$ and it is closely related to the Maxwell eigenvalue problem, which is formulated in terms of the curl-curl operator in $\mathrm{H}(\mathbf{curl})$ (see 101330). Although the two spectral problems have been initially studied in isolation from each other, a common framework now becomes clear thanks to the language of differential forms and the approach based on the finite element exterior calculus provided by 2 (see also 11, Part 4).

With $X$ representing either $\mathrm{H}(\text{div})$ or $\mathrm{H}(\mathbf{curl})$, we let $K$ be the kernel of the div or the curl operator in the corresponding space. If we denote by $X_{h}$ a conforming mixed finite element approximation of $X$ and consider $K_{h}:=K\cap X_{h}$, then the main tool in the proof of the correct spectral approximation relies on the construction of a projector $P$ onto $K^{\bot }$ such that

$$\begin{equation} \lim _{{h\to 0}}\|I - P\|_{\mathcal{L}(K_{h}^{\bot }, X)} = 0. {\tag{1.1}} \end{equation}$$

Here $K^\bot$ stands for $\mathrm{L}^{2}$-orthogonal complement of $K$ in $X$ and $K_{h}^\bot$ represents the $\mathrm{L}^{2}$-orthogonal complement of $K_{h}$ in $X_{h}$.

To our knowledge, this tool was first introduced in the $\mathrm{H}(\text{div})$-setting in 3 and more recently in 28, Lemma 4.2. The same tool appeared independently in the context of the Maxwell eigenvalue problem; see for example 18, Lemma 4.5 where the projector $P$ is named the Hodge mapping. A related abstract approach is considered in 12 where the so-called GAP property is introduced. It is shown in 12 that the GAP property is equivalent to 1.1 when $K_{h}\subset K$ (as in the cases considered here). It is also shown that the GAP property implies the discrete compactness property which was used in 13 as another tool to study edge element approximation of the Maxwell eigenproblem. See also the discussion given in 11, Section 19 regarding this topic in the language of differential forms.

The negative impact that material parameters have on the regularity of the solution of the boundary value problem complicates the analysis in this common framework (see 18, Remark 13). Here, we follow the lines of the methodology presented originally in 3 and use the information about interface singularities of solutions of the Neumann boundary value problem for $\mathop{\mathrm{div}}\nolimits (\kappa \mathop{\mathbf{grad}}\nolimits$), with $\kappa$ piecewise constant, to obtain (1.1) for the eigenproblem considered in this paper.

The paper is organized as follows: in Section 2, we introduce the spectral problem and the corresponding variational formulation, which leads to a quadratic eigenvalue problem. We introduce an auxiliary unknown to transform the quadratic eigenvalue problem into a linear one. Moreover, we introduce the corresponding solution operator for the spectral problem. In Section 3, we provide a thorough spectral characterization of the solution operator, based on the theory developed in 21. We also consider the limit problem (i.e., the case when the viscosity vanishes) and the relation between the solutions of the dissipative and non-dissipative problems. In Section 4, we introduce a finite element discretization using Raviart-Thomas elements for both fluids and imposing the continuity of the corresponding normal components on the interface. We analyze the discrete spectral problem analogously as in the continuous case and introduce the corresponding discrete solution operator. We use the abstract theory from 15 to prove the convergence. We also prove error estimates for our problem by adapting the arguments from 4. Finally, in Section 5, we report some numerical tests which allow us to asses the performance of the proposed method.

Throughout the paper, $\Omega$ is a generic Lipschitz bounded domain of $\mathbb{R}^d$ ($d=2,3$), with outer unit normal vector $\boldsymbol{n}$. We denote by $\mathcal{D}(\Omega )$ the space of infinitely smooth functions compactly supported in $\Omega$. For $r\geq 0$, $\left\|\cdot \right\|_{r,\Omega }$ stands indistinctly for the norm of the Hilbertian Sobolev spaces $\mathrm{H}^r(\Omega )$ or $\mathrm{H}^r(\Omega )^d$ with the convention $\mathrm{H}^0(\Omega ):=\mathrm{L}^2(\Omega )$. We also define the Hilbert space ${\mathrm{H}(\mathop{\mathrm{div}}\nolimits ;\Omega )}:=\{\boldsymbol{v}\in \mathrm{L}^2(\Omega )^d:\ \mathop{\mathrm{div}}\nolimits \boldsymbol{v}\in \mathrm{L}^2(\Omega )\}$, whose norm is given by $\left\|\boldsymbol{v}\right\|^2_{\mathop{\mathrm{div}}\nolimits ,\Omega } :=\left\|\boldsymbol{v}\right\|_{0,\Omega }^2+\left\|\mathop{\mathrm{div}}\nolimits \boldsymbol{v}\right\|^2_{0,\Omega }$, and its subspace $\mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega ):=\left\{\boldsymbol{v}\in {\mathrm{H}(\mathop{\mathrm{div}}\nolimits ;\Omega )}: \ \boldsymbol{v}\cdot \boldsymbol{n}=0\,\,\text{on}\,\,\partial \Omega \right\}$.

Finally, $C$ represents a generic constant independent of the discretization parameters, which may take different values at different places.

2. The model problem

We take as our model problem the case of two immiscible fluids within a rigid cavity. Let $\Omega _i$ with $i=1,2$ be the polygonal (in the 2D case) or polyhedral (in the 3D case) Lipschitz domains occupied by each of the fluids. Let $\rho _i$ be the corresponding densities, $\nu _i$ the fluid viscosities, and $c_i$ the acoustic speeds, which we consider all constant, $\rho _i$ and $c_i$ strictly positive, and $\nu _i$ non-negative. We denote by $\boldsymbol{n}_i$ the outward unit normal vectors corresponding to each subdomain. We define $\Omega :=(\overline{\Omega }_1\cup \overline{\Omega }_2)^{\circ }$, $\Gamma :=\partial \Omega _1\cap \partial \Omega _2$, and $\Gamma _i:=\partial \Omega _i\cap \partial \Omega$, $i=1,2$. We assume that each domain $\Omega _i$ as well as $\Omega$ are simply connected (see Figure 1).

We consider small displacements of a compressible viscous fluid at rest neglecting convective terms. The equation of motion derived from the Navier-Stokes equation reads

$$\begin{equation*} \rho _i\ddot{\boldsymbol{U}}_i=2\nu _i\Delta \dot{\boldsymbol{U}}_i-\nabla P_i\quad \mbox{in}\ \Omega _i, \end{equation*}$$

where $\boldsymbol{U}_i$ denotes the fluid displacement and $P_i$ the pressure fluctuation in the domain $\Omega _i$, $i=1,2$. The dot represents derivation with respect to time. (See 8 and 37 for a more detailed derivation.) Moreover, since the fluid is compressible, we consider the isentropic relation

$$\begin{equation} P_i+\rho _ic_i^2\mathop{\mathrm{div}}\nolimits \boldsymbol{U}_i=0\qquad \mbox{in}\ \Omega _i. {\tag{2.1}} \end{equation}$$

Since we are considering irrotational fluids, we assume that $\mathop{\mathbf{curl}}\nolimits \boldsymbol{U}_i=\boldsymbol{0}$. Hence, considering the identity $\Delta \dot{\boldsymbol{U}_i}=\nabla (\mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}_i}) -\mathop{\mathbf{curl}}\nolimits (\mathop{\mathbf{curl}}\nolimits \dot{\boldsymbol{U}_i})$, we conclude that $\Delta \dot{\boldsymbol{U}_i}=\nabla (\mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}_i})$. Then, the equations of our model problem are the following:

$$\begin{align} \rho _1\ddot{\boldsymbol{U}}_1 - 2\nu _1\nabla (\mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}}_1)+\nabla P_1&=\boldsymbol{0}\quad \mbox{in}\ \Omega _1\times (0,T),{\tag{2.2}}\\ P_1 + \rho _1 c_1^2\mathop{\mathrm{div}}\nolimits \boldsymbol{U}_1&=0\quad \mbox{in}\ \Omega _1\times [0,T],{\tag{2.3}}\\ \rho _2\ddot{\boldsymbol{U}}_2- 2\nu _2\nabla (\mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}}_2)+\nabla P_2&=\boldsymbol{0}\quad \mbox{in}\ \Omega _2\times (0,T),{\tag{2.4}}\\ P_2 + \rho _2 c_2^2\mathop{\mathrm{div}}\nolimits \boldsymbol{U}_2&=0\quad \mbox{in}\ \Omega _2\times [0,T],{\tag{2.5}}\\ \boldsymbol{U}_1\cdot \boldsymbol{n}_1+\boldsymbol{U}_2\cdot \boldsymbol{n}_2&=0\quad \mbox{on}\ \Gamma \times [0,T],\tag{2.6}\\ (2\nu _1\mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}}_1+P_1)-(2\nu _2\mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}}_2+P_2)&=0\quad \mbox{on}\ \Gamma \times (0,T),\tag{2.7}\\ \boldsymbol{U}_1\cdot \boldsymbol{n}_1&=0\quad \mbox{on}\hspace{0.1cm}\Gamma _1\times (0,T),\tag{2.8}\\ \boldsymbol{U}_2\cdot \boldsymbol{n}_2&=0\quad \mbox{on}\ \Gamma _2\times (0,T). {\tag{2.9}} \end{align}$$

Let us remark that a similar argument leads to exactly the same equations in 2D.

Multiplying equations 2.2 and 2.4 by a test function $\boldsymbol{v}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega )$, integrating by parts, and using the boundary conditions and the transmission conditions on $\Gamma$, we obtain

$$\begin{equation} \displaystyle \int _{\Omega } \rho \ddot{\boldsymbol{U}}\cdot \boldsymbol{v}+2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}}\mathop{\mathrm{div}}\nolimits \boldsymbol{v}-\int _{\Omega }P\mathop{\mathrm{div}}\nolimits \boldsymbol{v}=0\qquad \forall \boldsymbol{v}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ,\Omega ), {\tag{2.10}} \end{equation}$$

where

$$\begin{equation*} \boldsymbol{U}:= \left\{ \begin{array}{c l} \boldsymbol{U}_1 & \mbox{in}\ \Omega _1,\\\boldsymbol{U}_2 & \mbox{in}\ \Omega _2, \end{array} \right. \hspace{0.2cm} P := \left\{ \begin{array}{c l} P_1 & \mbox{in}\ \Omega _1,\\P_2 & \mbox{in}\ \Omega _2, \end{array} \right. \hspace{0.2cm} \nu := \left\{ \begin{array}{c l} \nu _1 & \mbox{in}\ \Omega _1,\\\nu _2 & \mbox{in}\ \Omega _2, \end{array} \right. \end{equation*}$$

$$\begin{equation*} \rho := \left\{ \begin{array}{c l} \rho _1 & \mbox{in}\ \Omega _1,\\\rho _2 & \mbox{in}\ \Omega _2, \end{array} \right. \hspace{0.2cm}\mbox{and}\hspace{0.2cm} c:= \left\{ \begin{array}{c l} c_1 & \mbox{in}\ \Omega _1,\\c_2 & \mbox{in}\ \Omega _2. \end{array} \right. \end{equation*}$$

Using 2.3 and 2.5 we eliminate $P$ in 2.10 and write

$$\begin{equation} \displaystyle \int _{\Omega } \rho \ddot{\boldsymbol{U}}\cdot \boldsymbol{v}+2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \dot{\boldsymbol{U}}\mathop{\mathrm{div}}\nolimits \boldsymbol{v}+\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{U}\mathop{\mathrm{div}}\nolimits \boldsymbol{v}=0\qquad \forall \boldsymbol{v}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ,\Omega ). {\tag{2.11}} \end{equation}$$

The vibration modes of this problem are complex solutions of the form $\boldsymbol{U}(\boldsymbol{x},t)=e^{\lambda t}\boldsymbol{u}(\boldsymbol{x})$ with $\lambda \in \mathbb{C}$. Looking for this kind of solutions leads to the following quadratic eigenvalue problem.

Problem 1.

Find $\lambda \in \mathbb{C}$ and $\boldsymbol{0}\ne \boldsymbol{u}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega )$ such that

$$\begin{equation*} \displaystyle \lambda ^2\int _{\Omega } \rho \boldsymbol{u}\cdot \overline{\boldsymbol{v}}+2\lambda \int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}+\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}=0\qquad \forall \boldsymbol{v}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega ). \end{equation*}$$

Let us remark that in the absence of viscosity (i.e., $\nu =0$) we are left with the free vibration problem of two inviscid fluids in contact (whose numerical approximation has not been analyzed either). The eigenvalues $\lambda ^2$ of such a problem are negative real numbers (as will be proved below), so that $\lambda$ are purely imaginary, namely, $\lambda =\pm i\omega$ with $\omega$ being the so-called natural vibration frequencies which correspond to periodic in time solutions $\boldsymbol{U}(\boldsymbol{x},t)=e^{-i\omega t}\boldsymbol{u}({\boldsymbol{x}})$ of the time domain problem. This is the reason why, for $\nu =0$, Problem 1 is usually written as follows: Find $\omega >0$ and $\boldsymbol{0}\ne \boldsymbol{u}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ; \Omega )$ such that

$$\begin{equation} \displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= \omega ^2\int _{\Omega }\rho \boldsymbol{u}\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega ). {\tag{2.12}} \end{equation}$$

In the applications, $\nu$ is typically very small. As we will show below, in such a case there are eigenvalues $\lambda$ of Problem 1 that lie close to $\pm i\omega$ with $\omega$ being a natural vibration frequency (i.e., a solution of 2.12). Actually, we will prove below that those $\lambda$ converge to $\pm i\omega$ as $\|\nu \|_{\infty ,\Omega }$ goes to zero. On solving Problem 1, the aim is to compute the eigenvalues $\lambda$ close to the smallest natural vibration frequencies $\omega >0$, which are the most relevant in the applications.

Remark 2.1.

In the case of a homogeneous viscous fluid, $\rho$, $c$, and $\nu$ are constant in the whole $\Omega$. Then, Problem 1 can be written as

$$\begin{equation*} \displaystyle \lambda ^2\int _{\Omega }\rho \boldsymbol{u}\cdot \overline{\boldsymbol{v}}+\frac{2\lambda \nu +\rho c^2}{\rho c^2}\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}=0\qquad \forall \boldsymbol{v}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ,\Omega ). \end{equation*}$$

Hence, in such a case, $(\lambda ,\boldsymbol{u})$ is an eigenpair of Problem 1 if and only if $-\frac{\lambda ^2\rho c^2}{2\lambda \nu +\rho c^2}=\omega ^2$ with $(\omega ,\boldsymbol{u})$ a solution to problem 2.12. Therefore, for a homogeneous viscous fluid, $\lambda$ can be algebraically computed from the solution of 2.12 as follows:

$$\begin{equation*} \displaystyle \lambda =\frac{-\nu \omega ^2\pm \sqrt {\nu ^2\omega ^4-\rho ^2c^4\omega ^2}}{\rho c^2}. \end{equation*}$$

We denote $\mathcal{H}:=\mathrm{L}^2(\Omega )^d$ endowed with the weighted inner product

$$\begin{equation*} \displaystyle (\boldsymbol{v},\boldsymbol{w})_{\mathcal{H}}:=\int _{\Omega }\rho \boldsymbol{v}\cdot \overline{\boldsymbol{w}}, \end{equation*}$$

and we denote $\mathcal{V}:=\mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega )$ with the inner product

$$\begin{equation*} \displaystyle (\boldsymbol{v},\boldsymbol{w})_{\mathcal{V}} :=\int _{\Omega }\rho \boldsymbol{v}\cdot \overline{\boldsymbol{w}}+\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{v}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{w}}. \end{equation*}$$

Notice that the inner products in $\mathcal{H}$ and $\mathcal{V}$ induce norms $\|\cdot \|_{\mathcal{H}}$ and $\|\cdot \|_{\mathcal{V}}$ on each of these spaces equivalent to the classical $\mathrm{L}^2(\Omega )^d$ and $\mathrm{H}(\mathop{\mathrm{div}}\nolimits ;\Omega )$ norms, respectively. Therefore, when it might be convenient, we will use these classical norms.

Clearly $\lambda =0$ is an eigenvalue of Problem 1 with associated eigenspace

$$\begin{equation*} \mathcal{K}=\mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ^0,\Omega ):= \left\{\boldsymbol{v}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega ):\ \mathop{\mathrm{div}}\nolimits \boldsymbol{v}=0\,\,\text{in}\,\,\Omega \right\}. \end{equation*}$$

We define

$$\begin{equation*} \mathcal{G}:=\mathcal{K}^{\bot _{\mathcal{V}}}=\{\boldsymbol{v}\in \mathcal{V}: (\boldsymbol{v},\boldsymbol{w})_{\mathcal{V}}=0\qquad \forall \boldsymbol{w}\in \mathcal{K}\}. \end{equation*}$$

Since $\mathcal{K}$ is a closed subspace of $\mathcal{V}$, clearly $\mathcal{V}=\mathcal{G}\oplus \mathcal{K}.$ Notice that $\mathcal{G}$ and $\mathcal{K}$ are also orthogonal in the $\mathcal{H}$ inner product. Hence,

$$\begin{equation*} \mathcal{G}=\{\boldsymbol{v}\in \mathcal{V}: (\boldsymbol{v},\boldsymbol{w})_{\mathcal{H}}=0\qquad \forall \boldsymbol{w}\in \mathcal{K}\}. \end{equation*}$$

The following result brings a characterization of the space $\mathcal{G}$.

Lemma 2.1.

There holds

$$\begin{equation*} \mathcal{G}=\displaystyle \frac{1}{\rho }\nabla (\mathrm{H}^1(\Omega ))\cap \mathcal{V}. \end{equation*}$$

Proof.

We will prove this result by checking the double inclusion. Let $\boldsymbol{v}\in \mathcal{G}.$ Then, for all $\boldsymbol{\psi }\in \mathcal{D}(\Omega )^d$, since $\mathop{\mathbf{curl}}\nolimits \boldsymbol{\psi }\in \mathcal{K}$, we have

$$\begin{equation*} \displaystyle 0=\int _{\Omega }\rho \mathop{\mathbf{curl}}\nolimits \boldsymbol{\psi }\cdot \overline{\boldsymbol{v}}= \int _{\Omega }\boldsymbol{\psi }\cdot \mathop{\mathbf{curl}}\nolimits \left(\rho \overline{\boldsymbol{v}}\right). \end{equation*}$$

Thus, $\mathop{\mathbf{curl}}\nolimits \left(\rho \boldsymbol{v}\right)=\boldsymbol{0}$ in $\Omega$. Since $\Omega$ is simply connected, this implies that there exists $\varphi \in \mathrm{H}^1(\Omega )$ such that $\rho \boldsymbol{v}=\nabla \varphi$. Hence, $\boldsymbol{v}\in \frac{1}{\rho }\nabla (\mathrm{H}^1(\Omega ))\cap \mathcal{V}$. Conversely, let $\boldsymbol{v}\in \frac{1}{\rho }\nabla (\mathrm{H}^1(\Omega ))\cap \mathcal{V}$ and $\boldsymbol{w}\in \mathcal{K}$. Let $\varphi \in \mathrm{H}^1(\Omega )$ be such that $\boldsymbol{v}=\frac{1}{\rho }\nabla \varphi$. Then,

$$\begin{eqnarray*} \displaystyle (\boldsymbol{v},\boldsymbol{w})_{\mathcal{H}}=\displaystyle \int _{\Omega }\rho \left(\frac{1}{\rho }\nabla \varphi \right) \cdot \overline{\boldsymbol{w}}=-\displaystyle \int _{\Omega }\varphi \mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{w}}+\int _{\partial \Omega }\varphi (\overline{\boldsymbol{w}}\cdot \boldsymbol{n})=0. \end{eqnarray*}$$

Therefore, $\boldsymbol{v}\in \mathcal{G}$. The proof is complete.

■

In what follows we prove additional regularity for the functions in $\mathcal{G}$ on each subdomain. From now on, $s$ will denote a positive number such that the following lemma holds true.

Lemma 2.2.

There exists $s>0$ (with $s$ depending on $\rho$, $\Omega _1$, and $\Omega _2$) such that, for all $\boldsymbol{v}\in \mathcal{G}$, $\boldsymbol{v}\in \mathrm{H}^{s}(\Omega _1\cup \Omega _2)^d$ and

$$\begin{equation} \|\boldsymbol{v}\|_{s,\Omega _1}+\|\boldsymbol{v}\|_{s,\Omega _2}\leq C\|\mathop{\mathrm{div}}\nolimits \boldsymbol{v}\|_{0,\Omega }, {\tag{2.13}} \end{equation}$$

where $C$ is a positive constant independent of $\boldsymbol{v}$.

Proof.

According to Lemma 2.1, there exists $\varphi \in \mathrm{H}^1(\Omega )$ such that $\boldsymbol{v}=\frac{1}{\rho }\nabla \varphi$. Consequently, $\varphi \in \mathrm{H}^1(\Omega )/ \mathbb{C}$ is the unique solution of the following well-posed Neumann problem:

$$\begin{align*} \displaystyle \mathop{\mathrm{div}}\nolimits \left(\frac{1}{\rho }\nabla \varphi \right)&=\mathop{\mathrm{div}}\nolimits \boldsymbol{v}\quad \mbox{in}\ \Omega ,\\ \displaystyle \frac{1}{\rho }\frac{\partial \varphi }{\partial \boldsymbol{n}}&=0\hspace{0.2cm}\qquad \mbox{on}\ \partial \Omega . \end{align*}$$

Hence, in the 3D case, the result follows from 34, Lemma 2.20, while in the 2D case, it follows by applying 34, Lemma 4.3. (See 36 for more details.)

■

Remark 2.2.

The above lemma establishes the existence of a regularity exponent $s>0$ that will play a role in the error estimates of the numerical method proposed in this paper. In spite of the fact that we refer to 34 in the proof of this lemma, the value of $s$ that arises from this reference is far from being optimal, since it is valid for global regularity in $\mathrm{H}^{s}(\Omega )^d$ and for any arbitrary geometrical setting of the subdomains $\Omega _1$ and $\Omega _2$. In most of the applications, the subdomains at which $\rho$ is constant are similar to those shown in Figure 1. In such a case, a detailed analysis of how this exponent depends on the geometry of the domain and on the coefficient $\rho$ can be found in 25 for the 2D case and in 24 for 3D problems (see also 933). Let us remark that, although the analysis of these references is for problems with Dirichlet boundary conditions, similar results hold true for Neumann boundary conditions as in our case (see 25, Remark III.5.2). In particular, for instance, Lemma 2.2 holds true for $s=1$ in the example reported as Test 1 in Section 5 (see Figure 2).

From the physical point of view, the time domain problem 2.11 is dissipative in the sense that its solution should decay as $t$ increases. The latter happens if and only if the so-called decay rate, $\mbox{Re}(\lambda )$, is negative. The following result shows that this is the case in our formulation.

Lemma 2.3.

Let $(\lambda ,\boldsymbol{u})\in \mathbb{C}\times \mathcal{V}$ be a solution of Problem 1. If $\lambda \ne 0$, then $\text{Re}(\lambda )<0.$

Proof.

The proof is similar to that of Lemma 2.1 of 4.

■

Remark 2.3.

Any eigenpair $(\lambda ,\boldsymbol{u})$ of Problem 1 satisfies

$$\begin{equation*} \displaystyle \lambda ^2\int _{\Omega }\rho \boldsymbol{u}\cdot \overline{\boldsymbol{v}} +\int _{\Omega }(2\lambda \nu +\rho c^2)\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}=0\qquad \forall \boldsymbol{v}\in \mathcal{V}. \end{equation*}$$

Since the coefficients are constant in each subdomain, if $2\lambda \nu +\rho c^2\ne 0$ in $\Omega _i$, by testing with $\boldsymbol{v}\in \mathcal{D}(\Omega _i)^d$ we obtain that $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}|_{\Omega _i}\in \mathrm{H}^1(\Omega _i)$, $i=1,2$. On the other hand, if $2\lambda \nu +\rho c^2=0$ in $\Omega _i$ ($i=1$ or $2$), then, for $\lambda \ne 0$, by testing again with $\boldsymbol{v}\in \mathcal{D}(\Omega _i)^d$, we obtain that $\boldsymbol{u}=\boldsymbol{0}$ in $\Omega _i$. Thus, in any case, $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}|_{\Omega _i}\in \mathrm{H}^1(\Omega _i)$, $i=1,2.$

For the theoretical analysis it is convenient to transform Problem 1 into a linear eigenvalue problem. With this aim we introduce the new variable $\widehat{\boldsymbol{u}}:=\lambda \boldsymbol{u}$, as usual in quadratic problems, and the space $\widetilde{\mathcal{V}}:=\mathcal{V}\times \mathcal{H}$ endowed with the corresponding product norm, which carry us to the following.

Problem 2.

Find $\lambda \in \mathbb{C}$ and $\boldsymbol{0}\ne (\boldsymbol{u},\widehat{\boldsymbol{u}})\in \widetilde{\mathcal{V}}$ such that

$$\begin{align} \displaystyle &\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= \lambda \left(-2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}- \int _{\Omega }\rho \widehat{\boldsymbol{u}}\cdot \overline{\boldsymbol{v}}\right)\qquad \forall \boldsymbol{v}\in \mathcal{V}, {\tag{2.14}}\\ \displaystyle &\int _{\Omega }\rho \widehat{\boldsymbol{u}}\cdot \overline{\widehat{\boldsymbol{v}}}= \lambda \int _{\Omega }\rho \boldsymbol{u}\cdot \overline{\widehat{\boldsymbol{v}}}\qquad \forall \widehat{\boldsymbol{v}}\in \mathcal{H}. {\tag{2.15}} \end{align}$$

We observe that $\lambda =0$ is an eigenvalue of Problem 2 and its associated eigenspace is $\widetilde{\mathcal{K}}:=\mathcal{K}\times \{0\}$. Let $\widetilde{\mathcal{G}}$ be the orthogonal complement of $\widetilde{\mathcal{K}}$ in $\mathcal{V}\times \mathcal{H}$. Notice that $\widetilde{\mathcal{G}}=\mathcal{G}\times \mathcal{H}$.

We introduce the sesquilinear continuous form $a:\mathcal{V}\times \mathcal{V}\rightarrow \mathbb{C}$ defined by

$$\begin{equation*} a(\boldsymbol{u},\boldsymbol{v}):=\displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}, \end{equation*}$$

and the sesquilinear continuous forms $\widetilde{a}, \widetilde{b}:\widetilde{\mathcal{V}}\rightarrow \widetilde{\mathcal{V}}$ defined as follows:

$$\begin{equation*} \widetilde{a}((\boldsymbol{u},\widehat{\boldsymbol{u}}),(\boldsymbol{v}, \widehat{\boldsymbol{v}})):=\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}+\int _{\Omega }\rho \widehat{\boldsymbol{u}}\cdot \overline{\widehat{\boldsymbol{v}}}, \end{equation*}$$

$$\begin{equation*} \widetilde{b}((\boldsymbol{u},\widehat{\boldsymbol{u}}),(\boldsymbol{v}, \widehat{\boldsymbol{v}})):=\displaystyle -2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}-\int _{\Omega }\rho \widehat{\boldsymbol{u}}\cdot \overline{\boldsymbol{v}}+ \int _{\Omega }\rho \boldsymbol{u}\cdot \overline{\widehat{\boldsymbol{v}}}. \end{equation*}$$ In what follows we prove that $a(\cdot ,\cdot )$ and $\widetilde{a}(\cdot ,\cdot )$ are elliptic in $\mathcal{G}$ and $\widetilde{\mathcal{G}}$, respectively.

Lemma 2.4.

The sesquilinear form $a: \mathcal{G}\times \mathcal{G}\rightarrow \mathbb{C}$ is $\mathcal{G}$-elliptic and, consequently, $\widetilde{a}:\widetilde{\mathcal{G}}\times \widetilde{\mathcal{G}}\rightarrow \mathbb{C}$ is $\widetilde{\mathcal{G}}$-elliptic.

Proof.

The proof is a direct consequence of Lemma 2.2.

■

Let $\boldsymbol{T}:\widetilde{\mathcal{V}}\rightarrow \widetilde{\mathcal{V}}$ be the bounded linear operator defined by $\boldsymbol{T}(\boldsymbol{f},\boldsymbol{g}):=(\boldsymbol{u},\widehat{\boldsymbol{u}})\in \widetilde{\mathcal{G}}$, where $(\boldsymbol{u},\widehat{\boldsymbol{u}})$ is the unique solution of the following problem:

$$\begin{equation*} \widetilde{a}((\boldsymbol{u},\widehat{\boldsymbol{u}}),(\boldsymbol{v},\widehat{\boldsymbol{v}}))= \widetilde{b}((\boldsymbol{f},\boldsymbol{g}),(\boldsymbol{v}, \widehat{\boldsymbol{v}})) \qquad \forall (\boldsymbol{v},\widehat{\boldsymbol{v}})\in \widetilde{\mathcal{G}}. \end{equation*}$$

It is easy to check that

$$\begin{equation} \widehat{\boldsymbol{u}}=\boldsymbol{f}\quad \text{in}\ \Omega , {\tag{2.16}} \end{equation}$$

and

$$\begin{equation} \displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}=-2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}} -\int _{\Omega }\rho \boldsymbol{g}\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. {\tag{2.17}} \end{equation}$$

As a consequence of the above equalities, we have that $\mu =0$ is an eigenvalue of $\boldsymbol{T}$ with associated eigenspace $\{\boldsymbol{0}\}\times \mathcal{G}^{\bot _{\mathcal{H}}}$, which is non-trivial since $\mathcal{G}^{\bot _{\mathcal{H}}}\supset \mathcal{K}$. The following lemma shows that the non-zero eigenvalues of $\boldsymbol{T}$ are exactly the reciprocals of the non-zero eigenvalues of Problem 2 with the same corresponding eigenfunctions.

Lemma 2.5.

There holds that $(\mu ,(\boldsymbol{u},\widehat{\boldsymbol{u}}))$ is an eigenpair of $\boldsymbol{T}$ $(i.e.,\,\,\boldsymbol{T}(\boldsymbol{u},\widehat{\boldsymbol{u}})=\mu (\boldsymbol{u},\widehat{\boldsymbol{u}}))$ with $\mu \ne 0$ if and only if $(\lambda ,(\boldsymbol{u},\widehat{\boldsymbol{u}}))$ is a solution of Problem 2 with $\lambda =1/\mu \ne 0.$

Proof.

The proof is similar to that of Lemma 2.3 from 4.

■

3. Spectral characterization

The goal of this section is to characterize the spectrum of the solution operator $\boldsymbol{T}$. Since the inclusion $\mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ;\Omega )\hookrightarrow \mathrm{L}^2(\Omega )^d$ is not compact, it is easy to check from 2.16 that $\boldsymbol{T}$ is not compact either. However, we will show that the essential spectrum has to lie in a small region of the complex plane, well-separated from the isolated eigenvalues which, according to Lemma 2.5, correspond to the solutions of Problem 2. With this aim, we will resort to the theory described in 21 to appropriately decompose $\boldsymbol{T}$. Let $\boldsymbol{T}_1, \boldsymbol{T}_2: \mathcal{G}\rightarrow \mathcal{G}$ be the operators given by

$$\begin{align} \displaystyle &\boldsymbol{T}_1 \boldsymbol{f}=\boldsymbol{u}_1\in \mathcal{G}:\quad a(\boldsymbol{u}_1,\boldsymbol{v})=2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G},{\tag{3.1}}\\ \displaystyle &\boldsymbol{T}_2 \boldsymbol{g}=\boldsymbol{u}_2\in \mathcal{G}:\quad a(\boldsymbol{u}_2,\boldsymbol{v})=\int _{\Omega }\rho \boldsymbol{g}\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. {\tag{3.2}} \end{align}$$

It is easy to check that these operators are self-adjoint with respect to $a(\cdot ,\cdot )$. Moreover $\boldsymbol{T}_1$ is non-negative and $\boldsymbol{T}_2$ is positive with respect to $a(\cdot ,\cdot )$ (namely, $a(\boldsymbol{T}_1\boldsymbol{v},\boldsymbol{v})\geq 0$ $\forall \boldsymbol{v}\in \mathcal{G}$ and $a(\boldsymbol{T}_2\boldsymbol{v},\boldsymbol{v})>0$ $\forall \boldsymbol{v}\in \mathcal{G}$, $\boldsymbol{v}\ne \boldsymbol{0}$). Moreover, we have the following result.

Lemma 3.1.

The operator $\boldsymbol{T}_2:\mathcal{G}\to \mathcal{G}$ is compact.

Proof.

Since $a(\cdot ,\cdot )$ is $\mathcal{G}$-elliptic (cf. Lemma 2.4), applying Lax-Milgram’s Lemma, we know that problem 3.2 is well-posed and has a unique solution $\boldsymbol{u}_2\in \mathcal{G}$. Moreover, according to Lemma 2.2, we know that there exists $s>0$ such that $\boldsymbol{u}_2\in \mathrm{H}^s(\Omega _1\cup \Omega _2)^d$. On the other hand, notice that 3.2 also holds for $\boldsymbol{v}\in \mathcal{K}$, since in such a case $a(\boldsymbol{u}_2,\boldsymbol{v})=0=\int _{\Omega }\rho \boldsymbol{g}\cdot \overline{\boldsymbol{v}}$ for $\boldsymbol{g}\in \mathcal{G}$. Hence, since $\mathcal{V}=\mathcal{G}\oplus \mathcal{K}$, we have that

$$\begin{equation*} \displaystyle a(\boldsymbol{u}_2,\boldsymbol{v})=\int _{\Omega }\rho \boldsymbol{g}\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{V}. \end{equation*}$$

Then, by testing this equation with $\boldsymbol{v}\in \mathcal{D}(\Omega )^d\subset \mathcal{V}$, we have that $-\nabla (\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2)=\rho \boldsymbol{g}$ in $\Omega$, so that $\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\in \mathrm{H}^1(\Omega )$. Therefore, since $\rho$ and $c$ are positive constants in each subdomain $\Omega _1$ and $\Omega _2$, we have that $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2|_{\Omega _i}\in \mathrm{H}^1(\Omega _i)$, $i=1,2.$ Since the inclusions $\{v\in \mathrm{L}^2(\Omega ): v|_{\Omega _i}\in \mathrm{H}^1(\Omega _i),\,i=1,2\}\subset \mathrm{L}^2(\Omega )$ and $\mathrm{H}^s(\Omega _1\cup \Omega _2)^{d}\subset \mathrm{L}^2(\Omega )^d$ are compact, we derive that $\boldsymbol{T}_2$ is compact too.

■

The operator $\boldsymbol{T}$ can be written in terms of the operators $\boldsymbol{T}_1$ and $\boldsymbol{T}_2$ given above as follows:

$$\begin{equation*} \boldsymbol{T}=\begin{pmatrix} -\boldsymbol{T}_1&-\boldsymbol{T}_2\\ \boldsymbol{I}&\boldsymbol{0}\end{pmatrix}. \end{equation*}$$

Moreover, by defining as in 21 the operators

$$\begin{equation*} \boldsymbol{S}:=\begin{pmatrix} \boldsymbol{I}&\boldsymbol{0}\\ \boldsymbol{0}&\boldsymbol{T}_2^{1/2} \end{pmatrix} \qquad \mbox{and}\qquad \boldsymbol{H}:=\begin{pmatrix} -\boldsymbol{T}_1&-\boldsymbol{T}_2^{1/2}\\ \boldsymbol{T}_2^{1/2}&\boldsymbol{0}\end{pmatrix}, \end{equation*}$$

we have that $\boldsymbol{S}\boldsymbol{T}=\boldsymbol{H}\boldsymbol{S}$. We note that the eigenvalues of $\boldsymbol{T}$ and $\boldsymbol{H}$ and their algebraic multiplicities coincide. Moreover the corresponding Jordan chains have the same length. In fact, let $\{\boldsymbol{x}_k\}_{k=1}^r$ be a Jordan chain associated with the eigenvalue $\mu$ of $\boldsymbol{T}$. Then, using the identities above, we observe that

$$\begin{equation*} \boldsymbol{H}\boldsymbol{S}\boldsymbol{x}_k=\boldsymbol{S}\boldsymbol{T}\boldsymbol{x}_k=\boldsymbol{S}(\mu \boldsymbol{x}_k+\boldsymbol{x}_{k-1})=\mu \boldsymbol{S}\boldsymbol{x}_k+\boldsymbol{S}\boldsymbol{x}_{k-1}, \quad k=1,\ldots ,r. \end{equation*}$$

This shows that $\{\boldsymbol{S}\boldsymbol{x}_k\}_{k=1}^r$ is a Jordan chain of $\boldsymbol{H}$ of the same length. Actually, the whole spectra of $\boldsymbol{T}$ and $\boldsymbol{H}$ coincide as is shown in the following result, which has been proved in Lemma 3.2 of 4.

Lemma 3.2.

There holds

$$\begin{equation*} \mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T})=\mathop{\mathrm{Sp}}\nolimits (\boldsymbol{H}). \end{equation*}$$

Moreover, $\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{T})=\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{H})$, too.

The operator $\boldsymbol{H}$ can be written as the sum of a self-adjoint operator $\boldsymbol{B}$ and a compact operator $\boldsymbol{C}$:

$$\begin{equation*} \boldsymbol{H}=\boldsymbol{B}+\boldsymbol{C}\quad \mbox{with}\quad \boldsymbol{B}:=\begin{pmatrix} -\boldsymbol{T}_1&\boldsymbol{0}\\ \boldsymbol{0}&\boldsymbol{0}\end{pmatrix}\quad \mbox{and}\quad \boldsymbol{C}:=\begin{pmatrix} \boldsymbol{0}&-\boldsymbol{T}_2^{1/2}\\ \boldsymbol{T}_2^{1/2}&\boldsymbol{0}\end{pmatrix}. \end{equation*}$$

Then, applying the classical Weyl’s Theorem (see 38), we have that $\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{H})=\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{B})$ and the rest of the spectrum $\mathop{\mathrm{Sp}_\mathrm{disc}}\nolimits (\boldsymbol{H}):=\mathop{\mathrm{Sp}}\nolimits (\boldsymbol{H})\backslash \mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{H})$ consists of isolated eigenvalues with finite algebraic multiplicity. Moreover, $\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{B})=\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (-\boldsymbol{T}_1)\cup \{0\}$.

Our next goal is to show that the essential spectrum of $\boldsymbol{T}_1$ must lie in a small region of the complex plane. Actually, we will localize the whole spectrum of $\boldsymbol{T}_1$. With this aim, we analyze separately for which values $\mu \in \mathbb{C}$ the operator $(\mu \boldsymbol{I}-\boldsymbol{T}_1)$ is not necessarily one-to-one and for which values it is not necessarily onto.

•

If $(\mu \boldsymbol{I}- \boldsymbol{T}_1)$ is not one-to-one, then there exists $\boldsymbol{f}\in \mathcal{G}$, $\boldsymbol{f}\neq \boldsymbol{0}$, such that $\boldsymbol{T}_1\boldsymbol{f}=\mu \boldsymbol{f}$, namely,$$\begin{equation*} \displaystyle \mu \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{f}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= 2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. \end{equation*}$$

Then, testing with $\boldsymbol{v}=\boldsymbol{f}$ and using that in each subdomain the coefficients $\rho$ and $c$ are positive, we deduce that$$\begin{equation*} \mu =\frac{2\int _{\Omega }\nu \left|\mathop{\mathrm{div}}\nolimits \boldsymbol{f}\right|^2}{\int _{\Omega }\rho c^2\left|\mathop{\mathrm{div}}\nolimits \boldsymbol{f}\right|^2} \end{equation*}$$

(we recall that for $\boldsymbol{0}\ne \boldsymbol{f}\in \mathcal{G}$, $\int _{\Omega }|\mathop{\mathrm{div}}\nolimits \boldsymbol{f}|^2>0$ because of Lemma 2.2). Hence,$$\begin{equation*} \mu \in \displaystyle \left[\frac{2\min _{\Omega }\{\nu \}}{\max _{\Omega }\{\rho c^2\}},\frac{2\max _{\Omega }\{\nu \}}{\min _{\Omega }\{\rho c^2\}}\right]. \end{equation*}$$

•

On the other hand, $(\mu \boldsymbol{I}- \boldsymbol{T}_1)$ is onto if and only if for any $\boldsymbol{g}\in \mathcal{G}$ there exists $\boldsymbol{f}\in \mathcal{G}$ such that $\boldsymbol{T}_1\boldsymbol{f}=\mu \boldsymbol{f}-\boldsymbol{g}$, which from 3.1 reads$$\begin{equation*} \displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{g}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= \int _{\Omega }(-2\nu +\mu \rho c^2)\mathop{\mathrm{div}}\nolimits \boldsymbol{f}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. \end{equation*}$$

By writing $\mu =\alpha +\beta i$ with $\alpha ,\beta \in \mathbb{R}$, the equation above reads:$$\begin{equation*} \int _{\Omega }(-2\nu +\alpha \rho c^2+\rho c^2\beta i)\mathop{\mathrm{div}}\nolimits \boldsymbol{f}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= \displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{g}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. \end{equation*}$$

We observe that for all $\beta \ne 0$ the problem above has a solution and hence the operator $(\mu \boldsymbol{I}-\boldsymbol{T}_1)$ is onto. On the other hand, if $\beta =0$, then $\mu$ has to be real. In such a case, the operator $\boldsymbol{T}_1$ will still be onto when $(-2\nu +\mu \rho c^2)$ has the same sign in the whole domain $\Omega$. This happens at least in two cases:

(i)

when $\mu >\frac{2\max _{\Omega }\{\nu \}}{\min _{\Omega }\{\rho c^2\}}$, in which case $-2\nu +\mu \rho c^2>0$,

(ii)

when $\mu <\frac{2\min _{\Omega }\{\nu \}}{\max _{\Omega }\{\rho c^2\}}$, in which case $-2\nu +\mu \rho c^2<0$.

Therefore, if $(\mu \boldsymbol{I}-\boldsymbol{T}_1)$ is not onto, then $\mu \in \left[\frac{2\min _{\Omega }\{\nu \}}{\max _{\Omega }\{\rho c^2\}},\frac{2\max _{\Omega }\{\nu \}}{\min _{\Omega }\{\rho c^2\}}\right]$, too.

Now we are in a position to write the following spectral characterization of the solution operator $\boldsymbol{T}$.

Theorem 3.1.

The spectrum of $\boldsymbol{T}$ consists of

$$\begin{equation*} \mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{T})=\mathop{\mathrm{Sp}}\nolimits (-\boldsymbol{T}_1)\cup \{0\} \end{equation*}$$

with

$$\begin{equation*} \mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T}_1)\subset \left[\frac{2\min _{\Omega }\{\nu \}}{\max _{\Omega }\{\rho c^2\}},\frac{2\max _{\Omega }\{\nu \}}{\min _{\Omega }\{\rho c^2\}}\right] \end{equation*}$$

and $\mathop{\mathrm{Sp}_\mathrm{disc}}\nolimits (\boldsymbol{T}):=\mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T})\setminus \mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{T})$, which is a set of isolated eigenvalues of finite algebraic multiplicity.

Proof.

As a consequence of the classical Weyl’s Theorem (see 38) and Lemma 3.2,

$$\begin{equation*} \mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{T})=\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{H})=\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (\boldsymbol{B})=\mathop{\mathrm{Sp}_\mathrm{ess}}\nolimits (-\boldsymbol{T}_1)\cup \{0\}, \end{equation*}$$

whereas the inclusion follows from the above analysis.

■

In what follows, we will show that for $\nu$ small enough some of the eigenvalues of $\boldsymbol{T}$ are well-separated from its essential spectrum. To this end, given $\boldsymbol{f}\in \mathcal{G}$, by testing 3.1 with $\boldsymbol{v}=\boldsymbol{u}_1\in \mathcal{G}$ and using the definition of $a(\cdot ,\cdot )$, we have that

$$\begin{align*} \displaystyle \min _{\Omega }\{\rho c^2\}\|\boldsymbol{u}_1\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }^2\leq a(\boldsymbol{u}_1,\boldsymbol{u}_1)&\leq 2\|\nu \|_{\infty ,\Omega }\|\mathop{\mathrm{div}}\nolimits \boldsymbol{f}\|_{0,\Omega }\|\boldsymbol{u}_1\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }. \end{align*}$$

Therefore $\|\boldsymbol{T}_1\|_{\mathcal{L}(\mathcal{G}\times \mathcal{G})}\rightarrow 0$ as $\|\nu \|_{\infty ,\Omega }$ goes to zero. Consequently, $\boldsymbol{H}$ converges in norm to the operator

$$\begin{equation*} \boldsymbol{H}_0:=\begin{pmatrix} \boldsymbol{0}&-\boldsymbol{T}_2^{1/2}\\ \boldsymbol{T}_2^{1/2}&\boldsymbol{0}\end{pmatrix} \end{equation*}$$

as $\|\nu \|_{\infty ,\Omega }$ goes to zero. Thus, from the classical spectral approximation theory (see 19), the isolated eigenvalues of $\boldsymbol{H}$ converge to those of $\boldsymbol{H}_0$.

Since the isolated eigenvalues of $\boldsymbol{H}$ and $\boldsymbol{T}$ coincide (cf. Lemma 3.2), in order to localize those of $\boldsymbol{T}$, we begin by characterizing those of $\boldsymbol{H}_0$. Let $\mu$ be an isolated eigenvalue of $\boldsymbol{H}_0$ and $(\boldsymbol{u},\widehat{\boldsymbol{u}})\in \mathcal{G}\times \mathcal{G}$ the corresponding eigenfunction. It is easy to check that

$$\begin{equation} \boldsymbol{H}_0\begin{pmatrix} \boldsymbol{u}\\ \widehat{\boldsymbol{u}}\end{pmatrix}=\mu \begin{pmatrix} \boldsymbol{u}\\ \widehat{\boldsymbol{u}}\end{pmatrix}\quad \Longleftrightarrow \quad \boldsymbol{T}_2\boldsymbol{u}=-\mu ^2\boldsymbol{u}\quad \text{and}\quad \boldsymbol{T}_2^{1/2}\boldsymbol{u}=\mu \widehat{\boldsymbol{u}}. {\tag{3.3}} \end{equation}$$

Since $\boldsymbol{T}_2$ is compact, self-adjoint, and positive, its spectrum consists of a sequence of positive eigenvalues that converge to zero and zero itself. Notice that the spectrum of $\boldsymbol{T}_2$ is related with the solution of the eigenvalue problem 2.12. In fact, this problem has $0$ as an eigenvalue with corresponding eigenspace $\mathcal{K}$. The rest of the eigenvalues $\omega ^2$ are strictly positive and the corresponding eigenfunctions $\boldsymbol{u}\in \mathcal{K}^{\bot _{\mathcal{V}}}=:\mathcal{G}$, so that they are also solutions of the following problem: Find $\omega >0$ and $\boldsymbol{u}\in \mathcal{G}$ such that

$$\begin{equation*} \displaystyle a(\boldsymbol{u},\boldsymbol{v})=\omega ^2\int _{\Omega }\rho \boldsymbol{u}\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. \end{equation*}$$

Clearly $(\omega ^2,\boldsymbol{u})$ is an eigenpair of the above problem with $\omega >0$ if and only if $\boldsymbol{T}_2\boldsymbol{u}=\frac{1}{\omega ^2}\boldsymbol{u}$. Thus, by virtue of 3.3, we have that the eigenvalues of $\boldsymbol{H}_0$ are given by $\pm i/\omega$ and hence they are purely imaginary.

Now we are in a position to establish the following result.

Theorem 3.2.

For each isolated eigenvalue $\pm i/\omega$ of $\boldsymbol{H}_0$ of algebraic multiplicity $m$, let $r>0$ be such that the disc $D_r:=\{ z\in \mathbb{C}:\hspace{0.2cm} |z\mp i/\omega |< r\}$ intersects $\mathop{\mathrm{Sp}}\nolimits (\boldsymbol{H}_0)$ only in $\pm i/\omega$. Then, there exists $\delta >0$ such that if $\|\nu \|_{\infty ,\Omega }<\delta$, there exist $m$ eigenvalues of $\boldsymbol{T}$, $\mu _1,\ldots , \mu _m$ (repeated according to their respective algebraic multiplicities), lying in the disc $D_r$. Moreover, $\mu _1,\ldots ,\mu _m\rightarrow \frac{i}{\omega }$ as $\|\nu \|_{\infty ,\Omega }$ goes to zero.

As claimed above, the eigenvalues of $\boldsymbol{T}$ that are relevant in the applications are those which are close to $\pm i/\omega$ for the smallest positive vibration frequencies $\omega$ of 2.12. According to the above theorem, these eigenvalues are well-separated from the real axis and, hence, from the essential spectrum of $\boldsymbol{T}$ (cf. Theorem 3.1).

4. Spectral approximation

In this section, we propose and analyze a finite element method to approximate the solutions of Problem 1. To this end, we introduce appropriate discrete spaces. Let $\{\mathcal{T}_h(\Omega )\}_{h>0}$ be a family of regular partitions of $\Omega$ such that $\mathcal{T}_h(\Omega _i):=\{T\in \mathcal{T}_h: T\subset \overline{\Omega }_i\}$ are partitions of $\Omega _i$, $i=1,2$. We introduce the lowest-order Raviart-Thomas finite element space:

$$\begin{equation*} \mathcal{V}_h:=\{\boldsymbol{v}\in \mathcal{V}: \boldsymbol{v}|_T(\boldsymbol{x})=\boldsymbol{a}+b\boldsymbol{x},\ \boldsymbol{a}\in \mathbb{R}^d,\ b\in \mathbb{R},\ \boldsymbol{x}\in T\}. \end{equation*}$$

The discretization of Problem 1 reads as follows.

Problem 3.

Find $\lambda _h\in \mathbb{C}$ and $\boldsymbol{0}\neq \boldsymbol{u}_h\in \mathcal{V}_h$ such that

$$\begin{equation*} \displaystyle \lambda _h^2\int _{\Omega }\rho \boldsymbol{u}_h\cdot \overline{\boldsymbol{v}}_h+ 2\lambda _h\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{u}_h\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h+\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_h\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h=0\qquad \forall \boldsymbol{v}_h\in \mathcal{V}_h. \end{equation*}$$

We proceed as we did in the continuous case and introduce a new discrete variable $\widehat{\boldsymbol{u}}_h:=\lambda _h\boldsymbol{u}_h$ to rewrite the problem above in the following equivalent form.

Problem 4.

Find $\lambda _h\in \mathbb{C}$ and $\boldsymbol{0}\ne (\boldsymbol{u}_h, \widehat{\boldsymbol{u}}_h)\in \mathcal{V}_h\times \mathcal{V}_h$ such that

$$\begin{align*} \displaystyle &\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_h\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h=\lambda _h\left(-2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{u}_h \mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h-\int _{\Omega }\rho \widehat{\boldsymbol{u}}_h\cdot \overline{\boldsymbol{v}}_h\right)\qquad \forall \boldsymbol{v}_h\in \mathcal{V}_h, \\ \displaystyle &\int _{\Omega }\rho \widehat{\boldsymbol{u}}_h\cdot \overline{\widehat{\boldsymbol{v}}}_h= \lambda _h\int _{\Omega }\rho \boldsymbol{u}_h\cdot \overline{\widehat{\boldsymbol{v}}}_h\qquad \forall \widehat{\boldsymbol{v}}_h\in \mathcal{V}_h. \end{align*}$$

We observe that $\lambda _h=0$ is an eigenvalue of this problem and its associated eigenspace is $\widetilde{\mathcal{K}}_h:=\mathcal{K}_h\times \{0\}$ with $\mathcal{K}_h:=\mathcal{K}\cap \mathcal{V}_h$ the eigenspace of $\lambda _h=0$ in Problem 3. At the beginning of Section 5, we will show that Problem 4 is well-posed, in the sense that it is equivalent to a generalized matrix eigenvalue problem with a symmetric positive definite right-hand side matrix.

We introduce the well-known Raviart-Thomas interpolation operator, $\Pi _h: \mathcal{V}\cap \mathrm{H}^r(\Omega _1\cup \Omega _2)^d\rightarrow \mathcal{V}_h$, $r\in (0,1]$ (see 29), for which the approximation result

$$\begin{equation} \displaystyle \|\boldsymbol{v}-\Pi _h\boldsymbol{v}\|_{0,\Omega }\leq Ch^r(\|\boldsymbol{v}\|_{r,\Omega _1}+\|\boldsymbol{v}\|_{r,\Omega _2}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{v}\|_{0,\Omega }) {\tag{4.1}} \end{equation}$$

and the commuting diagram property

$$\begin{equation} \mathop{\mathrm{div}}\nolimits (\Pi _h\boldsymbol{v})=\mathcal{P}_h(\mathop{\mathrm{div}}\nolimits \boldsymbol{v}) {\tag{4.2}} \end{equation}$$

hold, where

$$\begin{equation*} \mathcal{P}_h: \mathrm{L}^2(\Omega )\rightarrow \mathcal{U}_h:=\{v_h\in \mathrm{L}^2(\Omega ):\hspace{0.1cm}v_h|_T\in \mathcal{P}_0(T)\qquad \forall T\in \mathcal{T}_h\} \end{equation*}$$

is the standard $\mathrm{L}^2$-orthogonal projector. Then, for any $r\in (0,1]$ we have that

$$\begin{equation} \displaystyle \|q-\mathcal{P}_hq\|_{0,\Omega }\leq Ch^r\|q\|_{r,\Omega }\qquad \forall q\in \mathrm{H}^r(\Omega ). {\tag{4.3}} \end{equation}$$

Let $\mathcal{G}_h$ be the orthogonal complement of $\mathcal{K}_h$ in $\mathcal{V}_h$, and let $\widetilde{\mathcal{G}}_h:=\mathcal{G}_h\times \mathcal{G}_h\subset \widetilde{\mathcal{V}}=\mathcal{V}\times \mathcal{H}$ be endowed with the corresponding product norm. Note that $\mathcal{G}_h\nsubseteq \mathcal{G}$ and hence $\widetilde{\mathcal{G}}_h\nsubseteq \widetilde{\mathcal{G}}$.

The following result provides estimates for the terms in the Helmholtz decomposition of functions in $\mathcal{G}_h$. Let us recall that, here and thereafter, $s>0$ denotes the optimal regularity exponent such that Lemma 2.2 holds true.

Lemma 4.1.

For any $\boldsymbol{v}_h\in \mathcal{G}_h$, there exist $\xi \in \mathrm{H}^1(\Omega )$ and $\chi \in \mathcal{K}$ such that

$$\begin{equation*} \boldsymbol{v}_h=\frac{1}{\rho }\nabla \xi +\chi \end{equation*}$$

with $\frac{1}{\rho }\nabla \xi \in \mathrm{H}^{s}(\Omega _{1}\cup \Omega _2)^d$, and the following estimates hold:

$$\begin{equation*} \left\|\frac{1}{\rho }\nabla \xi \right\|_{s,\Omega _1} +\left\|\frac{1}{\rho }\nabla \xi \right\|_{s,\Omega _2}\leq C\|\mathop{\mathrm{div}}\nolimits \boldsymbol{v}_h\|_{0,\Omega }\qquad \mbox{and}\qquad \|\chi \|_{0,\Omega } \leq Ch^{s}\|\mathop{\mathrm{div}}\nolimits \boldsymbol{v}_h\|_{0,\Omega }. \end{equation*}$$

Proof.

The proof follows by repeating the arguments of the proof of Lemma 4.1 from 4, taking care of the presence of the discontinuous coefficient $\rho$.

■

Remark 4.1.

We notice that Lemma 4.1 provides (1.1) with $P \boldsymbol{v}_{h} = \frac{1}{\rho }\nabla \xi$. It is worthwhile to mention that, with the regularity results for $\xi$ at hand (see Lemma 2.2), Lemma 4.1 may also be deduced by considering the case of differential $2$-forms in Lemma 5.10 of 2.

The following result is a direct consequence of Lemma 4.1.

Lemma 4.2.

The sesquilinear form $a:\mathcal{G}_h\times \mathcal{G}_h\rightarrow \mathbb{C}$ is $\mathcal{G}_h$-elliptic, with ellipticity constant not depending on $h$. Consequently, $\widetilde{a}:\widetilde{\mathcal{G}}_h\times \widetilde{\mathcal{G}}_h\rightarrow \mathbb{C}$ is $\widetilde{\mathcal{G}}_h$-elliptic uniformly in $h$.

Now, we are in a position to introduce the discrete version of the operator $\boldsymbol{T}$. Let $\boldsymbol{T}_{h}:\widetilde{\mathcal{V}}\rightarrow \widetilde{\mathcal{V}}$ be defined by $\boldsymbol{T}_{h}(\boldsymbol{f},\boldsymbol{g}):=(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h)$ with $(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h)\in \widetilde{\mathcal{G}}_h$ the solution of

$$\begin{equation*} \widetilde{a}((\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h),(\boldsymbol{v}_h,\widehat{\boldsymbol{v}}_h))=\widetilde{b}((\boldsymbol{f},\boldsymbol{g}),(\boldsymbol{v}_h,\widehat{\boldsymbol{v}}_h))\qquad \forall (\boldsymbol{v}_h,\widehat{\boldsymbol{v}}_h)\in \widetilde{\mathcal{G}}_h. \end{equation*}$$

It is easy to check that $(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h)=\boldsymbol{T}_{h}(\boldsymbol{f},\boldsymbol{g})$ if and only if

$$\begin{equation} \widehat{\boldsymbol{u}}_h=\mathcal{P}_{\mathcal{G}_h}\boldsymbol{f}, {\tag{4.4}} \end{equation}$$

where $\mathcal{P}_{\mathcal{G}_h}$ is the $\mathcal{H}$-orthogonal projection onto $\mathcal{G}_h$, and $\boldsymbol{u}_h\in \mathcal{G}_h$ solves

$$\begin{equation} \displaystyle \int _{\Omega }\rho c^{2}\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_h\mathop{\mathrm{div}}\nolimits \boldsymbol{v}_h=-2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h-\int _{\Omega }\rho \boldsymbol{g}\cdot \overline{\boldsymbol{v}}_h\qquad \forall \boldsymbol{v}_h\in \mathcal{G}_h. {\tag{4.5}} \end{equation}$$

Since $\boldsymbol{T}_h(\widetilde{\mathcal{V}})\subset \widetilde{\mathcal{G}}_h$, $\mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T}_h)=\mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T}_h|_{\widetilde{\mathcal{G}}_h})\cup \{0\}$ holds (cf. 3, Lemma 4.1). Thus, we will restrict our attention to $\boldsymbol{T}_h|_{\widetilde{\mathcal{G}}_h}$.

As claimed above, at the beginning of Section 5, Problem 4 will be shown to be equivalent to a well-posed generalized matrix eigenvalue problem. This problem has $\lambda _h=0$ as an eigenvalue with corresponding eigenspace $\widetilde{\mathcal{K}}_h$. The rest of the eigenvalues are related with the spectrum of $\boldsymbol{T}_h|_{\widetilde{\mathcal{G}}_h}$ according to the following lemma.

Lemma 4.3.

There holds that $(\mu _h,(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h))$ is an eigenpair of $\boldsymbol{T}_{h}|_{\widetilde{\mathcal{G}}_h}$ with $\mu _h\ne 0$ if and only if $(\lambda _h,(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h))$ is a solution of Problem 4 with $\lambda _h=1/\mu _h$.

Proof.

The proof essentially follows that of Lemma 2.5, by using the fact that $\mathcal{V}_h\times \mathcal{V}_h=\widetilde{\mathcal{G}}_h\oplus (\mathcal{K}_h\times \mathcal{K}_h).$

■

Our next goal is to show that any isolated eigenvalue of $\boldsymbol{T}$ with algebraic multiplicity $m$ is approximated by exactly $m$ eigenvalues of $\boldsymbol{T}_{h}$ (repeated according to their respective algebraic multiplicities) and that spurious eigenvalues do not arise. To this end, we will adapt to our problem the theory from 4, which in turn uses arguments introduced in 1516 to deal with non-compact operators. From now on, let $\mu \in \mathop{\mathrm{Sp}_\mathrm{disc}}\nolimits (\boldsymbol{T})$, $\mu \ne 0$, be a fixed isolated eigenvalue of finite algebraic multiplicity $m$. Let $\mathcal{E}$ be the invariant subspace of $\boldsymbol{T}$ corresponding to $\mu$. Our analysis will be based on proving the following two properties:

$$\begin{align*} &\text{P1.}\quad \displaystyle \|\boldsymbol{T}-\boldsymbol{T}_{h}\|_h:=\sup _{\boldsymbol{0}\ne (\boldsymbol{f}_h,\boldsymbol{g}_h)\in \widetilde{\mathcal{G}}_h} \frac{\|(\boldsymbol{T}-\boldsymbol{T}_{h})(\boldsymbol{f}_h,\boldsymbol{g}_h)\|_{\widetilde{\mathcal{V}}}}{\Vert (\boldsymbol{f}_h,\boldsymbol{g}_h)\Vert _{\widetilde{\mathcal{V}}}}\rightarrow 0\hspace{0.3cm} \text{as} \hspace{0.15cm}h\rightarrow 0;\\ \newline &\text{P2.}\quad \forall (\boldsymbol{v},\widehat{\boldsymbol{v}})\in \mathcal{E}, \displaystyle \inf _{(\boldsymbol{v}_h,\widehat{\boldsymbol{v}}_h)\in \widetilde{\mathcal{G}}_h}\|(\boldsymbol{v},\widehat{\boldsymbol{v}}) -(\boldsymbol{v}_h,\widehat{\boldsymbol{v}}_h)\|_{\widetilde{\mathcal{V}}}\rightarrow 0 \hspace{0.3cm} \text{as} \hspace{0.15cm} h\rightarrow 0. \end{align*}$$

Let $(\boldsymbol{f}_h,\boldsymbol{g}_h)\in \widetilde{\mathcal{G}}_h$ and $(\boldsymbol{u},\widehat{\boldsymbol{u}}):=\boldsymbol{T}(\boldsymbol{f}_h,\boldsymbol{g}_h)$. From 2.17, we can write $\boldsymbol{u}=\boldsymbol{u}_1+\boldsymbol{u}_2$ with $\boldsymbol{u}_1,\boldsymbol{u}_2\in \mathcal{G}$ satisfying

$$\begin{equation} \boldsymbol{u}_1\in \mathcal{G}:\hspace{0.3cm}\displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}=-2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}_h\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}, {\tag{4.6}} \end{equation}$$

and

$$\begin{equation} \boldsymbol{u}_2\in \mathcal{G}:\hspace{0.3cm}\displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}=-\int _{\Omega }\rho \boldsymbol{g}_h\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. {\tag{4.7}} \end{equation}$$

The following result states some properties of the solutions of the problems above.

Lemma 4.4.

For $(\boldsymbol{f}_h,\boldsymbol{g}_h)\in \widetilde{\mathcal{G}}_h$, let $(\boldsymbol{u},\widehat{\boldsymbol{u}}):=\boldsymbol{T}(\boldsymbol{f}_h,\boldsymbol{g}_h)$ and consider the decomposition $\boldsymbol{u}=\boldsymbol{u}_1+\boldsymbol{u}_2$ as above. Hence, $\boldsymbol{u}_1,\boldsymbol{u}_2\in \mathrm{H}^{s}(\Omega _1\cup \Omega _2)^d$, $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1\in \mathcal{U}_h$, $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\in \mathrm{H}^{1+s}(\Omega _1\cup \Omega _2)$, and the following estimates hold true:

$$\begin{equation} \|\boldsymbol{u}_1\|_{s,\Omega _1}+\|\boldsymbol{u}_1\|_{s,\Omega _2}\leq C\|\boldsymbol{f}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }, {\tag{4.8}} \end{equation}$$

$$\begin{equation} \|\boldsymbol{u}_2\|_{s,\Omega _1}+\|\boldsymbol{u}_2\|_{s,\Omega _2}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\|_{1+s,\Omega _1}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\|_{1+s,\Omega _2} \leq C\|\boldsymbol{g}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }. {\tag{4.9}} \end{equation}$$

Proof.

Since $\boldsymbol{u}_1\in \mathcal{G}$, due to Lemma 2.2 we have that $\boldsymbol{u}_1\in \mathrm{H}^s(\Omega _1\cup \Omega _2)^d$ and $\|\boldsymbol{u}_1\|_{s,\Omega _1}+\|\boldsymbol{u}_1\|_{s,\Omega _2}\leq C\|\mathop{\mathrm{div}}\nolimits \boldsymbol{f}_h\|_{0,\Omega }$. Moreover, note that 4.6 also holds for $\boldsymbol{v}\in \mathcal{K}$ and hence for all $\boldsymbol{v}\in \mathcal{V}$. Then, we write

$$\begin{equation*} \displaystyle \int _{\Omega }(\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1+2 \nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}_h)\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}=0\qquad \forall \boldsymbol{v}\in \mathcal{V}. \end{equation*}$$

Thus, taking test functions in $\mathcal{D}(\Omega )^d\subset \mathcal{V}$ we have $\nabla (\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1+2\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}_h)=0$. Since $\rho , c, \nu$, and $\mathop{\mathrm{div}}\nolimits \boldsymbol{f}_h$ are piecewise constant, we have that $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1$ is piecewise constant as well; namely, $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1\in \mathcal{U}_h$.

On the other hand, since $\boldsymbol{u}_2\in \mathcal{G}$, by applying Lemma 2.2 again we have that $\boldsymbol{u}_2\in \mathrm{H}^{s}(\Omega _1\cup \Omega _2)^d$ and $\|\boldsymbol{u}_2\|_{s,\Omega _1}+\|\boldsymbol{u}_2\|_{s,\Omega _2}\leq C\|\boldsymbol{g}_h\|_{0,\Omega }$. To prove additional regularity for $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2$, we use Lemma 4.1 to write $\boldsymbol{g}_h=\frac{1}{\rho }\nabla \xi +\chi$ with $\chi \in \mathcal{K}$, $\frac{1}{\rho }\nabla \xi \in \mathrm{H}^{s}(\Omega _1\cup \Omega _2)^d$, and $\|\frac{1}{\rho }\nabla \xi \|_{s,\Omega _1} +\|\frac{1}{\rho }\nabla \xi \|_{s,\Omega _2}\leq C\|\mathop{\mathrm{div}}\nolimits \boldsymbol{g}_h\|_{0,\Omega }$. Moreover, since $\rho$ is constant in each subdomain $\Omega _i$, also $\nabla \xi |_{\Omega _i}\in \mathrm{H}^{s}(\Omega _i)^d$, $i=1,2$. Then, from 4.7 we have that

$$\begin{equation*} \displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= -\int _{\Omega }\nabla \xi \cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. \end{equation*}$$

Since the above equation trivially holds for $\boldsymbol{v}\in \mathcal{K}$ too, it holds for all $\boldsymbol{v}\in \mathcal{V}$. Then, by testing it with $\boldsymbol{v}\in \mathcal{D}(\Omega )^d$ we have that $\nabla (\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2)=-\nabla \xi \in \Omega$. Therefore, by restricting to $\Omega _i$, $i=1,2$, we have that $\nabla (\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2|_{\Omega _i})=-\nabla \xi |_{\Omega _i}\in \mathrm{H}^{s}(\Omega _i)^d$. Since $\rho$ and $c$ are piecewise constant, we conclude that $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2|_{\Omega _i}\in \mathrm{H}^{1+s}(\Omega _1\cup \Omega _2)$, and

$$\begin{equation*} \|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\|_{1+s,\Omega _1}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\|_{1+s,\Omega _2}\leq C\|\nabla \xi \|_{0,\Omega }\leq C\|\mathop{\mathrm{div}}\nolimits \boldsymbol{g}_h\|_{0,\Omega }. \end{equation*}$$

Hence, we conclude the proof.

■

We consider a similar decomposition in the discrete case. For $(\boldsymbol{f}_h,\boldsymbol{g}_h)\in \widetilde{\mathcal{G}}_h$, let $(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h):=\boldsymbol{T}_h(\boldsymbol{f}_h,\boldsymbol{g}_h)$. We write $\boldsymbol{u}_h=\boldsymbol{u}_{1h}+\boldsymbol{u}_{2h}$ with $\boldsymbol{u}_{1h}$ and $\boldsymbol{u}_{2h}$ satisfying

$$\begin{equation} \boldsymbol{u}_{1h}\in \mathcal{G}_h:\hspace{0.3cm}\displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{1h}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h=-2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{f}_h\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h\qquad \forall \boldsymbol{v}_h\in \mathcal{G}_h, {\tag{4.10}} \end{equation}$$

and

$$\begin{equation} \boldsymbol{u}_{2h}\in \mathcal{G}_h:\hspace{0.3cm}\displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{2h}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h=-\int _{\Omega }\rho \boldsymbol{g}_h\cdot \overline{\boldsymbol{v}}_h\qquad \forall \boldsymbol{v}_h\in \mathcal{G}_h. {\tag{4.11}} \end{equation}$$

These are the finite element discretizations of problems 4.6 and 4.7, respectively, and the following error estimates hold true.

Lemma 4.5.

Let $(\boldsymbol{f}_h,\boldsymbol{g}_h)\in \widetilde{\mathcal{G}}_h$. Let $\boldsymbol{u}_1,\boldsymbol{u}_2$ be the solutions of problems 4.6 and 4.7, respectively, and $\boldsymbol{u}_{1h},\boldsymbol{u}_{2h}$ those of problems 4.10 and 4.11, respectively. Then, the following estimates hold true:

$$\begin{equation} \|\boldsymbol{u}_{1}-\boldsymbol{u}_{1h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq Ch^s\|\boldsymbol{f}_h\|_{\mathcal{V}}, {\tag{4.12}} \end{equation}$$

$$\begin{equation} \|\boldsymbol{u}_{2}-\boldsymbol{u}_{2h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq Ch^s\|\boldsymbol{g}_h\|_{\mathcal{H}}. {\tag{4.13}} \end{equation}$$

Proof.

Since $\mathcal{G}_h\nsubseteq \mathcal{G}$, we will resort to the second Strang Lemma, which for problems 4.6 and 4.10 reads as follows:

$$\begin{equation} \displaystyle \|\boldsymbol{u}_1-\boldsymbol{u}_{1h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq C\left[\inf _{\boldsymbol{v}_{h}\in \mathcal{G}_h}\|\boldsymbol{u}_1-\boldsymbol{v}_{h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }+\sup _{\boldsymbol{0}\neq \boldsymbol{v}_h\in \mathcal{G}_h}\frac{a(\boldsymbol{u}_1-\boldsymbol{u}_{1h},\boldsymbol{v}_h)}{\|\boldsymbol{v}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }}\right]. {\tag{4.14}} \end{equation}$$

Because of Lemma 4.4, $\Pi _h\boldsymbol{u}_1$ is well-defined. Since $\Pi _{h}\boldsymbol{u}_1\in \mathcal{V}_h=\mathcal{G}_h\oplus \mathcal{K}_h$, there exist $\widetilde{\boldsymbol{u}}_{1h}\in \mathcal{G}_h$ and $\breve{\boldsymbol{u}}_h\in \mathcal{K}_h$ such that $\Pi _{h}\boldsymbol{u}_1=\widetilde{\boldsymbol{u}}_{1h}+\breve{\boldsymbol{u}}_h$. Then, since $\boldsymbol{u}_1-\widetilde{\boldsymbol{u}}_{1h}$ is orthogonal to $\breve{\boldsymbol{u}}_h$, we observe that

$$\begin{align} \|\boldsymbol{u}_1-\widetilde{\boldsymbol{u}}_{1h}\|_{\mathcal{V}}^2 &\leq \|\boldsymbol{u}_1-\widetilde{\boldsymbol{u}}_{1h}\|_{\mathcal{V}}^2 +\|\breve{\boldsymbol{u}}_h\|_{\mathcal{V}}^2\\ &=\|(\widetilde{\boldsymbol{u}}_{1h}-\boldsymbol{u}_1)+\breve{\boldsymbol{u}}_h\|_{\mathcal{V}}^2 =\|\boldsymbol{u}_1-\Pi _h\boldsymbol{u}_1\|_{\mathcal{V}}^2\\ &\leq C\left(\|\boldsymbol{u}_1-\Pi _h\boldsymbol{u}_1\|_{0,\Omega }^2 +\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1-\mathop{\mathrm{div}}\nolimits (\Pi _h\boldsymbol{u}_1)\|_{0,\Omega }^2\right). \end{align}$$

The first term on the right-hand side above is bounded as follows:

$$\begin{equation*} \|\boldsymbol{u}_1-\Pi _h\boldsymbol{u}_1\|_{0,\Omega }\leq Ch^{s}(\|\boldsymbol{u}_1\|_{s,\Omega _1}+\|\boldsymbol{u}_1\|_{s,\Omega _2}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1\|_{0,\Omega })\leq Ch^{s}\|\boldsymbol{f}_h\|_{\mathcal{V}}, \end{equation*}$$

where we have used 4.1, 4.8, and the fact that $\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1\|_{0,\Omega }\leq C\|\mathop{\mathrm{div}}\nolimits \boldsymbol{f}_h\|_{0,\Omega }$, which in turn follows from 4.6 by taking $\boldsymbol{v}=\boldsymbol{f}_h$. On the other hand, the second term vanishes because of 4.2 since $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_1\in \mathcal{U}_h$ (cf. Lemma 4.4). Hence, $\|\boldsymbol{u}_1-\widetilde{\boldsymbol{u}}_{1h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq Ch^{s} \|\boldsymbol{f}_h\|_{\mathcal{V}}$, which allows us to control the approximation term in 4.14.

For the consistency term, it is enough to recall that 4.6 holds for all $\boldsymbol{v}\in \mathcal{V}$. Then, by using 4.10, it is easy to check that $a(\boldsymbol{u}_1-\boldsymbol{u}_{1h},\boldsymbol{v}_h)=0$ for all $\boldsymbol{v}_h\in \mathcal{G}_h\subset \mathcal{V}$. From this, the Strang estimate for $\|\boldsymbol{u}_1-\boldsymbol{u}_{1h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }$ reads as follows:

$$\begin{equation*} \displaystyle \|\boldsymbol{u}_1-\boldsymbol{u}_{1h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq C\inf _{\boldsymbol{v}_{h}\in \mathcal{G}_h}\|\boldsymbol{u}_1-\boldsymbol{v}_{h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq Ch^{s}\|\boldsymbol{f}_h\|_{\mathcal{V}}. \end{equation*}$$

Thus 4.12 holds true.

To prove 4.13, we resort again to the second Strang Lemma:

$$\begin{equation} \displaystyle \|\boldsymbol{u}_2-\boldsymbol{u}_{2h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq C\left[\inf _{\boldsymbol{v}_{h}\in \mathcal{G}_h}\|\boldsymbol{u}_2-\boldsymbol{v}_{h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega } +\sup _{\boldsymbol{0}\neq \boldsymbol{v}_h\in \mathcal{G}_h}\frac{a(\boldsymbol{u}_2-\boldsymbol{u}_{2h},\boldsymbol{v}_h)}{\|\boldsymbol{v}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }}\right]. {\tag{4.15}} \end{equation}$$

Since $\boldsymbol{u}_2\in \mathrm{H}^{s}(\Omega _1\cup \Omega _2)^d$ (cf. Lemma 4.4), we have that $\Pi _h\boldsymbol{u}_2$ is well-defined. We proceed as above and write $\Pi _{h}\boldsymbol{u}_2=\widetilde{\boldsymbol{u}}_{2h}+\check{\boldsymbol{u}}_h$ with $\widetilde{\boldsymbol{u}}_{2h}\in \mathcal{G}_h$ and $\check{\boldsymbol{u}}_h\in \mathcal{K}_h$ to obtain

$$\begin{equation} \|\boldsymbol{u}_2-\widetilde{\boldsymbol{u}}_{2h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq C\left[\|\boldsymbol{u}_2-\Pi _h\boldsymbol{u}_2\|_{0,\Omega }+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2-\mathop{\mathrm{div}}\nolimits (\Pi _h\boldsymbol{u}_2)\|_{0,\Omega }\right]. \tag{4.16} \end{equation}$$

For the first term on the right-hand side above, 4.1 and Lemma 4.4 yield

$$\begin{equation*} \|\boldsymbol{u}_2-\Pi _h\boldsymbol{u}_2\|_{0,\Omega }\leq Ch^{s}(\|\boldsymbol{u}_2\|_{s,\Omega _1}+\|\boldsymbol{u}_2\|_{s,\Omega _2} +\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\|_{0,\Omega })\leq Ch^{s}\|\boldsymbol{g}_h\|_{\mathcal{H}}. \end{equation*}$$

For the second term, we have from 4.3 and from Lemma 4.4 again

$$\begin{align*} \displaystyle \|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2-\mathop{\mathrm{div}}\nolimits \Pi _h\boldsymbol{u}_2\|_{0,\Omega }^2&=\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2-\mathcal{P}_h(\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2)\|_{0,\Omega }^2,\\ &\leq Ch(\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\|_{1,\Omega _1}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\|_{1,\Omega _2})\leq Ch\|\boldsymbol{g}_h\|_{\mathcal{H}}. \end{align*}$$

Hence, $\|\boldsymbol{u}_2-\widetilde{\boldsymbol{u}}_{2h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq Ch^{s}\|\boldsymbol{g}_h\|_{\mathcal{H}}$, which allows us to bound the approximation term in 4.15.

For the consistency term, given $\boldsymbol{v}_h\in \mathcal{G}_h$ we apply Lemma 4.1 to write $\boldsymbol{v}_h=\frac{1}{\rho }\nabla \xi +\chi$ with $\frac{1}{\rho }\nabla \xi \in \mathrm{H}^{s}(\Omega _1\cup \Omega _2)^d$, $\chi \in \mathcal{K}$, and $\|\chi \|_{0,\Omega }\leq Ch^{s}\|\mathop{\mathrm{div}}\nolimits \boldsymbol{v}_h\|_{0,\Omega }$. Then, from 4.7 we have

$$\begin{equation*} \displaystyle a(\boldsymbol{u}_2,\boldsymbol{v}_h)=\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h= \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_2\mathop{\mathrm{div}}\nolimits \left(\frac{1}{\rho }\nabla \overline{\xi }\right) =\int _{\Omega }\boldsymbol{g}_h\cdot \nabla \overline{\xi }. \end{equation*}$$

On the other hand, from 4.11,

$$\begin{equation*} \displaystyle a(\boldsymbol{u}_{2h},\boldsymbol{v}_h)=\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{2h}\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}_h=\int _{\Omega }\rho \boldsymbol{g}_h\cdot \overline{\boldsymbol{v}}_h =\int _{\Omega }\boldsymbol{g}_h\cdot \nabla \overline{\xi }+\int _{\Omega }\rho \boldsymbol{g}_h\cdot \overline{\chi }. \end{equation*}$$

Therefore,

$$\begin{equation*} a(\boldsymbol{u}_2-\boldsymbol{u}_{2h},\boldsymbol{v}_h)=-\int _{\Omega }\rho \boldsymbol{g}_h\cdot \overline{\chi } \leq Ch^{s}\|\boldsymbol{g}_h\|_{0,\Omega }\|\boldsymbol{v}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega } \end{equation*}$$

and, hence,

$$\begin{equation*} \displaystyle \sup _{\boldsymbol{0}\ne \boldsymbol{v}_h\in \mathcal{G}_h}\frac{a(\boldsymbol{u}_2-\boldsymbol{u}_{2h},\boldsymbol{v}_h)}{\|\boldsymbol{v}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }}\leq Ch^{s}\|\boldsymbol{g}_h\|_{0,\Omega }, \end{equation*}$$

which allows us to complete the proof.

■

Now, we are in a position to establish the following result.

Lemma 4.6.

Property P1 holds true. Moreover,

$$\begin{equation*} \|\boldsymbol{T}-\boldsymbol{T}_{h}\|_h\leq Ch^{s}. \end{equation*}$$

Proof.

For $(\boldsymbol{f}_h,\boldsymbol{g}_h)\in \widetilde{\mathcal{G}}_h$, let $(\boldsymbol{u},\widehat{\boldsymbol{u}}):=\boldsymbol{T}(\boldsymbol{f}_h,\boldsymbol{g}_h)$ and $(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h):=\boldsymbol{T}_{h}(\boldsymbol{f}_h,\boldsymbol{g}_h)$. From 2.16 and 4.4 we have that $\widehat{\boldsymbol{u}}-\widehat{\boldsymbol{u}}_h=\boldsymbol{f}_h-\mathcal{P}_{\mathcal{G}_h}\boldsymbol{f}_h=0$. Hence, by writing $\boldsymbol{u}=\boldsymbol{u}_{1}+\boldsymbol{u}_{2}$ and $\boldsymbol{u}_h=\boldsymbol{u}_{1h}+\boldsymbol{u}_{2h}$ as in Lemma 4.5, we have from this lemma

$$\begin{equation*} \displaystyle \|\boldsymbol{T}-\boldsymbol{T}_{h}\|_h \leq \sup _{\boldsymbol{0}\ne (\boldsymbol{g}_h,\boldsymbol{f}_h)\in \widetilde{\mathcal{G}}_h} \frac{\|\boldsymbol{u}_1-\boldsymbol{u}_{1h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }+\|\boldsymbol{u}_2-\boldsymbol{u}_{2h}\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }}{\|(\boldsymbol{f}_h,\boldsymbol{g}_h)\|_{\widetilde{\mathcal{V}}}} \leq Ch^{s}. \end{equation*}$$

Thus, we conclude the proof.

■

Our next goal is to prove property P2. With this aim, first we will prove the following additional regularity result.

Lemma 4.7.

Let $(\boldsymbol{u},\widehat{\boldsymbol{u}})\in \mathcal{E}$. Then, $\boldsymbol{u},\widehat{\boldsymbol{u}}\in \mathcal{G}\subset \mathrm{H}^s(\Omega _1\cup \Omega _2)^d$, $\mathop{\mathrm{div}}\nolimits \boldsymbol{u},\mathop{\mathrm{div}}\nolimits \widehat{\boldsymbol{u}}\in \mathrm{H}^{1+s}(\Omega _1\cup \Omega _2)$, and

$$\begin{align} \|\boldsymbol{u}\|_{s,\Omega _1}+\|\boldsymbol{u}\|_{s,\Omega _2}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\|_{1+s,\Omega _1}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\|_{1+s,\Omega _2}\leq & C \|(\boldsymbol{u},\widehat{\boldsymbol{u}})\|_{\widetilde{\mathcal{V}}},{\tag{4.17}}\\ \|\widehat{\boldsymbol{u}}\|_{s,\Omega _1}+\|\widehat{\boldsymbol{u}}\|_{s,\Omega _2}+\|\mathop{\mathrm{div}}\nolimits \widehat{\boldsymbol{u}}\|_{1+s,\Omega _1}+\|\mathop{\mathrm{div}}\nolimits \widehat{\boldsymbol{u}}\|_{1+s,\Omega _2}\leq & C \|(\boldsymbol{u},\widehat{\boldsymbol{u}})\|_{\widetilde{\mathcal{V}}}. {\tag{4.18}} \end{align}$$

Proof.

We prove the above inequalities for all the generalized eigenfunctions of $\boldsymbol{T}$. Let $\{(\boldsymbol{u}_k,\widehat{\boldsymbol{u}}_k)\}_{k=1}^p$ be a Jordan chain of the operator $\boldsymbol{T}$ associated with $\mu$. Then, $\boldsymbol{T}(\boldsymbol{u}_k,\widehat{\boldsymbol{u}}_k)=\mu (\boldsymbol{u}_k,\widehat{\boldsymbol{u}}_k)+(\boldsymbol{u}_{k-1}, \widehat{\boldsymbol{u}}_{k-1})$, $k=1,\ldots ,p$, with $(\boldsymbol{u}_0,\widehat{\boldsymbol{u}}_0)=\boldsymbol{0}$. We will use an induction argument on $k$. Assume that $\boldsymbol{u}_{k-1}$ and $\widehat{\boldsymbol{u}}_{k-1}$ belong to $\mathcal{G}$ and satisfy 4.17 and 4.18, respectively (which obviously hold for $k=1$). First note that, because of the boundedness of $\boldsymbol{T}$, we have

$$\begin{equation} \|(\boldsymbol{u}_{k-1},\widehat{\boldsymbol{u}}_{k-1})\|_{\widetilde{\mathcal{V}}}\leq C\|(\boldsymbol{u}_k,\widehat{\boldsymbol{u}}_k)\|_{\widetilde{\mathcal{V}}}. {\tag{4.19}} \end{equation}$$

On the other hand, by using 2.16 and 2.17 we have that

$$\begin{equation} \mu \widehat{\boldsymbol{u}}_k+\widehat{\boldsymbol{u}}_{k-1}=\boldsymbol{u}_k\hspace{0.2cm}\mbox{ in}\ \Omega {\tag{4.20}} \end{equation}$$

and that $\mu \boldsymbol{u}_k+\boldsymbol{u}_{k-1}\in \mathcal{G}$ satisfies

$$\begin{equation*} \displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits (\mu \boldsymbol{u}_k+\boldsymbol{u}_{k-1})\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= -2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{u}_k\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}-\int _{\Omega }\rho \widehat{\boldsymbol{u}}_k\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{G}. \end{equation*}$$

Hence, $\boldsymbol{u}_k,\widehat{\boldsymbol{u}}_k\in \mathcal{G}$.

We observe that the equation above also holds for any $\boldsymbol{v}\in \mathcal{K}$. Then,

$$\begin{equation*} \displaystyle \int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits (\mu \boldsymbol{u}_k+\boldsymbol{u}_{k-1})\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}= -2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \boldsymbol{u}_k\mathop{\mathrm{div}}\nolimits \overline{\boldsymbol{v}}-\int _{\Omega }\rho \widehat{\boldsymbol{u}}_k\cdot \overline{\boldsymbol{v}}\qquad \forall \boldsymbol{v}\in \mathcal{V}. \end{equation*}$$

Thus, considering test functions in $\mathcal{D}(\Omega )^d\subset \mathcal{V}$ we obtain

$$\begin{equation} \displaystyle \nabla ((\mu \rho c^2+2\nu )\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_k)=\rho \widehat{\boldsymbol{u}}_k-\nabla (\rho c^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{k-1}). {\tag{4.21}} \end{equation}$$

Let us assume that $\mu \rho c^2+2\nu \ne 0$ in both $\Omega _1$ and $\Omega _2$ (we discuss the other case at the end of the proof). Hence, since $\rho$, $c$, and $\nu$ are constant in each $\Omega _i$, $\rho _i\widehat{\boldsymbol{u}}_k-\nabla (\rho _i c_i^2\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{k-1})\in \mathrm{L}^2(\Omega _i)^d$, $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_k|_{\Omega _i}\in \mathrm{H}^1(\Omega _i)$, and

$$\begin{equation*} \|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_k\|_{1,\Omega _i}\leq C\left(\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{k-1}\|_{1,\Omega _i} +\|(\boldsymbol{u}_k,\widehat{\boldsymbol{u}}_k)\|_{\widetilde{\mathcal{V}}}\right),\qquad i=1,2. \end{equation*}$$

Now, since $\boldsymbol{u}_k\in \mathcal{G}$, due to Lemma 2.2 we have that $\boldsymbol{u}_k\in \mathrm{H}^s(\Omega _1\cup \Omega _2)^d$. Then, from 2.13 and the previous estimate we have

$$\begin{equation} \|\boldsymbol{u}_k\|_{s,\Omega _1}+\|\boldsymbol{u}_k\|_{s,\Omega _2}\leq C\left(\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{k-1}\|_{1,\Omega _1} +\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{k-1}\|_{1,\Omega _2}+\|(\boldsymbol{u}_k,\widehat{\boldsymbol{u}}_k)\|_{\widetilde{\mathcal{V}}}\right). {\tag{4.22}} \end{equation}$$

On the other hand, from 4.20 we obtain

$$\begin{equation} \displaystyle \|\widehat{\boldsymbol{u}}_k\|_{s,\Omega _i}\leq \frac{1}{\mu } \left(\|\widehat{\boldsymbol{u}}_{k-1}\|_{s,\Omega _i}+\|\boldsymbol{u}_k\|_{s,\Omega _i}\right),\qquad i=1,2, {\tag{4.23}} \end{equation}$$

and, from 4.21,

$$\begin{equation} \|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_k\|_{1+s,\Omega _i}\leq C(\|\widehat{\boldsymbol{u}}_k\|_{s,\Omega _i}+\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_{k-1}\|_{1+s,\Omega _i}),\qquad i=1,2. {\tag{4.24}} \end{equation}$$

Finally, from 4.20 again,

$$\begin{equation} \displaystyle \|\mathop{\mathrm{div}}\nolimits \widehat{\boldsymbol{u}}_k\|_{1+s,\Omega _i}\leq \frac{1}{\mu }(\|\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_k\|_{1+s,\Omega _i} +\|\mathop{\mathrm{div}}\nolimits \widehat{\boldsymbol{u}}_{k-1}\|_{1+s,\Omega _i}),\qquad i=1,2. {\tag{4.25}} \end{equation}$$

Hence, from inequalities 4.22–4.25, the inductive assumption, and 4.19, we derive 4.17 and 4.18 provided $\mu \rho c^2+2\nu \ne 0$ in both $\Omega _1$ and $\Omega _2$.

In the case that $\mu \rho c^2+2\nu$ vanishes in $\Omega _i$, $i=1$ or $2$, arguing as in Remark 2.3 we obtain that $\boldsymbol{u}_1|_{\Omega _i}=\widehat{\boldsymbol{u}}_1|_{\Omega _i}=\boldsymbol{0}$ and, once again, an induction argument allows us to conclude that $\boldsymbol{u}_k, \widehat{\boldsymbol{u}}_k=\boldsymbol{0}$ in $\Omega _i$, $k=1,\ldots ,p$. The proof is complete.

■

Now, we are in a position to establish property P2.

Lemma 4.8.

Property P2 holds true. Moreover, for any $(\boldsymbol{u},\widehat{\boldsymbol{u}})\in \mathcal{E}$, there exist $\widetilde{\boldsymbol{u}}_h$, $\widetilde{\widehat{\boldsymbol{u}}}_h\in \mathcal{G}_h$ such that

$$\begin{equation*} \|\boldsymbol{u}-\widetilde{\boldsymbol{u}}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }+\|\widehat{\boldsymbol{u}}-\widetilde{\widehat{\boldsymbol{u}}}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq Ch^s\|(\boldsymbol{u},\widehat{\boldsymbol{u}})\|_{\widetilde{\mathcal{V}}}. \end{equation*}$$

Proof.

Let $(\boldsymbol{u},\widehat{\boldsymbol{u}})\in \mathcal{E}$. According to Lemma 4.7, we have that $\boldsymbol{u},\widehat{\boldsymbol{u}}\in \mathrm{H}^s(\Omega _1\cup \Omega _2)^d$ and $\mathop{\mathrm{div}}\nolimits \boldsymbol{u},\mathop{\mathrm{div}}\nolimits \widehat{\boldsymbol{u}}\in \mathrm{H}^{1+s}(\Omega _1\cup \Omega _2)$. Let $\Pi _h\boldsymbol{u}\in \mathcal{V}_h$ be the Raviart-Thomas interpolant of $\boldsymbol{u}$. Since $\mathcal{V}_h=\mathcal{G}_h\oplus \mathcal{K}_h$, we decompose $\Pi _h\boldsymbol{u}=\widetilde{\boldsymbol{u}}_h+\breve{\boldsymbol{u}}_h$ with $\widetilde{\boldsymbol{u}}_h\in \mathcal{G}_h$ and $\breve{\boldsymbol{u}}_h\in \mathcal{K}_h$. The same arguments from the proof of Lemma 4.5 that lead to 4.13 apply in this case and combined with Lemma 4.7 allow us to prove that $\|\boldsymbol{u}-\widetilde{\boldsymbol{u}}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega } \leq Ch^s\|(\boldsymbol{u},\widehat{\boldsymbol{u}})\|_{\widetilde{\mathcal{V}}}$. A similar procedure can be used to define $\widetilde{\widehat{\boldsymbol{u}}}_h$ and to prove that $\|\widehat{\boldsymbol{u}}-\widetilde{\widehat{\boldsymbol{u}}}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }\leq Ch^s\|(\boldsymbol{u},\widehat{\boldsymbol{u}})\|_{\widetilde{\mathcal{V}}}$.

■

We also have the following auxiliary result when the source terms are in $\mathcal{E}$, whose proof follows the arguments in Case 2 of Lemma 4.7 of 4.

Lemma 4.9.

For $(\boldsymbol{f},\boldsymbol{g})\in \mathcal{E}$, let $(\boldsymbol{u},\widehat{\boldsymbol{u}}):=\boldsymbol{T}(\boldsymbol{f},\boldsymbol{g})$ and $(\boldsymbol{u}_h,\widehat{\boldsymbol{u}}_h):=\boldsymbol{T}_{h}(\boldsymbol{f},\boldsymbol{g})$. Then,

$$\begin{equation*} \|\boldsymbol{u}-\boldsymbol{u}_h\|_{\mathop{\mathrm{div}}\nolimits ,\Omega }+\|\widehat{\boldsymbol{u}}-\widehat{\boldsymbol{u}}_h\|_{0,\Omega }\leq Ch^s \|(\boldsymbol{f},\boldsymbol{g})\|_{\widetilde{\mathcal{V}}}. \end{equation*}$$

The above lemmas are the ingredients to prove spectral convergence and to obtain error estimates. Our first result is the following theorem which has been proved in 15 as a consequence of property P1 (cf. Lemma 4.6) and which shows that the proposed method is free of spurious modes.

Theorem 4.1.

Let $K\subset \mathbb{C}$ be a compact set such that $K\cap \mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T})=\emptyset$. Then, there exists $h_0>0$ such that, for all $h\leq h_0$, $K\cap \mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T}_{h})=\emptyset .$

Let $D\subset \mathbb{C}$ be a closed disk centered at $\mu$, such that $D\cap \mathop{\mathrm{Sp}}\nolimits (\boldsymbol{T})=\{\mu \}$. Let $\mu _{1h},\ldots ,\mu _{m(h)h}$ be the eigenvalues of $\boldsymbol{T}_h$ contained in $D$ (repeated according to their algebraic multiplicities). Under assumptions P1 and P2, it is proved in 15 that $m(h)=m$ for $h$ small enough and that $\lim _{h\rightarrow 0}\mu _{kh}=\mu$ for $k=1,\ldots , m$.

On the other hand the arguments used in Section 5 of 4 can be readily adapted to our problem, to obtain error estimates. We recall the definition of the gap between two closed subspaces $\mathcal{W}$ and $\mathcal{Y}$ of $\widetilde{\mathcal{V}}$:

$$\begin{equation*} \widehat{\delta }(\mathcal{W},\mathcal{Y}):=\max \{\delta (\mathcal{W},\mathcal{Y}), \delta (\mathcal{Y},\mathcal{W})\}, \end{equation*}$$

with

$$\begin{equation*} \displaystyle \delta (\mathcal{W},\mathcal{Y}):=\sup _{\stackrel{\scriptstyle (\boldsymbol{\phi },\boldsymbol{\psi })\in \mathcal{W}}{\|(\boldsymbol{\phi },\boldsymbol{\psi })\|_{\widetilde{\mathcal{V}}}=1}}\left[\inf _{(\widehat{\boldsymbol{\phi }},\widehat{\boldsymbol{\psi }})\in \mathcal{Y}}\|(\boldsymbol{\phi },\boldsymbol{\psi })-(\widehat{\boldsymbol{\phi }},\widehat{\boldsymbol{\psi }})\|_{\widetilde{\mathcal{V}}}\right]. \end{equation*}$$

Let $\mathcal{E}_h$ be the invariant subspace of $\boldsymbol{T}_h$ relative to the eigenvalues $\mu _{1h},\ldots ,\mu _{mh}$ converging to $\mu$. From Lemmas 4.6–4.9, we derive the following results for which we do not include proofs to avoid repeating step by step those of 4, Section 5.

Theorem 4.2.

There exist constants $h_0>0$ and $C>0$ such that, for all $h\leq h_0$,

$$\begin{equation*} \widehat{\delta }\left(\mathcal{E}_h,\mathcal{E}\right)\leq Ch^s. \end{equation*}$$

Theorem 4.3.

There exist constants $h_0>0$ and $C>0$ such that, for all $h\leq h_0$,

$$\begin{equation*} \displaystyle \left|\mu -\frac{1}{m}\sum _{k=1}^m\mu _{kh}\right|\leq Ch^{2s}, \end{equation*}$$

$$\begin{equation*} \displaystyle \left|\frac{1}{\mu }-\frac{1}{m}\sum _{k=1}^m\frac{1}{\mu _{kh}}\right|\leq Ch^{2s}, \end{equation*}$$

$$\begin{equation*} \displaystyle \max _{k=1,\ldots , m}|\mu -\mu _{kh}|\leq Ch^{2s/p}, \end{equation*}$$

where $p$ is the ascent of the eigenvalue $\mu$ of $\boldsymbol{T}$.

5. Numerical results

We implemented the proposed method in a MATLAB code. We report in this section the results of some numerical tests, in order to assess its performance. To this end, first we introduce a convenient matrix form of the discrete problem which allows us to use standard eigensolvers. As a byproduct, this matrix form also allows us to prove that Problems 3 and 4 are well-posed.

Let $\{\phi _j\}_{j=1}^N$ be a nodal basis of $\mathcal{V}_h$. We define the matrices $\mathbf{K}_1:=(\mathbf{K}^{(1)}_{ij})$, $\mathbf{K}_2:=(\mathbf{K}^{(2)}_{ij})$, and $\mathbf{M}:=(\mathbf{M}_{ij})$ as follows:

$$\begin{equation*} \displaystyle \mathbf{K}^{(1)}_{ij}:=2\int _{\Omega }\nu \mathop{\mathrm{div}}\nolimits \phi _i\mathop{\mathrm{div}}\nolimits \phi _j,\ \ \mathbf{K}^{(2)}_{ij}:=\int _{\Omega }\rho c^2\mathop{\mathrm{div}}\nolimits \phi _i\mathop{\mathrm{div}}\nolimits \phi _j,\ \ and\ \ \mathbf{M}_{ij}:=\int _{\Omega }\rho \phi _i\cdot \phi _j. \end{equation*}$$

The matrix form of Problem 3 reads

$$\begin{equation} (\lambda ^2_h\mathbf{M}+\lambda _h\mathbf{K}_1+\mathbf{K}_2)\vec{\boldsymbol{u}}_h=\boldsymbol{0}, {\tag{5.1}} \end{equation}$$

where we denote by $\vec{\boldsymbol{u}}_h$ the vector of components of $\boldsymbol{u}_h$ in the nodal basis of $\mathcal{V}_h$.

Analogously, the matrix form of Problem 4 reads

$$\begin{equation*} \begin{pmatrix} \mathbf{K}_2&\boldsymbol{0}\\ \boldsymbol{0}&\mathbf{M}\end{pmatrix} \begin{pmatrix} \vec{\boldsymbol{u}}_h\\ \vec{\widehat{\boldsymbol{u}}}_h \end{pmatrix} =\lambda _h\begin{pmatrix} -\mathbf{K}_1&-\mathbf{M}\\ \mathbf{M}&\boldsymbol{0}\end{pmatrix} \begin{pmatrix} \vec{\boldsymbol{u}}_h\\ \vec{\widehat{\boldsymbol{u}}}_h \end{pmatrix}, \end{equation*}$$

with $\vec{\widehat{\boldsymbol{u}}}_h$ the vector of components of $\widehat{\boldsymbol{u}}_h$. However, this problem is not suitable to be solved with standard eigensolvers, since neither the right-hand side nor the left-hand side matrix are Hermitian and positive definite.

Alternatively, for $\lambda _h\neq 0$, let $\mu _h:=\frac{1}{\lambda _h}$. Then, problem 5.1 is equivalent to

$$\begin{equation*} (\mathbf{M}+2\mu _h\mathbf{K}_1+\mu _h^2\mathbf{K}_2)\vec{\boldsymbol{u}}_h=\boldsymbol{0}. \end{equation*}$$

Introducing $\vec{\boldsymbol{w}}_h:=\mu _h\vec{\boldsymbol{u}}_h$, the problem above is equivalent to

$$\begin{equation*} \begin{pmatrix} \mathbf{M}&\boldsymbol{0}\\ \boldsymbol{0}&\mathbf{M}\end{pmatrix} \begin{pmatrix} \vec{\boldsymbol{u}}_h\\ \vec{\boldsymbol{w}}_h \end{pmatrix} =\mu _h \begin{pmatrix} -\mathbf{K}_1&-\mathbf{K}_2\\ \mathbf{M}&\boldsymbol{0}\end{pmatrix} \begin{pmatrix} \vec{\boldsymbol{u}}_h\\ \vec{\boldsymbol{w}}_h \end{pmatrix}, \end{equation*}$$ which in turn is equivalent to

$$\begin{equation*} \begin{pmatrix} -\mathbf{K}_1&-\mathbf{K}_2\\ \mathbf{M}&\boldsymbol{0}\end{pmatrix} \begin{pmatrix} \vec{\boldsymbol{u}}_h\\ \vec{\boldsymbol{w}}_h \end{pmatrix} =\lambda _h \begin{pmatrix} \mathbf{M}&\boldsymbol{0}\\ \boldsymbol{0}&\mathbf{M}\end{pmatrix} \begin{pmatrix} \vec{\boldsymbol{u}}_h\\ \vec{\boldsymbol{w}}_h \end{pmatrix}. \end{equation*}$$

Thus, the last problem is equivalent to Problem 3 except for $\lambda _h=0$ and the matrix in its right-hand side is Hermitian and positive definite. Hence, it is well-posed and can be safely solved by standard eigensolvers.

Test 1.

We applied the proposed method to a 2D rectangular rigid cavity filled with two fluids with different physical parameters as shown in Figure 2. The domain occupied by the fluids are $\Omega _1:=(0,A)\times (0,H)$ and $\Omega _2:=(0,A)\times (H,B)$. For such a simple geometry, it is possible to calculate an analytical solution which will be used to validate our method.

Let $\boldsymbol{u}\in \mathrm{H}_0(\mathop{\mathrm{div}}\nolimits ,\Omega )$ be a solution of Problem 1. Testing it with $\boldsymbol{v}\in \mathcal{D}(\Omega )^2$ we have $\nabla ((2\lambda \nu +\rho c^2)\mathop{\mathrm{div}}\nolimits \boldsymbol{u})=-\lambda ^2\rho \boldsymbol{u}\in \mathrm{L}^2(\Omega )^2$. Then, $\widehat{p}:=-(2\nu \lambda +\rho c^2)\mathop{\mathrm{div}}\nolimits \boldsymbol{u}\in \mathrm{H}^1(\Omega )$. Hence, $\widehat{p}_1|_{\Gamma }=\widehat{p}_2|_{\Gamma }$. Moreover, $\boldsymbol{u}=-\frac{1}{\lambda ^2\rho }\nabla \widehat{p}$, which implies that $\frac{1}{\rho _1}\frac{\partial \widehat{p}_1}{\partial \nu } =\frac{1}{\rho _2}\frac{\partial \widehat{p}_2}{\partial \nu }$ on $\Gamma$. Then, we write problem 2.2–2.9 in terms of $\widehat{p}_i$ as follows:

$$\begin{align*} \Delta \widehat{p}_i&=\frac{\lambda ^2\rho _i}{\rho _i c_i^2+2\nu _i\lambda }\widehat{p}_i\hspace{0.6cm}\mbox{in}\ \Omega _i,\hspace{0.2cm}i=1,2,\\ \frac{\partial \widehat{p}_i}{\partial \boldsymbol{n}_i}&=0\hspace{2.5cm}\mbox{on}\ \Gamma _i,\hspace{0.2cm}i=1,2,\\ \widehat{p}_1&=\widehat{p}_2\hspace{2.35cm}\mbox{on}\ \Gamma ,\\ \frac{1}{\rho _1}\frac{\partial \widehat{p}_1}{\partial \boldsymbol{n}}&=\frac{1}{\rho _2}\frac{\partial \widehat{p}_2}{\partial \boldsymbol{n}}\hspace{1.65cm}\mbox{on}\ \Gamma . \end{align*}$$

We proceed by separation of variables. Assuming that $\widehat{p}_i(x,y)=X_i(x)Y_i(y)$, we are left with the following problem:

$$\begin{align} \displaystyle \frac{X_i''(x)}{X_i(x)}+\frac{Y_i''(y)}{Y_i(y)}&= \frac{\lambda ^2\rho _i}{\rho _i c_i^2+2\nu _i\lambda }\quad \mbox{in}\ \Omega _i,{\tag{5.2}}\\ X_i'(0)=X_i'(A)&=0,\hspace{1.1cm} \qquad i=1,2,{\tag{5.3}}\\ Y_1'(0)=Y_2'(B)&=0,{\tag{5.4}}\\ \frac{1}{\rho _1}X_1(x)Y_1'(H)&=\frac{1}{\rho _2}X_2(x)Y_2'(H),\hspace{0.2cm}0<x<A,{\tag{5.5}}\\ X_1(x)Y_1(H)&=X_2(x)Y_2(H),\quad 0<x<A. {\tag{5.6}} \end{align}$$

From 5.2 we have that $X_i(x)''/X_i(x)$ and $Y_i(y)''/Y_i(y)$ are constant. Moreover, from 5.5 and 5.6, it is easy to check that $Y_i(H)$ and $Y'_i(H)$ cannot vanish simultaneously and $X_1(x)=X_2(x)$ (actually, it is derived that $X_1(x)=CX_2(x)$, but the constant $C$ can be chosen equal to one without loss of generality).

From the fact that $X_i(x)''/X_i(x)$ is constant and 5.3, we have that

$$\begin{equation*} X_1(x)=X_2(x)=\cos \left(\frac{m\pi x}{A}\right),\qquad m=0,1,2,\ldots . \end{equation*}$$

On the other hand, from the fact that $Y_i(y)''/Y_i(y)$ is also constant and 5.4 we derive

$$\begin{equation} Y_1(y)=C_1\cosh (r_m^{(1)}(\lambda )y)\qquad \text{and}\qquad Y_2(y)=C_2\cosh (r_m^{(2)}(\lambda )(y-B)), {\tag{5.7}} \end{equation}$$

where $C_1$ and $C_2$ are constants and

$$\begin{equation*} r_m^{(i)}:=\sqrt {\frac{\lambda ^2\rho _i}{\rho _i c_i^2+2\nu _i\lambda }+\frac{m^2\pi ^2}{A^2}}, \qquad m=0,1,2,\ldots , \qquad i=1,2. \end{equation*}$$

Since $Y_i(H)$ and $Y_i'(H)$ cannot vanish simultaneously, 5.5 and 5.6 lead to

$$\begin{equation*} \frac{1}{\rho _1}Y_1'(H)=\frac{1}{\rho _{2}}Y_2'(H)\qquad \mbox{and}\qquad Y_1(H)=Y_2(H), \end{equation*}$$

respectively. Thus, substituting 5.7 into these equations yields the following linear system for the coefficients $C_1$ and $C_2$:

$$\begin{align*} \displaystyle C_1\cosh (r_{m}^{(1)}(\lambda )H)&=C_2\cosh (r_{m}^{(2)}(\lambda )(H-B)),\\ \displaystyle \frac{C_1r_{m}^{(1)}(\lambda )}{\rho _1}\sinh (r_{m}^{(1)}(\lambda )H)&=\frac{C_2r_{m}^{(2)}(\lambda )}{\rho _2}\sinh (r_{m}^{(2)}(\lambda )(H-B)). \end{align*}$$

For this system to have non-trivial solutions, its determinant must vanish, which yields the following non-linear equation in $\lambda$ for $m=0,1,2,\ldots$ whose roots are the eigenvalues of Problem 1:

$$\begin{multline*} f_m(\lambda ):=\frac{r_{m}^{(1)}(\lambda )}{\rho _1}\sinh (r_{m}^{(1)}(\lambda )H)\cosh (r_{m}^{(2)}(\lambda )(H-B))\\ -\frac{r_{m}^{(2)}(\lambda )}{\rho _2}\sinh (r_{m}^{(2)}(\lambda )(H-B))\cosh (r_{m}^{(1)}(\lambda )H)=0. \end{multline*}$$

We have computed some roots of the above equation and used these roots as exact eigenvalues to compare those obtained with the method proposed in this paper. For the geometrical parameters, we have taken $A=1\, \mbox{m}$, $B=2\,\mbox{m}$, and $H=1.25\,\mbox{m}$.

We have used physical parameters of water and air for the density and acoustic speed of the fluids in $\Omega _1$ and $\Omega _2$, respectively: $c_1=1430 \hspace{0.1cm}\mbox{m/s}$, $\rho _1=1000\hspace{0.1cm}\mbox{ kg/m}^3$, $c_2=340\hspace{0.1cm}\mbox{m/s}$, and $\rho _2=1\hspace{0.1cm}\mbox{kg/m}^3$. We have used uniform meshes as those shown in Figure 3. The refinement parameter $N$ refers to the number of elements per width of the rectangle.

In the presence of dissipation ($\nu \ne 0$), the eigenvalues $\lambda$ are complex numbers $\lambda =\eta +i\omega$, with $\eta <0$ the decay rate and $\omega$ the vibration frequency. In the absence of dissipation ($\nu =0$), the eigenvalues $\lambda$ are purely imaginary ($\eta =0$). The same holds for the computed eigenvalues $\lambda _h$.

In our first test, we neglected the viscosity damping effects by taking $\nu _1=\nu _2=0$. In this case, the eigenvalues $\lambda$ are actually purely imaginary. Table 1 shows the four smallest eigenvalues computed with the proposed method on successively refined meshes. Accurate values of the zeros of $f_m(\lambda )$ obtained with the MATLAB command $\mathtt{{fminsearch}}$ applied to $|f_m(\lambda )|$ are also reported on the last line of the table as exact eigenvalues.

As predicted by the theory, these eigenvalues are purely imaginary. The high accuracy of the computed eigenvalues can be observed from Table 1 even for the coarsest mesh. We have used a least squares fitting to estimate the convergence rate for each eigenvalue, which are also reported in Table 1. A clear order $\mathcal{O}(h^2)$ can be seen in all cases.

Secondly, we have used the same physical parameters as above for both fluids, but considering now non-vanishing viscosities. In order to make the dissipation effects more visible, we have used unrealistically large viscosity values: $\nu _1=9\,\mbox{N/m$s^2$}$ and $\nu _2=1\,\mbox{N/m$s^2$}$. We have repeated the scheme used above. We report in Table 2 the computed and “exact” eigenvalues and the estimated convergence rates, which are in accordance with the theory once again. Notice that now all $\lambda$ have negative real parts (decay rate) as predicted by the theory.

It can be seen from Table 2 that even in the coarsest mesh the vibration frequencies are computed with at least three correct significant digits. In turn, the decay rates are computed with at least two correct significant digits, in spite of the fact that they are much smaller than the vibration frequencies. Moreover, a least square fitting of the computed decay rates shows that they also converge quadratically. As a consequence, the decay rates computed with the finest mesh attain three or four significant correct digits.

Finally, Figure 4 shows the real and imaginary parts of the computed pressure as defined in 2.1 for the smallest eigenvalue.

Test 2.

As a second test, we have applied our code to solve a problem with a curved interface. In spite of the fact that such a problem does not lie in our theoretical framework, we will report experimental results which will allow us to assess the performance of the method in this case.

The whole domain $\Omega$ is the same as in the previous experiment, but with a curved interface $\Gamma$. In particular, we have taken as $\Gamma$ an arc of a circle with center at the point $(0.5,2.0)$ and endpoints $(0.0,1.25)$ and $(1.0,1.25)$. We report in Table 3 the computed eigenvalues. In this case, there is no analytical solution available. Therefore, we have obtained more accurate approximations of the exact eigenvalues by means of a least square fitting. These more accurate values are reported on the last line of Table 3 as “Exact”. We also report in this table the estimated order of convergence, which once more is clearly quadratic.

Finally, Figure 5 shows the real and imaginary parts of the computed pressure as defined in 2.1 for the smallest eigenvalue. The curved interface can be clearly appreciated in the imaginary part of the pressure.

Figure 1.

2D sketch of the polygonal domains for the fluids.

Figure 2.

Test 1. Two fluids in a rectangular rigid cavity.

Figure 3.

Test 1. Meshes for $N=4$ (left) and $N=8$ (right).

Table 1.

Test 1. Computed and exact eigenvalues for dissipative fluids in a rigid cavity.

$m$	$1$	$0$	$1$	$0$
$N=8$	$1066.07\,i$	$1418.42\,i$	$1784.37\,i$	$1796.61\,i$
$N=16$	$1067.78\,i$	$1422.52\,i$	$1781.49\,i$	$1797.09\,i$
$N=32$	$1068.21\,i$	$1423.54\,i$	$1780.73\,i$	$1797.21\,i$
$N=64$	$1068.33\,i$	$1423.79\,i$	$1780.55\,i$	$1797.23\,i$
Order	$2.00$	$2.00$	$1.99$	$2.00$
Exact	$1068.36\,i$	$1423.87\,i$	$1780.49\,i$	$1797.24\,i$

Table 2.

Test 1. Computed and exact eigenvalues for dissipative fluids in a rigid cavity.

$m$	$1$	$0$	$1$	$0$
$N=8$	$-9.83127 + 1066.03\,i$