# Numerical solution of the cavity scattering problem for flexural waves on
thin plates: linear finite element methods

Junhong Yue  College of Computer Science and Technology (College of Data
Science), Taiyuan University of Technology, Shanxi 030024, China
<EMAIL_ADDRESS>and  Peijun Li  Department of Mathematics, Purdue
University, West Lafayette, Indiana 47907, USA<EMAIL_ADDRESS>
###### Abstract.

Flexural wave scattering plays a crucial role in optimizing and designing
structures for various engineering applications. Mathematically, the flexural
wave scattering problem on an infinite thin plate is described by a fourth-
order plate-wave equation on an unbounded domain, making it challenging to
solve directly using the regular linear finite element method (FEM). In this
paper, we propose two numerical methods, the interior penalty FEM (IP-FEM) and
the boundary penalty FEM (BP-FEM) with a transparent boundary condition (TBC),
to study flexural wave scattering by an arbitrary-shaped cavity on an infinite
thin plate. Both methods decompose the fourth-order plate-wave equation into
the Helmholtz and modified Helmholtz equations with coupled conditions at the
cavity boundary. A TBC is then constructed based on the analytical solutions
of the Helmholtz and modified Helmholtz equations in the exterior domain,
effectively truncating the unbounded domain into a bounded one. Using linear
triangular elements, the IP-FEM and BP-FEM successfully suppress the
oscillation of the bending moment of the solution at the cavity boundary,
demonstrating superior stability and accuracy compared to the regular linear
FEM when applied to this problem.

###### Key words and phrases:

Biharmonic wave equation, flexural wave scattering problem, transparent
boundary condition, linear finite element method

###### 2010 Mathematics Subject Classification:

65N30, 74J20, 74S05

The first author is supported by the NSFC grant 11901423. The second author is
supported in part by the NSF grant DMS-2208256.

##  1\. Introduction

Flexural wave scattering is a widespread phenomenon with practical engineering
applications, including the design of lightweight mechanical structures with
low noise [1], ultra-broadband elastic cloaking devices [13, 18], Platonic
diffraction gratings and arrays [22, 23], massive floating concrete runways
offshore, and health monitoring of thin-walled structures such as aircraft
wings, oil tank walls, and pressure vessels [41]. This phenomenon arises from
the interaction of incident waves with different types of scatterers (e.g.,
voids, rigid, and elastic scatterers) on thin-wall structures [33], attracting
considerable attention in the research community. For instance, Akrucci et al.
[1] investigated the effect of acoustic black holes on flexural wave
scattering on infinite thin plates, effectively reducing plate oscillations
without increasing structure mass. Liu et al. [36] realized a broadband
cylindrical cloak for flexural waves in elastic thin plates using nonlinear
transformation, guiding flexural waves more effectively outside the cloak
region. Haslinger et al. [24] studied scattering and transmission of flexural
waves in a thin plate with a semi-infinite array of point scatterers,
demonstrating dynamically anisotropic wave effects in semi-infinite platonic
crystals. Evans et al. [17] explored flexural wave scattering in an elastic
thin plate floating on water. Wang et al. [42] investigated the scattering
behavior of extensional and flexural plate waves by a cylindrical
inhomogeneity for structural health monitoring, characterizing the interaction
of plate waves with structural damage. Consequently, studying flexural wave
scattering in thin-walled structures holds great importance for optimizing and
designing new structures.

Mathematically, the scattering problem of flexural waves on infinite thin
plates can be described by a fourth-order plate-wave equation on an unbounded
domain. Analytical solutions for such problems are only attainable for
isotropic thin plates containing scatterers with simple geometries, posing
challenges for complex geometries and media. Apparently, numerical methods are
necessary to simulate and solve these problems. Several numerical approaches
have been proposed for solving flexural wave scattering problems. For
instance, Norris et al. [38] explored the energy flux conservation and the
optical theorem in the context of flexural wave scattering, applying them to
flexural scattering by circular scatterers. Matus et al. [37] employed the
transfer matrix (T-matrix) method for flexural wave scattering by a single
noncircular scatterer. Climente et al. [9, 10] utilized T-matrix and impedance
matrix methods to study flexural wave scattering by a hole containing beam
resonators in an infinite thin plate for vibration control. Lee et al. [30,
31, 32] investigated flexural wave scattering in a thin plate with multiple
circular inclusions using the multipole method, the multipole Trefftz method,
and the null-field integral equation approach. Cai et al. [6] applied the
T-matrix method for the multiple scattering of flexural waves by different
types of circular scatterers on thin plates. Wang et al. [43] studied the
multiple scattering of flexural waves by varying-thickness annular inclusions
on infinite thin plates using a semi-analytical method. Wang et al. [44]
developed a semi-analytical model for a novel plate-harvester system,
combining multiple flexural scattering theories of thin plates with coupled
electroelastic dynamics of piezoelectric composite beams. Dong et al. [14]
proposed a novel formulation of boundary integral equations for the scattering
of flexural waves by obstacles on infinite thin plates.

However, the finite element method (FEM), known for its efficiency and
stability in numerical algorithms, is rarely employed to solve the scattering
problem of flexural waves due to the challenges posed by the unbounded domain
of the problem and the presence of high-order partial differential equations
(PDEs). To address the unbounded domain issue, common truncated techniques,
such as absorbing boundary condition (ABC) [11, 16], perfectly matched layer
(PML) [4], and transparent boundary condition (TBC) [27, 28, 34], are used to
truncate the unbounded domain into bounded computational domains when applying
FEM. However, these techniques are primarily constructed for second-order PDEs
in exterior domains, such as acoustic, electromagnetic, and elastic wave
scattering problems. Directly applying these techniques to flexural wave
scattering problems, which involve a fourth-order wave equation, is
challenging. Therefore, the development of an effective truncation technique
specifically for flexural wave scattering is crucial. Moreover, to solve the
fourth-order wave equation within the truncation domain, various numerical
methods have been proposed, including meshless methods [40] and different
types of FEMs. In the literature, classical FEM approaches include $C^{1}$
conforming FEMs like the Argyris element [8] with 21 degrees of freedom, as
well as nonconforming FEMs like the Adini element [25] and Morley element
[45]. However, these elements are rarely practical due to either their high
number of degrees of freedom or the complexity involved in their
implementation. It remains a challenge to find a simpler and more efficient
approach to solve the fourth-order wave equation in the truncation domain.

As a result, numerous FEMs based on linear triangular elements have been
proposed to solve fourth-order problems, including mixed methods [2, 3, 5, 19,
20] and the recovery-based linear FEM [7, 21, 26, 29]. While mixed methods are
widely used, they require careful treatment of essential and natural boundary
conditions. In particular, the Ciarlet–Raviart mixed FEM exhibits oscillation
behavior of the middle variable on the boundary, necessitating the addition of
corresponding penalty terms [2]. On the other hand, the recovery-based linear
FEM is a nonconforming FEM that discretizes the Laplace operator by utilizing
the gradient recovery operator acting on the gradient of the $C^{0}$ linear
element.

This paper first reduces the biharmonic plate wave equation to the Helmholtz
and modified Helmholtz equations with coupled conditions at the cavity
boundary using two auxiliary functions. To truncate the unbounded domain into
a bounded one, the transparent boundary conditions (TBCs) are constructed
using Fourier series solutions of the Helmholtz and modified Helmholtz
equations, satisfying the Sommerfeld radiation conditions. The finite element
approximation with linear triangular elements [12] is then utilized to solve
the decomposed problem with coupled boundary conditions. However, the
solutions obtained by this method exhibit oscillatory behavior at the cavity
boundary. Motivated by recent works [2, 15, 35, 46], we introduce the interior
penalty term or the boundary penalty term to the variational formulations of
the Helmholtz and modified Helmholtz equations, effectively suppressing the
oscillation of the bending moment of the solution at the cavity boundary.

This paper focuses on the numerical computation of flexural wave scattering by
an arbitrary shaped cavity with the clamped boundary and contributes in the
following four aspects:

  1. (1)

Construction of TBCs for flexural wave cavity scattering in two dimensions,
which is equivalent to satisfying the Sommerfeld radiation conditions.

  2. (2)

Deduction of a decomposed problem of the biharmonic plate wave equation by
introducing two auxiliary functions, and proof of the uniqueness of its
solution.

  3. (3)

Proposal of the linear finite element method with interior penalty term (IP-
FEM) and boundary penalty term (BP-FEM) for solving the coupled boundary
Helmholtz and modified Helmholtz equations, providing stable numerical
solutions.

  4. (4)

Construction of an analytical solution for flexural wave scattering by a
circular cavity with a clamped boundary, facilitating comparative analysis.

In this work, we propose the IP-FEM and the BP-FEM as numerical methods to
simulate cavity scattering in an infinite thin plate. The paper is outlined as
follows. Section 2 describes the cavity scattering problems using the
biharmonic plate wave equation with the clamped boundary condition. In Section
3, we construct the transparent boundary conditions (TBCs) to truncate the
unbounded domain for plate wave scattering. Section 4 presents a decomposed
problem by introducing two auxiliary functions for the plate wave equation,
along with the proof of the uniqueness of this decomposition problem. In
Section 5, we provide the variational formulations and discretized systems of
the IP-FEM and the BP-FEM. Section 6 showcases numerical experiments to
validate the effectiveness of the proposed methods, comparing them with
analytical solutions or reference solutions. Finally, in Section 7, we draw
conclusions from this study.

##  2\. Problem formulation

Let us consider a cavity denoted as $D$, located within an infinite thin
plate, with a Lipschitz continuous boundary $\partial D$. The structure is
illuminated by a time-harmonic plane wave represented by $u^{\rm
inc}(\boldsymbol{x})=e^{{\rm i}\kappa\boldsymbol{x}\cdot\boldsymbol{d}}$,
where $\kappa>0$ is the wavenumber determined by $\kappa^{4}=\omega^{2}\rho
h/D_{c}$. Here, $\omega$ denotes the angular frequency, while $\rho$, $h$, and
$D_{c}$ refer to the mass density, thickness, and flexural rigidity of the
plate, respectively. The incident direction is given by
$\boldsymbol{d}=(\cos(\alpha),\sin(\alpha))$, with $\alpha\in[0,2\pi)$
representing the incident angle.

It can be verified that the incident wave field $u^{\rm inc}$ satisfies the
two-dimensional biharmonic wave equation given by:

$\Delta^{2}u^{\rm inc}-\kappa^{4}u^{\rm
inc}=0\quad\text{in}~{}\mathbb{R}^{2}.$ (2.1)

The out-of-plane displacement of the plate, denoted as $u$, satisfies the two-
dimensional biharmonic wave equation in the exterior of $D$, which is
expressed as:

$\Delta^{2}u-\kappa^{4}u=0\quad\text{in}~{}\mathbb{R}^{2}\setminus\overline{D}.$
(2.2)

The total field $u$ is assumed to satisfy the following clamped boundary
condition on $\partial D$:

$u=0,\quad\partial_{\nu}u=0,$ (2.3)

where $\nu$ is the unit normal vector on $\partial D$.

It can be observed from (2.1)–(2.3) that the scattered field $v=u-u^{\rm inc}$
satisfies:

$\Delta^{2}v-\kappa^{4}v=0\quad\text{in}~{}\mathbb{R}^{2}\setminus\overline{D},$
(2.4)

and the following boundary conditions on $\partial D$:

$v=-u^{\rm inc},\quad\partial_{\nu}v=-\partial_{\nu}u^{\rm inc}.$ (2.5)

In addition, the scattered field $v$ and its Laplacian $\Delta v$ are required
to satisfy the Sommerfeld radiation condition:

$\lim_{r\rightarrow\infty}r^{1/2}(\partial_{r}v-{\rm i}\kappa
v)=0,\quad\lim_{r\rightarrow\infty}r^{1/2}(\partial_{r}\Delta v-{\rm
i}\kappa\Delta v)=0,\quad r=|x|.$ (2.6)

We introduce standard notations used in this paper. Let
$B_{R}=\\{\boldsymbol{x}\in\mathbb{R}^{2}:|\boldsymbol{x}|<R\\}$ represent a
disc with boundary
$\Gamma_{R}=\\{\boldsymbol{x}\in\mathbb{R}^{2}:|\boldsymbol{x}|=R\\}$. The
radius $R$ is chosen to be sufficiently large such that the cavity $D$ is
completely contained within $B_{R}$. We denote
$\Omega=B_{R}\setminus\overline{D}$. The inner product and the norm in
$L^{2}(\Omega)$ are

$(\phi,\psi)=\int_{\Omega}\phi\bar{\psi}d\boldsymbol{x},\quad\|\phi\|_{0,\Omega}=(\phi,\phi)^{1/2}.$

Define $H_{\partial D}^{1}=\\{\phi\in
H^{1}(\Omega):\phi=0~{}\text{on}~{}\partial D\\}$, $H_{\Omega}^{1}=\\{\phi\in
H^{1}(\Omega):\phi=0~{}\text{on}~{}\overline{\Omega}\backslash\partial D\\}$,
and $H^{1}=H_{\partial D}^{1}\cup H_{\Omega}^{1}$. It is known that
$H^{-1/2}(\Gamma_{R})$ is the dual space of $H^{1/2}(\Gamma_{R})$ with respect
to the inner product

$\langle\phi,\psi\rangle_{\Gamma_{R}}=\int_{\Gamma_{R}}\phi\bar{\psi}ds.$

##  3\. Transparent boundary conditions

In this section, we introduce the transparent boundary condition (TBC) on
$\Gamma_{R}$ to transform the cavity scattering problem from the open domain
$\mathbb{R}^{2}\setminus\overline{D}$ to the bounded domain $\Omega$. This
allows us to truncate the unbounded domain into a bounded one for numerical
simulation and analysis.

Following [14], we consider two auxiliary functions $v_{H}$ and $v_{M}$,
defined as:

$v_{H}=-\frac{1}{2\kappa^{2}}(\Delta v-\kappa^{2}v),\quad
v_{M}=\frac{1}{2\kappa^{2}}(\Delta v+\kappa^{2}v).$

It can be observed that the quantities $v$, $\Delta v$, $v_{H}$, and $v_{M}$
are related through the equations

$v=v_{H}+v_{M},\quad\Delta v=\kappa^{2}(-v_{H}+v_{M}),$ (3.1)

and

$v_{H}=\frac{1}{2}(v-\kappa^{-2}\Delta v),\quad
v_{M}=\frac{1}{2}(v+\kappa^{-2}\Delta v).$ (3.2)

The biharmonic wave equation (2.4) can be written as

$(\Delta^{2}-\kappa^{4})v=(\Delta-\kappa^{2})(\Delta+\kappa^{2})(v_{H}+v_{M})=0\quad\text{in}~{}\mathbb{R}^{2}\setminus\overline{D},$

which implies that $v_{H}$ and $v_{M}$ satisfy the Helmholtz equation and the
modified Helmholtz equation, respectively:

$\Delta v_{H}+\kappa^{2}v_{H}=0,\quad\Delta v_{M}-\kappa^{2}v_{M}=0.$ (3.3)

Combining (2.6) and (3.3), we deduce that the functions $v_{H}$ and $v_{M}$
satisfy the Sommerfeld radiation condition:

$\lim_{r\rightarrow\infty}r^{1/2}(\partial_{r}v_{H}-{\rm i}\kappa
v_{H})=0,\quad\lim_{r\rightarrow\infty}r^{1/2}(\partial_{r}v_{M}-{\rm i}\kappa
v_{M})=0.$ (3.4)

From (3.3) and (3.4), we can conclude that $v_{H}$ and $v_{M}$ have the
following Fourier series expansions in
$\mathbb{R}^{2}\setminus\overline{B_{R}}$:

$v_{H}(r,\theta)=\sum_{n\in\mathbb{Z}}\frac{H_{n}^{(1)}(\kappa
r)}{H_{n}^{(1)}(\kappa R)}v_{H}^{(n)}(R)e^{{\rm i}n\theta},\quad
v_{M}(r,\theta)=\sum_{n\in\mathbb{Z}}\frac{K_{n}(\kappa r)}{K_{n}(\kappa
R)}v_{M}^{(n)}(R)e^{{\rm i}n\theta},$

where $H_{n}^{(1)}$ is the Hankel function of the first kind with order $n$,
$K_{n}$ is the modified Bessel function with order $n$, $v_{H}^{(n)}$ and
$v_{M}^{(n)}$ are the Fourier coefficients given by

$v_{H}^{(n)}(R)=\frac{1}{2\pi}\int_{0}^{2\pi}v_{H}(R,\theta)e^{-{\rm
i}n\theta}d\theta,\quad
v_{M}^{(n)}(R)=\frac{1}{2\pi}\int_{0}^{2\pi}v_{M}(R,\theta)e^{-{\rm
i}n\theta}d\theta.$

For any function $\phi\in L^{2}(\Gamma_{R})$ with the Fourier series expansion

$\phi(R,\theta)=\sum_{n\in\mathbb{Z}}\phi^{(n)}(R)e^{{\rm
i}n\theta},\quad\phi^{(n)}(R)=\frac{1}{2\pi}\int_{0}^{2\pi}\phi(R,\theta)e^{-{\rm
i}n\theta}d\theta,$

we define two Dirichlet-to-Neumann (DtN) operators

$\displaystyle T_{1}\phi=\frac{1}{R}\sum_{n\in\mathbb{Z}}h_{n}(\kappa
R)\phi^{(n)}(R)e^{{\rm i}n\theta},\quad
T_{2}\phi=\frac{1}{R}\sum_{n\in\mathbb{Z}}k_{n}(\kappa R)\phi^{(n)}(R)e^{{\rm
i}n\theta},$ (3.5)

where

$h_{n}(z)=z\frac{H_{n}^{(1)^{\prime}}(z)}{H_{n}^{(1)}(z)},~{}~{}~{}k_{n}(z)=z\frac{K^{\prime}_{n}(z)}{K_{n}(z)}.$

Using (3.5), we deduce the TBC on $\Gamma_{R}$:

$\partial_{r}v_{H}=T_{1}v_{H},\quad\partial_{r}v_{M}=T_{2}v_{M}.$ (3.6)

Given that $v$ and $\Delta v$ are periodic functions with respect to $\theta$
on $\Gamma_{R}$, they can be represented by their Fourier series expansions:

$\displaystyle v(R,\theta)=\sum_{n\in\mathbb{Z}}v^{(n)}(R)e^{{\rm
i}n\theta},\quad\Delta v(R,\theta)=\sum_{n\in\mathbb{Z}}(\Delta
v)^{(n)}(R)e^{{\rm i}n\theta}.$

By combining (3.6) and (3.1)–(3.2), we obtain the TBC for the scattered
fields:

$\displaystyle\partial_{r}v$ $\displaystyle=T_{1}v_{H}+T_{2}v_{M}$
$\displaystyle=\frac{1}{2}T_{1}(v-\kappa^{-2}\Delta
v)+\frac{1}{2}T_{2}(v+\kappa^{-2}\Delta v)$
$\displaystyle=\frac{1}{2}(T_{1}+T_{2})v+\frac{1}{2\kappa^{2}}(T_{2}-T_{1})\Delta
v$

and

$\displaystyle\partial_{r}\Delta v$
$\displaystyle=\kappa^{2}(-T_{1}v_{H}+T_{2}v_{M})$
$\displaystyle=\kappa^{2}(-\frac{1}{2}T_{1}(v-\kappa^{-2}\Delta
v)+\frac{1}{2}T_{2}(v+\kappa^{-2}\Delta v))$
$\displaystyle=\frac{\kappa^{2}}{2}(T_{2}-T_{1})v+\frac{1}{2}(T_{1}+T_{2})\Delta
v.$

Using the above equations, we deduce the TBC for the total fields:

$\partial_{r}u=\frac{1}{2}(T_{1}+T_{2})u+\frac{1}{2\kappa^{2}}(T_{2}-T_{1})\Delta
u+g_{1}$ (3.7)

and

$\partial_{r}\Delta
u=\frac{\kappa^{2}}{2}(T_{2}-T_{1})u+\frac{1}{2}(T_{1}+T_{2})\Delta u+g_{2},$
(3.8)

where

$g_{1}=\partial_{r}u^{\rm inc}-T_{1}u^{\rm inc},\quad g_{2}=-\kappa^{2}g_{1}.$

##  4\. The decomposed problem

In this section, we present a decomposed formulation for the cavity scattering
problem, consisting of (2.2)–(2.3) and (3.7)–(3.8).

Let us consider two auxiliary functions defined as follows:

$\displaystyle p=\frac{1}{2\kappa^{2}}(\Delta u-\kappa^{2}u),\quad
q=\frac{1}{2\kappa^{2}}(\Delta u+\kappa^{2}u).$ (4.1)

It is evident that $p$ satisfies the Helmholtz equation, while $q$ satisfies
the modified Helmholtz equation. Moreover, it can be verified that

$\displaystyle u=q-p,\quad\Delta u=\kappa^{2}(p+q),$ (4.2)

and

$\displaystyle p=-u_{H},\quad q=u_{M},$ (4.3)

where $u_{H}$ and $u_{M}$ are the Helmholtz and modified Helmholtz wave
components of $u$, respectively.

Using (4.1)–(4.3), we obtain the following boundary value problem for $p$ and
$q$:

$\left\\{\begin{aligned} &\Delta p+\kappa^{2}p=0,\quad\Delta
q-\kappa^{2}q=0&\text{in}~{}\Omega,\\\
&p-q=0,\quad\partial_{\nu}p-\partial_{\nu}q=0&\text{on}~{}\partial D,\\\
&\partial_{r}p=T_{1}p-g_{1},\quad\partial_{r}q=T_{2}q&\text{on}~{}\Gamma_{R}.\end{aligned}\right.$
(4.4)

Equivalently, we may consider two auxiliary functions for the scattered field:

$p^{s}=\frac{1}{2\kappa^{2}}(\Delta v-\kappa^{2}v),\quad
q^{s}=\frac{1}{2\kappa^{2}}(\Delta v+\kappa^{2}v).$

Hence, we have

$v=q^{s}-p^{s},\quad\Delta v=\kappa^{2}(p^{s}+q^{s}),$

and

$p^{s}=-v_{H},\quad q^{s}=v_{M}.$

It can be verified that $p^{s}$ and $q^{s}$ satisfy the following boundary
value problem:

$\left\\{\begin{aligned} &\Delta p^{s}+\kappa^{2}p^{s}=0,\quad\Delta
q^{s}-\kappa^{2}q^{s}=0&\text{in}~{}\Omega,\\\ &p^{s}-q^{s}=u^{\rm
inc},\quad\partial_{\nu}p^{s}-\partial_{\nu}q^{s}=\partial_{\nu}u^{\rm
inc}&\text{on}~{}\partial D,\\\
&\partial_{r}p^{s}=T_{1}p^{s},\quad\partial_{r}q^{s}=T_{2}q^{s}&\text{on}~{}\Gamma_{R}.\end{aligned}\right.$
(4.5)

######  Lemma 4.1.

Let $z$ be a positive real number. Then

$\displaystyle\Re(h_{n}(z))<0,\quad\Im(h_{n}(z))>0,\quad\Re(k_{n}(z))<0,\quad\Im(k_{n}(z))=0.$

###### Proof.

Using the definition $H_{n}^{(1)}(z)=J_{n}(z)+{\rm i}Y_{n}(z)$, we can express
$h_{n}(z)$ into

$\displaystyle
h_{n}(z)=\frac{z}{|H_{n}^{(1)}(z)|^{2}}\left(J_{n}^{\prime}(z)+{\rm
i}Y_{n}^{\prime}(z)\right)\left(J_{n}(z)-{\rm i}Y_{n}(z)\right),$

where the real-valued functions $J_{n}(z)$ and $Y_{n}(z)$ are the Bessel
functions of the first kind and second kind with order $n$, respectively.

First, we consider the real part of $h_{n}(z)$:

$\displaystyle\Re(h_{n}(z))$
$\displaystyle=\frac{z}{|H_{n}^{(1)}(z)|^{2}}\left(J_{n}^{\prime}(z)J_{n}(z)+Y_{n}^{\prime}(z)Y_{n}(z)\right)$
$\displaystyle=\frac{z}{2|H_{n}^{(1)}(z)|^{2}}\frac{d}{dz}\left(J_{n}(z)^{2}+Y_{n}(z)^{2}\right).$

By the Nicholson’s integral [39, $(10.9.30)$],

$J_{n}(z)^{2}+Y_{n}(z)^{2}=\frac{8}{\pi^{2}}\int_{0}^{\infty}\cosh(2nt)K_{0}(2z\sinh(t))dt,$

where $\sinh(t)$ and $\cosh(t)$ are the hyperbolic sine and hyperbolic cosine
functions, respectively, we have

$\frac{d}{dz}\left(J_{n}(z)^{2}+Y_{n}(z)^{2}\right)=\frac{16}{\pi^{2}}\int_{0}^{\infty}\cosh(2nt)\sinh(t)K_{0}^{\prime}(2z\sinh(t))dt.$

Given that $\cosh(t)>0$ for $t\in\mathbb{R}$ and $\sinh(t)>0$ for $t>0$, along
with the fact that $K_{n}(z)$ is positive and decreasing throughout the
interval $0<z<\infty$ for $n\geq 0$ (cf. [39, Section 10.37]), and
$K_{0}^{\prime}(z)=-K_{1}(z)<0$, we can conclude that $\Re(h_{n}(z))<0$.

Second, we consider the imaginary part of $h_{n}(z)$:

$\displaystyle\Im(h_{n}(z))=\frac{z}{|H_{n}^{(1)}(z)|^{2}}\left(Y_{n}^{\prime}(z)J_{n}(z)-J_{n}^{\prime}(z)Y_{n}(z)\right).$

Using the identities

$2J_{n}^{\prime}=J_{n-1}(z)-J_{n+1}(z),\quad
2Y_{n}^{\prime}(z)=Y_{n-1}(z)-Y_{n+1}(z),$

and the Wronskian [39]:

$W\left\\{J_{n}(z),Y_{n}(z)\right\\}=J_{n+1}(z)Y_{n}(z)-J_{n}(z)Y_{n+1}(z)=2/(\pi
z),$

we have

$\displaystyle\Im(h_{n}(z))$
$\displaystyle=\frac{z}{|H_{n}^{(1)}(z)|^{2}}\left((Y_{n-1}(z)-Y_{n+1}(z))J_{n}(z)-(J_{n-1}(z)-J_{n+1}(z))Y_{n}(z)\right).$
$\displaystyle=\frac{2}{\pi}\frac{1}{|H_{n}^{(1)}(z)|^{2}}>0.$

Next, we examine the properties of $\Re(k_{n}(z))$ and $\Im(k_{n}(z))$. Since
$K_{n}(z)$ is a real-valued function, $K_{n}^{\prime}(z)$ is also real,
implying that $\Im(k_{n}(z))=0$. For a given $n$ and $z>0$, $K_{n}(z)>0$, and
$K_{n}(z)$ is monotonically decreasing with respect to $z$, i.e.,
$K^{\prime}_{n}(z)<0$. Consequently, we have $\Re(k_{n}(z))<0$ for $z>0$. ∎

######  Theorem 4.2.

The coupled boundary value problem (4.4) has at most one solution for
$\kappa>0$.

###### Proof.

It suffices to show that $p=0$ and $q=0$ in $\Omega$ when $g_{1}=0$. Applying
Green’s theorem in $\Omega$ and the boundary condition, we obtain

$\displaystyle(\nabla p,\nabla p)-\kappa^{2}(p,p)-\langle
T_{1}p,p\rangle_{\Gamma_{R}}-\langle\partial_{\nu}p,p\rangle_{\partial D}=0,$
$\displaystyle(\nabla q,\nabla q)+\kappa^{2}(q,q)-\langle
T_{2}q,q\rangle_{\Gamma_{R}}-\langle\partial_{\nu}q,q\rangle_{\partial D}=0.$

Since $\langle\partial_{\nu}p,p\rangle_{\partial
D}=\langle\partial_{\nu}q,q\rangle_{\partial D}$ on $\partial D$, we have

$\displaystyle(\nabla p,\nabla p)-\kappa^{2}(p,p)-\langle
T_{1}p,p\rangle_{\Gamma_{R}}=(\nabla q,\nabla q)+\kappa^{2}(q,q)-\langle
T_{2}q,q\rangle_{\Gamma_{R}}.$

A simple calculation yields

$\displaystyle\langle
T_{1}p,p\rangle_{\Gamma_{R}}=2\pi\sum_{n\in\mathbb{Z}}h_{n}(\kappa
R)|p^{(n)}|^{2},\quad\langle
T_{2}q,q\rangle_{\Gamma_{R}}=2\pi\sum_{n\in\mathbb{Z}}k_{n}(\kappa
R)|q^{(n)}|^{2},$

where $p^{(n)}$ and $q^{(n)}$ are the Fourier coefficients of $q$ and $q$ on
$\Gamma_{R}$. Taking the imaginary part of the above equation gives

$\displaystyle\Im\\{-\langle T_{1}p,p\rangle_{\Gamma_{R}}+\langle
T_{2}q,q\rangle_{\Gamma_{R}}\\}$
$\displaystyle=-2\pi\sum_{n\in\mathbb{Z}}\Im(h_{n}(\kappa R))|p^{(n)}|^{2}$
$\displaystyle\quad+2\pi\sum_{n\in\mathbb{Z}}\Im(k_{n}(\kappa
R))|q^{(n)}|^{2}=0,$

which gives $p^{(n)}=0$ for $n\in\mathbb{Z}$ using Lemma 4.1. Thus we have
$p=0$ and $\partial_{r}p=0$ on $\Gamma_{R}$. According to the Holmgren
uniqueness theorem, we obtain $p=0$ and $\partial_{\nu}p=0$ in
$\mathbb{R}^{2}\setminus\overline{B_{R}}$. Furthermore, a unique continuation
result implies that $p=0$ and $\partial_{\nu}p=0$ in $\Omega$. Considering the
boundary conditions on $\partial D$, we find that $q=0$ and
$\partial_{\nu}q=0$ on $\partial D$. Consequently, by applying the Holmgren
uniqueness theorem, we can deduce that $q=0$ in $\Omega$. ∎

To solve the decomposed problem (4.4) by using the linear FEM, we introduce
its variational formulation. Using the test functions ${\phi,\psi,\varphi}\in
H_{\partial D}^{1}(\Omega)\times H_{\partial D}^{1}(\Omega)\times
H_{\Omega}^{1}$, the weak formulation of (4.4) aims to find ${p,q}\in
H^{1}(\Omega)\times H^{1}(\Omega)$ that satisfy the following equations:

$\left\\{\begin{aligned} &b_{1}(p,\phi)=-\langle
g_{1},\phi\rangle_{\Gamma_{R}}&\forall~{}\phi\in H_{\partial
D}^{1}(\Omega),\\\ &b_{2}(q,\psi)=0&\forall~{}\psi\in H_{\partial
D}^{1}(\Omega),\\\ \end{aligned}\right.$ (4.6)

and

$(\nabla(p-q),\nabla\varphi)-\kappa^{2}((p+q),\varphi)=0\quad\forall~{}\varphi\in
H_{\Omega}^{1},$ (4.7)

where $p=p_{0}+p_{D}$ and $q=q_{0}+p_{D}$ with $p_{0},q_{0}\in H_{\partial
D}^{1}(\Omega)$ and $p_{D}\in H_{\Omega}^{1}$. Here the sesquilinears
$b_{1}:H^{1}\times H_{\partial D}^{1}\rightarrow\mathbb{C}$ and
$b_{2}:H^{1}\times H_{\partial D}^{1}\rightarrow\mathbb{C}$ are defined by

$\displaystyle b_{1}(\phi,\psi)$
$\displaystyle=(\nabla\phi,\nabla\psi)-\kappa^{2}(\phi,\psi)-\langle
T_{1}\phi,\psi\rangle_{\Gamma_{R}},$ $\displaystyle b_{2}(\phi,\psi)$
$\displaystyle=(\nabla\phi,\nabla\psi)+\kappa^{2}(\phi,\psi)-\langle
T_{2}\phi,\psi\rangle_{\Gamma_{R}}.$

##  5\. The linear finite element methods

In this section, we introduce the IP-FEM and BP-FEM methods for solving the
problem (4.6)–(4.7). First, we define the linear finite element spaces and the
corresponding symbols. Next, we construct the variational formulations by
incorporating an interior penalty term and a boundary penalty term,
respectively. Finally, we present the discretized systems using the linear
FEM.

###  5.1. Finite element spaces

Let $\mathcal{M}_{h}$ be a triangulation of $\Omega$ such that
$\overline{\Omega}=\cup_{K\in\mathcal{M}_{h}}K$, where $K$ denotes a
triangular element. Let $\mathcal{C}_{h}^{I}$ and $\mathcal{C}_{h}^{B}$ be the
set of all interior and boundary edges of mesh $\mathcal{M}_{h}$,
respectively.

We define the finite element space using piecewise linear functions, denoted
as $\mathbb{P}_{1}$, associated with $\mathcal{M}_{h}$. We consider the
following discrete spaces:

$S_{h}=\\{\phi_{h}\in
C(\overline{\Omega}):\phi_{h}|_{K}\in\mathbb{P}_{1}(K)~{}\forall
K\in\mathcal{M}_{h}\\},$

where $S_{h}^{0}=S_{h}\cap H_{\partial D}^{1}$ and $S_{h}^{\Omega}=S_{h}\cap
H_{\Omega}^{1}$. Both $S_{h}^{0}$ and $S_{h}^{\Omega}$ are subspaces of
$S_{h}$ that have vanishing degrees of freedom (DoFs) on $\partial D$ and
$\Omega\cup\Gamma_{R}$, respectively.

###  5.2. The variational formulation for IP-FEM

We derive the variational formulation with an interior penalty term for the
problem (4.6)–(4.7). To facilitate the formulation, we assign a unique index
$i_{K}\in\mathbb{N}$ to each element $K\in\mathcal{M}_{h}$. Furthermore, we
define the jump of a function $\phi$ across an interior edge $e=\partial
K\cap\partial K^{\prime}$ as follows:

$[\phi]_{e}:=\left\\{\begin{aligned}
&\phi|_{K}-\phi|_{K^{\prime}}&\text{if}~{}~{}i_{K}>i_{K^{\prime}},\\\
&\phi|_{K^{\prime}}-\phi|_{K}&\text{if}~{}~{}i_{K}<i_{K^{\prime}}.\end{aligned}\right.$

For any functions $\phi,\psi\in S_{h}$, we define the sesquilinear form of the
interior Neumann penalty by

$J(\phi,\psi):=\sum_{e\in\mathcal{C}_{h}^{I}}\gamma_{e}h_{e}\langle[\partial_{\nu}\phi],[\partial_{\nu}\psi]\rangle_{e},$

where $h_{e}$ is the length of interior edge $e$ and $\gamma_{e}$ is a real
positive parameter.

The sesquilinear forms $b^{h}_{1}:S_{h}\times S_{h}^{0}\rightarrow\mathbb{C}$
and $b^{h}_{2}:S_{h}\times S_{h}^{0}\rightarrow\mathbb{C}$ are defined by

$\displaystyle b^{h}_{1}(\phi,\psi)=b_{1}(\phi,\psi)-J(\phi,\psi),\quad
b^{h}_{2}(\phi,\psi)=b_{2}(\phi,\psi)+J(\phi,\psi).$

It is important to note that the sign of the penalty term should be consistent
with that of the lower-order term (i.e., the mass matrix term) in the
variational formulation. This consistency ensures enhanced stability of the
solution for discrete systems constructed using linear finite elements, from a
numerical computational perspective.

The variational formulation with an interior penalty term for problem
(4.6)–(4.7) is to find $\\{p_{h},q_{h}\\}\in S_{h}\times S_{h}$ such that

$\left\\{\begin{aligned} &b^{h}_{1}(p_{h},\phi_{h})=-\langle
g_{1},\phi_{h}\rangle_{\Gamma_{R}}&\forall~{}\phi_{h}\in S_{h}^{0},\\\
&b^{h}_{2}(q_{h},\psi_{h})=0&\forall~{}\psi_{h}\in
S_{h}^{0},\end{aligned}\right.$ (5.1)

and

$(\nabla(p_{h}-q_{h}),\nabla\varphi_{h})-\kappa^{2}(p_{h}+q_{h},\varphi_{h})-J(p_{h}+q_{h},\varphi_{h})=0\quad\forall~{}\varphi_{h}\in
S_{h}^{\Omega},$ (5.2)

where $p_{h}=p^{h}_{0}+p^{h}_{D}$ and $q_{h}=q^{h}_{0}+p^{h}_{D}$ with
$p^{h}_{0},q^{h}_{0}\in S_{h}^{0}$ and $p^{h}_{D}\in S_{h}^{\Omega}$.

###  5.3. The variational formulation for BP-FEM

We now establish the variational formulation with a boundary penalty term for
the problem (4.6)–(4.7). Consider any functions $\phi,\psi\in S_{h}^{\Omega}$,
the sesquilinear form of the boundary penalty term is defined as:

$G(\phi,\psi):=\sum_{e\in\mathcal{C}_{h}^{B}}\eta_{e}h_{e}\langle\partial_{\tau}\phi,\partial_{\tau}\psi\rangle_{e},$

where $h_{e}$ represents the length of the boundary edge $e$, $\eta_{e}$ is a
positive real parameter, and $\tau$ is the unit tangent vector on the boundary
edge $e$.

The variational formulation with a boundary penalty term for problem
(4.6)–(4.7) is defined as follows: find $\\{p_{h},q_{h}\\}\in S_{h}\times
S_{h}$ such that

$\left\\{\begin{aligned} &b_{1}(p_{h},\phi_{h})=-\langle
g_{1},\phi_{h}\rangle_{\Gamma_{R}}&\forall~{}\phi_{h}\in S_{h}^{0},\\\
&b_{2}(q_{h},\psi_{h})=0&\forall~{}\psi_{h}\in S_{h}^{0},\end{aligned}\right.$
(5.3)

and

$(\nabla(p_{h}-q_{h}),\nabla\varphi_{h})-\kappa^{2}(p_{h}+q_{h},\varphi_{h})-G(p^{h}_{D},\varphi_{h})=0\quad\forall~{}\varphi_{h}\in
S_{h}^{\Omega},$ (5.4)

where $p_{h}=p^{h}_{0}+p^{h}_{D}$ and $q_{h}=q^{h}_{0}+p^{h}_{D}$ with
$p^{h}_{0},q^{h}_{0}\in S_{h}^{0}$ and $p^{h}_{D}\in S_{h}^{\Omega}$.

Similarly, it is essential to ensure the consistency of the sign of the
penalty term $G(p^{h}_{D},\varphi_{h})$ with that of the lower-order term
$(p_{h}+q_{h},\varphi_{h})$ in the variational formulation with a boundary
penalty.

###  5.4. The discretized problems

Next, we proceed to discretize the variational problem with the interior
penalty term (5.1)–(5.2) and the boundary penalty term (5.3)–(5.4) using
linear FEM. Subsequently, we express these equations in matrix form.

Let $\\{\alpha_{j}\\}_{j=1}^{N_{h}^{I}}$ and
$\\{\beta_{j}\\}_{j=1}^{N_{h}^{T}}$ be sets of bases in the space $S_{h}^{0}$.
In the case of piecewise linear triangular elements, $N_{h}^{I}$ and
$N_{h}^{T}$ correspond to the number of mesh nodes in the interior of $\Omega$
and on the boundary $\Gamma_{R}$, respectively. Let
$\\{\zeta_{j}\\}_{j=1}^{N_{h}^{D}}$ represent the set of basis functions for
the space $S_{h}^{\Omega}$, where $N_{h}^{D}$ denotes the number of mesh nodes
on the boundary $\partial D$.

The discretized formulations of (5.1)–(5.2) and (5.3)–(5.4) for the IP-FEM and
BP-FEM, using linear triangular elements, can be expressed as:

$\displaystyle\big{(}\mathbf{K}-\kappa^{2}\mathbf{M}-\gamma\mathbf{K}_{J}-\mathbf{K}^{tb}\big{)}\mathbf{W}$
$\displaystyle=\mathbf{F},$ (5.5)
$\displaystyle\big{(}\mathbf{K}-\kappa^{2}\mathbf{M}-\eta\mathbf{K}_{G}-\mathbf{K}^{tb}\big{)}\mathbf{W}$
$\displaystyle=\mathbf{F}.$ (5.6)

In these equations, the penalty parameters $\gamma_{e}$ and $\eta_{e}$ are
selected as $\gamma_{e}=\gamma$ for all interior edges and $\eta_{e}=\eta$ for
all boundary edges, respectively. The unknown nodal vector $\mathbf{W}$ has a
dimension of $2N_{h}^{I}+2N_{h}^{T}+N_{h}^{D}$, given by

$\displaystyle\mathbf{W}=\begin{bmatrix}\mathbf{W}_{I}\\\ \mathbf{W}_{T}\\\
\mathbf{P}_{D}\end{bmatrix},\quad\mathbf{W}_{I}=\begin{bmatrix}\mathbf{P}_{I}\\\
\mathbf{Q}_{I}\end{bmatrix},\quad\mathbf{W}_{T}=\begin{bmatrix}\mathbf{P}_{T}\\\
\mathbf{Q}_{T}\end{bmatrix},$

where $\mathbf{P}_{I}$ and $\mathbf{Q}_{I}$ represent the values of $p$ and
$q$ at the interior nodes, $\mathbf{P}_{T}$ and $\mathbf{Q}_{T}$ denote the
values of $p$ and $q$ at the nodes on $\Gamma_{R}$, and $\mathbf{P}_{D}$
corresponds to the unknown nodal vector associated with the cavity boundary
$\partial D$.

The stiffness matrix $\mathbf{K}$ and the mass matrix $\mathbf{M}$ are given
in blockwise form as follows:

$\displaystyle\mathbf{K}=\begin{bmatrix}\mathbf{K}_{II}&\mathbf{K}_{IT}&\mathbf{K}_{ID}\\\
\mathbf{K}_{TI}&\mathbf{K}_{TT}&\mathbf{0}\\\
\mathbf{K}_{DI}&\mathbf{0}&\mathbf{0}\end{bmatrix},\quad\mathbf{M}=\begin{bmatrix}\mathbf{M}_{II}&\mathbf{M}_{IT}&\mathbf{M}_{ID}\\\
\mathbf{M}_{TI}&\mathbf{M}_{TT}&\mathbf{0}\\\
\mathbf{M}_{DI}&\mathbf{0}&2\overline{\mathbf{M}}_{DD}\end{bmatrix},$

where

$\displaystyle\mathbf{K}_{ij}=\begin{bmatrix}\overline{\mathbf{K}}_{ij}&\mathbf{0}\\\
\mathbf{0}&\overline{\mathbf{K}}_{ij}\end{bmatrix},\quad\mathbf{M}_{ij}=\begin{bmatrix}\overline{\mathbf{M}}_{ij}&\mathbf{0}\\\
\mathbf{0}&-\overline{\mathbf{M}}_{ij}\end{bmatrix},\quad i,j=I,T,$
$\displaystyle\mathbf{K}_{ID}=\begin{bmatrix}\overline{\mathbf{K}}_{ID}\\\
\overline{\mathbf{K}}_{ID}\end{bmatrix},\quad\mathbf{M}_{ID}=\begin{bmatrix}\overline{\mathbf{M}}_{ID}\\\
-\overline{\mathbf{M}}_{ID}\end{bmatrix},$
$\displaystyle\mathbf{K}_{DI}=\begin{bmatrix}\overline{\mathbf{K}}_{DI}&-\overline{\mathbf{K}}_{DI}\end{bmatrix},\quad\mathbf{M}_{DI}=\begin{bmatrix}\overline{\mathbf{M}}_{DI}&\overline{\mathbf{M}}_{DI}\end{bmatrix}.$

Specifically, the stiffness and mass matrices associated with the IP-FEM and
the BP-FEM are given in Table 1, where we have the relationships
$\overline{\mathbf{K}}_{DI}=(\overline{\mathbf{K}}_{ID})^{\top}$,
$\overline{\mathbf{K}}_{TI}=(\overline{\mathbf{K}}_{IT})^{\top}$,
$\overline{\mathbf{M}}_{DI}=(\overline{\mathbf{M}}_{ID})^{\top}$, and
$\overline{\mathbf{M}}_{TI}=(\overline{\mathbf{M}}_{IT})^{\top}$. Here $I,T$,
and $D$ stand for the interior node in $\Omega$, the boundary node on
$\Gamma_{R}$, and the boundary node on $\partial D$, respectively.

Table 1. The stiffness and mass matrices for the IP-FEM and BP-FEM. Matrix |  Dimension |  Matrix entries
---|---|---
$\overline{\mathbf{K}}_{II},\overline{\mathbf{M}}_{II}$ |  $N_{h}^{I}\times N_{h}^{I}$ |  $(\overline{\mathbf{K}}_{II})_{j,l}=\displaystyle\int_{\Omega}\nabla\alpha_{j}\cdot\nabla\alpha_{l}d\boldsymbol{x},\quad(\overline{\mathbf{M}}_{II})_{j,l}=\displaystyle\int_{\Omega}\alpha_{j}\alpha_{l}d\boldsymbol{x}$
$\overline{\mathbf{K}}_{IT},\overline{\mathbf{M}}_{IT}$ |  $N_{h}^{I}\times N_{h}^{T}$ |  $(\overline{\mathbf{K}}_{IT})_{j,l}=\displaystyle\int_{\Omega}\nabla\alpha_{j}\cdot\nabla\beta_{l}d\boldsymbol{x},\quad(\overline{\mathbf{M}}_{IT})_{j,l}=\displaystyle\int_{\Omega}\alpha_{j}\beta_{l}d\boldsymbol{x}$
$\overline{\mathbf{K}}_{ID},\overline{\mathbf{M}}_{ID}$ |  $N_{h}^{I}\times N_{h}^{D}$ |  $(\overline{\mathbf{K}}_{ID})_{j,l}=\displaystyle\int_{\Omega}\nabla\alpha_{j}\cdot\nabla\zeta_{l}d\boldsymbol{x},\quad(\overline{\mathbf{M}}_{ID})_{j,l}=\displaystyle\int_{\Omega}\alpha_{j}\zeta_{l}d\boldsymbol{x}$
$\overline{\mathbf{K}}_{TT},\overline{\mathbf{M}}_{TT}$ |  $N_{h}^{T}\times N_{h}^{T}$ |  $(\overline{\mathbf{K}}_{TT})_{j,l}=\displaystyle\int_{\Omega}\nabla\beta_{j}\cdot\nabla\beta_{l}d\boldsymbol{x},\quad(\overline{\mathbf{M}}_{TT})_{j,l}=\displaystyle\int_{\Omega}\beta_{j}\beta_{l}d\boldsymbol{x}$
$\overline{\mathbf{M}}_{DD}$ |  $N_{h}^{D}\times N_{h}^{D}$ |  $(\overline{\mathbf{M}}_{DD})_{j,l}=\displaystyle\int_{\Omega}\zeta_{j}\zeta_{l}d\boldsymbol{x}$

The matrix $\mathbf{K}^{tb}$ is associated with the TBC and is given by

$\displaystyle\mathbf{K}^{tb}=\begin{bmatrix}\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{K}_{L}^{tb}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}\end{bmatrix},$

where the matrix $\mathbf{K}_{L}^{tb}$ can be given by

$\displaystyle\mathbf{K}_{L}^{tb}=\sum_{|n|\leq
N}\begin{bmatrix}a_{n}\mathbf{K}_{n}^{tb}&0\\\
0&b_{n}\mathbf{K}_{n}^{tb}\end{bmatrix}.$

Here the truncation parameter $N$ is a positive integer, $a_{n}=h_{n}(\kappa
R)/2\pi$, $b_{n}=k_{n}(\kappa R)/2\pi$, and $\mathbf{K}_{n}^{tb}$ is evaluated
as follows:

$\displaystyle\mathbf{K}_{n}^{tb}=\left(\int_{0}^{2\pi}\mathbf{N}^{tb}e^{{\rm
i}n\theta}d\theta\right)_{N_{h}^{T}\times
1}\left(\int_{0}^{2\pi}(\mathbf{N}^{tb})^{T}e^{-{\rm
i}n\theta^{\prime}}d\theta^{\prime}\right)_{1\times N_{h}^{T}},$

where $\mathbf{N}^{tb}$ is a vector consisting of base functions
$\\{\beta_{j}\\}_{j=1}^{N_{h}^{T}}$ on the boundary $\Gamma_{R}$. The matrix
$\overline{\mathbf{K}}_{J}$ associated with the interior penalty term
$J(\phi,\psi)$ is given by

$\overline{\mathbf{K}}_{J}=\begin{bmatrix}\overline{\mathbf{J}}_{II}&\overline{\mathbf{J}}_{IT}&\overline{\mathbf{J}}_{ID}\\\
\overline{\mathbf{J}}_{TI}&\overline{\mathbf{J}}_{TT}&0\\\
\overline{\mathbf{J}}_{DI}&0&\overline{\mathbf{J}}_{DD}\end{bmatrix}=\sum_{j=\partial
K\cap\partial K^{\prime}\in\mathcal{C}_{h}}h_{j}^{2}\mathbf{k}_{j},$

where $\overline{\mathbf{J}}_{TI}=(\overline{\mathbf{J}}_{IT})^{\top}$,
$\overline{\mathbf{J}}_{DI}=(\overline{\mathbf{J}}_{ID})^{\top}$,
$\mathbf{k}_{j}=\mathbf{g}_{j}\mathbf{g}_{j}^{\top}$ with $\mathbf{g}_{j}$
being the discretized vector associated with the jump $[\partial_{\nu}\phi]$.
Specifically, the jumps $[\partial_{\nu}\phi]$ and $[\partial_{\nu}\psi]$ on
the interior edge $j=\partial K\cap\partial K^{\prime}$ in the interior
penalty term $J(\phi,\psi)$ can be written as

$[\partial_{\nu}\phi]=\partial_{\nu}^{K}\phi^{K}+\partial_{\nu}^{K^{\prime}}\phi^{K^{\prime}}=\mathbf{g}_{j}^{\top}\boldsymbol{\phi},\quad[\partial_{\nu}\psi]=\partial_{\nu}^{K}\psi^{K}+\partial_{\nu}^{K^{\prime}}\psi^{K^{\prime}}=\mathbf{g}_{j}^{\top}\boldsymbol{\psi},$

where $\boldsymbol{\phi}$ and $\boldsymbol{\psi}$ are vectors composed of the
function values of $\phi$ and $\psi$ at all nodes in the domain $\Omega$,
respectively. Additionally, the normal direction in $\partial_{\nu}^{K}$ is
opposite to that in $\partial_{\nu}^{K^{\prime}}$.

Therefore, the interior penalty stiffness matrix $\mathbf{K}_{J}$ can be
expressed as

$\displaystyle\mathbf{K}_{J}=\begin{bmatrix}\mathbf{J}_{II}&\mathbf{J}_{IT}&\mathbf{J}_{ID}\\\
\mathbf{J}_{TI}&\mathbf{J}_{TT}&\mathbf{0}\\\
\mathbf{J}_{DI}&\mathbf{0}&2\overline{\mathbf{J}}_{DD}\end{bmatrix},$

where

$\displaystyle\mathbf{J}_{ij}=\begin{bmatrix}\overline{\mathbf{J}}_{ij}&\mathbf{0}\\\
\mathbf{0}&-\overline{\mathbf{J}}_{ij}\end{bmatrix},i,j=I,T,\quad\mathbf{J}_{ID}=\begin{bmatrix}\overline{\mathbf{J}}_{ID}\\\
-\overline{\mathbf{J}}_{ID}\end{bmatrix},\quad\mathbf{J}_{DI}=\begin{bmatrix}\overline{\mathbf{J}}_{DI}&\overline{\mathbf{J}}_{DI}\end{bmatrix}.$

Let us assume that the boundary $\partial D$ of the cavity is divided into $K$
segments $\Gamma_{k}$, where $k=1,\cdots,K$, in the mesh $\mathcal{M}_{h}$.
The boundary penalty stiffness matrix $\mathbf{K}_{G}$ is associated with the
boundary penalty term $G(p_{h}^{D},\varphi_{h})$ and can be obtained by
mapping $\mathbf{K}_{G}^{L}$ from local to global numbering. The matrix
$\mathbf{K}_{G}^{L}$ can be evaluated as follows:

$\displaystyle\mathbf{K}_{G}^{L}=\sum_{k=1}^{K}\begin{bmatrix}1&-1\\\
-1&1\end{bmatrix}.$

##  6\. Numerical experiments

In this section, we present numerical results obtained using the IP-FEM and
BP-FEM for three examples: a circular-shaped cavity, an ellipse-shaped cavity,
and a kite-shaped cavity. In the experiments, we investigate the out-of-plane
displacement of the scattered field $v$ and its bending moment
$w=\kappa^{-2}\Delta v$ by solving the boundary value problem (4.5) and using
the relationships $(p^{s},q^{s})$ and $(v,\Delta v)$. The relative errors in
the $\emph{L}^{2}$ norm and the $\emph{H}^{1}$ semi-norm are employed to
assess the numerical solutions. For the circular-shaped cavity, we compare the
results against the analytic solution, while for the ellipse-shaped and kite-
shaped cavities, we use reference solutions, i.e., the numerical solutions
obtained with fine meshes. The relative errors of the $\emph{L}^{2}$ norm and
the $\emph{H}^{1}$ semi-norm of any function $\phi$ are defined as follows:

$\mathrm{E}_{L^{2}}=\frac{\|\phi^{e}-\phi^{n}\|_{0,\Omega}}{\|\phi^{e}\|_{0,\Omega}},\quad\mathrm{E}_{H^{1}}=\frac{\|\nabla\phi^{e}-\nabla\phi^{n}\|_{0,\Omega}}{\|\nabla\phi^{e}\|_{0,\Omega}},$

where $\phi^{e}$ and $\phi^{n}$ represent the analytical or reference solution
and the numerical solution, respectively.

###  6.1. A circular-shaped cavity

Consider a circular-shaped cavity $D=B_{\hat{R}}$, which is illuminated by a
plane wave

$u^{\rm inc}(\boldsymbol{x})=e^{{\rm
i}\kappa\boldsymbol{x}\cdot\boldsymbol{d}},$

where $\kappa>0$ is the wavenumber and
$\boldsymbol{d}=(\cos\alpha,\sin\alpha)$ is the incident direction with
$\alpha$ being the incident angle. The parameter equation of the circular-
shaped cavity with radius $\hat{R}$ is

$\displaystyle x(\theta)=\hat{R}\cos\theta,\quad
y(\theta)=\hat{R}\sin\theta,\quad\theta\in[0,2\pi).$

####  6.1.1. The analytical solution

The Helmholtz and modified Helmholtz wave components $v_{H}$ and $v_{M}$ of
the out-of-plane displacement of the scattered field $v$ satisfy the coupled
boundary value problem

$\left\\{\begin{aligned} &\Delta v_{H}+\kappa^{2}v_{H}=0,\quad\Delta
v_{M}-\kappa^{2}v_{M}=0&\text{in}~{}\mathbb{R}^{2}\setminus\overline{B_{\hat{R}}},\\\
&v_{H}+v_{M}=f(\theta),\quad\partial_{r}v_{H}+\partial_{r}v_{M}=g(\theta)&\text{on}~{}\partial
B_{\hat{R}},\\\ \end{aligned}\right.$ (6.1)

where $f(\theta)=-u^{\rm inc}$ and $g(\theta)=-\partial_{r}u^{\rm inc}$. The
analytical solution of (6.1) has the Fourier series expansion in polar
coordinates:

$v_{H}(r,\theta)=\sum_{n\in\mathbb{Z}}\frac{H_{n}^{(1)}(\kappa
r)}{H_{n}^{(1)}(\kappa\hat{R})}v_{H}^{(n)}(\hat{R})e^{{\rm i}n\theta},\quad
v_{M}(r,\theta)=\sum_{n\in\mathbb{Z}}\frac{K_{n}(\kappa
r)}{K_{n}(\kappa\hat{R})}v_{M}^{(n)}(\hat{R})e^{{\rm i}n\theta},$ (6.2)

where the Fourier coefficients $v_{H}^{(n)}(\hat{R})$ and
$v_{M}^{(n)}(\hat{R})$ are given by

$v_{H}^{(n)}(\hat{R})=\frac{1}{2\pi}\int_{0}^{2\pi}v_{H}(\hat{R},\theta)e^{-{\rm
i}n\theta}d\theta,\quad
v_{M}^{(n)}(\hat{R})=\frac{1}{2\pi}\int_{0}^{2\pi}v_{M}(\hat{R},\theta)e^{-{\rm
i}n\theta}d\theta.$

Since $f(\theta)$ and $g(\theta)$ are periodic functions with period $2\pi$,
we have

$\displaystyle f(\theta)=\sum\limits_{n\in\mathbb{Z}}f^{(n)}e^{{\rm
i}n\theta},\quad g(\theta)=\sum\limits_{n\in\mathbb{Z}}g^{(n)}e^{{\rm
i}n\theta},$ (6.3)

where the Fourier coefficients $f^{(n)}$ and $g^{(n)}$ are

$f^{(n)}=\frac{1}{2\pi}\int_{0}^{2\pi}f(\theta)e^{-{\rm
i}n\theta}d\theta,\quad g^{(n)}=\frac{1}{2\pi}\int_{0}^{2\pi}g(\theta)e^{-{\rm
i}n\theta}d\theta.$

Substituting (6.2)–(6.3) into the boundary condition on $\partial B_{\hat{R}}$
yields a linear system of algebraic equations

$\left\\{\begin{aligned}
&v_{H}^{(n)}(\hat{R})+v_{M}^{(n)}(\hat{R})=f^{(n)},\\\
&\kappa\frac{H_{n}^{(1)^{\prime}}(\kappa\hat{R})}{H_{n}^{(1)}(\kappa\hat{R})}v_{H}^{(n)}+\kappa\frac{K_{n}^{{}^{\prime}}(\kappa\hat{R})}{K_{n}(\kappa\hat{R})}v_{M}^{(n)}=g^{(n)},\end{aligned}\right.$

which has a matrix form

$\displaystyle A\begin{bmatrix}v_{H}^{(n)}\\\
v_{M}^{(n)}\end{bmatrix}=\begin{bmatrix}1&1\\\
\kappa\frac{H_{n}^{(1)^{\prime}}(\kappa\hat{R})}{H_{n}^{(1)}(\kappa\hat{R})}&\kappa\frac{K_{n}^{{}^{\prime}}(\kappa\hat{R})}{K_{n}(\kappa\hat{R})}\end{bmatrix}\begin{bmatrix}v_{H}^{(n)}\\\
v_{M}^{(n)}\end{bmatrix}=\begin{bmatrix}f^{(n)}\\\ g^{(n)}\end{bmatrix}.$
(6.4)

We can obtain the solution of (6.4) using Cramer’s rule that

$\displaystyle\left\\{\begin{aligned}
&v_{H}^{(n)}=\frac{1}{b_{n}}\left(\kappa^{-1}g^{(n)}-\frac{K_{n}^{{}^{\prime}}(\kappa\hat{R})}{K_{n}(\kappa\hat{R})}f^{(n)}\right),\\\
&v_{M}^{(n)}=\frac{1}{b_{n}}\left(\frac{H_{n}^{(1)^{\prime}}(\kappa\hat{R})}{H_{n}^{(1)}(\kappa\hat{R})}f^{(n)}-\kappa^{-1}g^{(n)}\right),\end{aligned}\right.$
(6.5)

where $b_{n}$ is the determinant of the coefficient matrix $A$ and is given by

$b_{n}=\left(\frac{H_{n}^{(1)^{\prime}}(\kappa\hat{R})}{H_{n}^{(1)}(\kappa\hat{R})}-\frac{K_{n}^{{}^{\prime}}(\kappa\hat{R})}{K_{n}(\kappa\hat{R})}\right).$

From (6.2) and (6.5), we can obtain the analytical solutions $v_{H}$ and
$v_{M}$. Then, using the following relationships, the scattered field $v$ and
its bending moment $w$ can be expressed as follows:

$v=v_{H}+v_{M},\quad w=v_{M}-v_{H}.$

######  Theorem 6.1.

The linear system (6.4) has a unique solution.

###### Proof.

It suffices to show that the coefficient matrix $A$ of (6.4) is nonsingular,
i.e., $\det(A)\neq 0$. A simple calculation gives

$\displaystyle\det(A)=\frac{\kappa
K_{n}{{}^{\prime}}(\kappa\hat{R})}{K_{n}(\kappa\hat{R})}-\frac{\kappa
H_{n}^{(1)^{\prime}}(\kappa\hat{R})}{H_{n}^{(1)}(\kappa\hat{R})}=\frac{1}{\hat{R}}\left(k_{n}(\kappa\hat{R})-h_{n}(\kappa\hat{R})\right).$

Taking the imaginary part of $\det(A)$ and using Lemma 4.1, we have

$\Im(\det(A))=-\frac{1}{\hat{R}}\Im\\{h_{n}(\kappa\hat{R})\\}=-\frac{2}{\pi\hat{R}}\frac{1}{|H_{n}^{(1)}(\kappa\hat{R})|^{2}}\neq
0,$

which implies that the coefficient matrix $A$ of (6.4) is nonsingular and
there exists a unique solution to the system of equations (6.4). ∎

In the experiments, we set $\hat{R}=0.3$ for the circular-shaped cavity and
the radius $R=0.6$ for the TBC. The incident angle $\alpha=\pi/3$ and the
wavenumber $\kappa=\pi$, corresponding to a wavelength $\lambda=2$. We choose
the DtN operator truncation number $N$ to be 15.

####  6.1.2. The influence of $\gamma$

In the IP-FEM, the penalty parameter $\gamma$ plays a crucial role. In this
subsection, we investigate the influence of the penalty parameter $\gamma$ on
the accuracy of the IP-FEM. If $\gamma$ is too large, it introduces artificial
dissipation in the numerical results. On the other hand, if $\gamma$ is too
small, we observe an oscillation behavior of the bending moment $w$ on the
cavity boundary, similar to what is seen in the regular linear FEM
($\gamma=0$), as shown in the left part of Figure 1.

Figure 1. Example 1: The bending moment $w$ on the boundary of cavity: (left)
The regular linear FEM ($\gamma=0$); (right) The IP-FEM ($\gamma=4.2\times
10^{-3}$) and the BP-FEM ($\eta=2.5\kappa\times 10^{-3}$).

Figure 2 presents the relative errors of the $L^{2}$ norm and the $H^{1}$
semi-norm of $v$ and $w$ with different values of the parameter $\gamma$ at
the mesh size $h=0.05$. The relative $L^{2}$ and $H^{1}$ errors of $v$
increase as $\gamma$ becomes larger, but they remain at levels of $10^{-3}$
and $10^{-2}$, respectively. However, the relative $L^{2}$ and $H^{1}$ errors
of $w$ first decrease and then increase with increasing $\gamma$, with the
minimum relative $L^{2}$ error of $w$ occurring at $\gamma=4.2\times 10^{-3}$.
The solution of $w$ on the cavity boundary for the IP-FEM with the optimal
parameter $\gamma=4.2\times 10^{-3}$ is shown in the right part of Figure 1.
It is evident that the boundary oscillation behavior of $w$ is mitigated
compared to the regular linear FEM, i.e., the linear FEM without any penalty
term ($\gamma=0$). Based on these observations, we conclude that there exists
a range of values for $\gamma$ that significantly improves the results for the
bending moment $w$ while maintaining good results for the displacement $v$.

Figure 2. Example 1: The relative errors are plotted against the parameter
$\gamma$ for $\kappa=\pi$. Figure 3. The optimal parameter $\gamma$ with
different wavenumbers and mesh sizes $h$ for the IP-FEM.

Figure 3 presents the optimal penalty parameter $\gamma$ at different
wavenumbers ($\kappa=0.5\pi,\pi,1.5\pi,2\pi$) and mesh sizes
($h=0.04,0.06,0.08,0.1$). The figure contains sixteen cases, and each
subfigure shows the variation of the relative $L^{2}$ error of $w$ with the
parameter $\gamma$. In this paper, relative $L^{2}$ errors of $w$ within
$1\times 10^{-3}$ (shown by two horizontal red dashed lines from the smallest
error $e_{\rm low}$ to $e_{\rm low}+1\times 10^{-3}$) are considered
acceptable, and the corresponding penalty parameters $\gamma$ are considered
acceptably optimal. From these subfigures, we observe that for
$\kappa=0.5\pi$, the optimal parameters are located at $\gamma=1.57\times
10^{-3}$ for $h=0.04,0.06,0.08,0.1$, and similar results are obtained for
other cases ($\kappa=\pi,1.5\pi,2\pi$). This implies that the optimal
parameter $\gamma$ is directly proportional to the wavenumber $\kappa$ and is
less affected by the mesh size $h$. For the discretized problem (5.5), a
suitable choice for the penalty parameter is $\gamma=\kappa\times 10^{-3}$.

####  6.1.3. The influence of $\eta$

In this subsection, we explore the influence of the penalty parameter $\eta$
in the BP-FEM. In this method, $\eta$ plays an important role as a penalty
parameter. An appropriate parameter value $\eta$ can effectively suppress the
oscillation of $w$ on the cavity boundary, as shown in the right part of
Figure 1. However, if $\eta$ is too large, it introduces artificial
dissipation in the numerical results. Conversely, if $\eta$ is too small, we
observe oscillation behavior of the bending moment $w$ on the cavity boundary,
similar to what is seen in the regular linear FEM ($\eta=0$), as shown in the
left part of Figure 1.

Figure 4 displays the relative $L^{2}$ and $H^{1}$ errors of $v$ and $w$ for
different values of the parameter $\eta$ at a mesh size of $h=0.05$. These
figures reveal that the relative $L^{2}$ and $H^{1}$ errors of $v$ increase as
$\eta$ increases, but they remain at levels of $10^{-3}$ and $10^{-2}$,
respectively. Conversely, the relative $L^{2}$ errors of $w$ first decrease
and then increase as $\eta$ increases, with the smallest error obtained for
$\eta\in(7.0\times 10^{-3},1.7\times 10^{-2})$. Based on these observations,
we conclude that there exists a range of values for $\eta$, where the results
for $w$ are significantly improved while maintaining good results for $v$.

Figure 4. Example 1: The relative errors are plotted against the parameter
$\eta$ for $\kappa=\pi$.

Figure 5 presents the optimal penalty parameter $\eta$ at different
wavenumbers ($\kappa=0.5\pi,\pi,1.5\pi,2\pi$) and mesh sizes
($h=0.04,0.06,0.08,0.1$). The figure contains sixteen cases, and each
subfigure shows the variation of the relative $L^{2}$ error of $w$ with the
parameter $\eta$. In this paper, relative $L^{2}$ errors of $w$ within
$1\times 10^{-3}$ (shown by two horizontal red dashed lines from the smallest
error $e_{\rm low}$ to $e_{\rm low}+1\times 10^{-3}$) are considered
acceptable, and the corresponding penalty parameters $\eta$ are considered
acceptably optimal. From these subfigures, we observe similar results to those
in Figure 3. This also implies that the optimal parameter $\eta$ is directly
proportional to the wavenumber $\kappa$ and is less affected by the mesh size
$h$. Consequently, for convenience, the penalty parameter $\eta$ for the
discretized problem (5.6) can be chosen as $\eta=2.5\kappa\times 10^{-3}$.

Figure 5. The optimal parameter $\eta$ with different wavenumbers and mesh
sizes $h$ in the BP-FEM.

####  6.1.4. The influence of $\kappa$

In this subsection, we consider the effects of the wavenumber on the solution
accuracy with a fixed mesh size of $h=0.05$. Figure 6 illustrates the relative
$L^{2}$ and $H^{1}$ errors of solutions $v$ and $w$ obtained using the regular
linear FEM ($\gamma=0$), the IP-FEM ($\gamma=\kappa\times 10^{-3}$), and the
BP-FEM ($\eta=2.5\kappa\times 10^{-3}$). From these figures, we observe that
the behavior of $v$ is similar for all three methods. However, for the
solution $w$, both the IP-FEM and the BP-FEM show significant improvements
compared to the regular linear FEM. Additionally, the errors in both $L^{2}$
and $H^{1}$ norms increase as the wavenumber $\kappa$ increases, regardless of
the method used.

Figure 6. Example 1: The relative errors are plotted against the wavenumber
with $\gamma=\kappa\times 10^{-3}$ and $\eta=2.5\kappa\times 10^{-3}$.

####  6.1.5. Convergence

In this subsection, we examine the convergence of the IP-FEM and the BP-FEM.
Figure 7 displays the relative errors of the $L^{2}$ norm and the $H^{1}$
semi-norm for the scattered field $v$ and its bending moment $w$ using
different methods. From these figures, we observe that the convergence rates
of the relative $L^{2}$ and $H^{1}$ errors of $v$ and $w$ for both the IP-FEM
and the BP-FEM achieve the optimal convergence order.

Figure 7. Example 1: The relative errors are plotted against the degrees of
freedom (Dofs) for the IP-FEM ($\gamma=\pi\times 10^{-3}$) and the BP-FEM
($\eta=2.5\pi\times 10^{-3}$).

###  6.2. An ellipse-shaped cavity

In this example, we study the flexural wave scattering by an elliptical cavity
with the clamped boundary. The boundary of the ellipse is described by the
following parametric equations:

$x(t)=a\cos(\theta),\quad y(t)=b\sin(\theta),$

where the major semi-axis $a=0.4$ and the minor semi-axis $b=0.2$. The
parameter $\theta$ ranges from $0$ to $2\pi$. In the experiments, the open
domain is truncated by a circle with a radius $R=0.6$, and we choose the
penalty parameters as follows: $\gamma=\kappa\times 10^{-3}$ for the IP-FEM
and $\eta=2.5\kappa\times 10^{-3}$ for the BP-FEM. All other related
parameters remain the same as in the first example. For the sake of
comparison, we obtain the reference solution using the IP-FEM with
$\gamma=\kappa\times 10^{-3}$ on a very fine mesh.

####  6.2.1. Accuracy

In this subsection, we consider the effectiveness of the IP-FEM and the BP-
FEM. The mesh size and the wavenumber are set as $h=0.05$ and $\kappa=\pi$,
respectively. Figure 8 shows the solutions $w$ obtained using the regular FEM
($\gamma=0$), the IP-FEM ($\gamma=\kappa\times 10^{-3}$), and the BP-FEM
($\eta=2.5\kappa\times 10^{-3}$) on the cavity boundary. For the results on
the entire domain, we only present the regular linear FEM and the IP-FEM, as
the BP-FEM yields similar outcomes to the IP-FEM. From these figures, we
observe that both the IP-FEM and the BP-FEM effectively suppress the
oscillations of the bending moment $w$ on the cavity boundary when compared
with the regular linear FEM.

Figure 8. Example 2: The bending moment $w$ on the boundary of cavity and the
entire domain: (left) The regular linear FEM ($\gamma=0$); (right) the IP-FEM
($\gamma=\kappa\times 10^{-3}$ and the BP-FEM ($\eta=2.5\kappa\times
10^{-3}$).

####  6.2.2. Convergence

The convergence of the IP-FEM and the BP-FEM is investigated in this
subsection. Figure 9 illustrates the convergence of the relative errors of the
$L^{2}$ norm and the $H^{1}$ semi-norm for the scattered field $v$ and its
bending moment $w$ using various methods. From these figures, it is evident
that the convergence rates of the relative $L^{2}$ and $H^{1}$ errors for $v$
and $w$ with the IP-FEM and the BP-FEM achieve optimal convergence orders.

Figure 9. Example 2: The relative errors are plotted against the degrees of
freedom (Dofs) for the elliptical cavity with $\gamma=\kappa\times 10^{-3}$
for the IP-FEM and with $\eta=2.5\kappa\times 10^{-3}$ for the BP-FEM.

###  6.3. A kite-shaped cavity

The subsection investigates the scattering by a kite-shaped cavity. The
parametric equations for the kite-shaped cavity boundary are given as follows:

$x(t)=a\cos(\theta)+b\cos(2\theta)-c,\quad y(t)=a\sin(\theta),$

where the parameters are defined as $\theta\in[0,2\pi)$, $a=0.3$, $b=0.2$, and
$c=0.1$. In the experiments, the open domain is truncated by a circle with a
radius of $R=0.6$. Figure 10 displays the truncated domain and the mesh used
for solving the kite-shaped cavity scattering problem. The remaining
parameters for this problem are the same as those used in the second example.
For the sake of comparison, the reference solution is obtained using the IP-
FEM with $\gamma=\kappa\times 10^{-3}$ on a very fine mesh.

Figure 10. Geometry of a kite-shaped cavity: (a) The truncated domain; (b) A
FEM mesh.

####  6.3.1. Accuracy

In this subsection, we present the regular linear FEM, the IP-FEM, and the BP-
FEM to validate the effectiveness of the proposed method. The mesh size and
the wavenumber are set as $h=0.05$ and $\kappa=\pi$, respectively. Figure 11
displays the solutions $w$ obtained using the regular linear FEM ($\gamma=0$),
the IP-FEM ($\gamma=\kappa\times 10^{-3}$), and the BP-FEM
($\eta=2.5\kappa\times 10^{-3}$) on the cavity boundary. Additionally, the
corresponding results of the regular linear FEM and the IP-FEM on the entire
domain are also presented in the figure. It is noted that the result of the
BP-FEM on the entire domain is similar to that of the IP-FEM. From these
figures, we observe that both the IP-FEM and BP-FEM effectively suppress the
oscillation of the bending moment on the cavity boundary compared with the
regular linear FEM.

Figure 11. Example 3: The bending moment $w$ on the boundary of cavity and the
entire domain: (left) the regular linear FEM ($\gamma=0$); (right) the IP-FEM
($\gamma=\kappa\times 10^{-3}$) and the BP-FEM ($\eta=2.5\kappa\times
10^{-3}$).

####  6.3.2. Convergence

In this subsection, we investigate the convergence of the IP-FEM and the BP-
FEM. Figure 12 illustrates the convergence of the relative errors of the
$L^{2}$ norm and the $H^{1}$ semi-norm for the scattered field $v$ and the
bending moment $w$. From these figures, we observe that the convergence rates
of $v$ and $w$ also achieve good convergence orders.

Figure 12. Example 3: The relative errors are plotted against the degrees of
freedom (Dofs) for the kite-shaped cavity with $\gamma=\kappa\times 10^{-3}$
for the IP-FEM and with $\eta=2.5\kappa\times 10^{-3}$ for the BP-FEM.

##  7\. Conclusion

In this paper, we have introduced and applied the IP-FEM and the BP-FEM to
investigate the flexural scattering by a clamped cavity in an infinite thin
plate. The proposed model utilizes the decomposition technique and the TBC
technique to transform a fourth-order problem on an unbounded domain into two
second-order equations with coupled boundary conditions and TBCs on a bounded
domain. To effectively suppress the oscillation of the bending moment on the
cavity boundary, we have incorporated the interior penalty (IP) and boundary
penalty (BP) techniques into the original variational formulation of the
decomposed problem. The resulting two new variational formulations, augmented
with penalty terms, are discretized using linear triangular elements.

To verify the effectiveness of the proposed method, we conducted a numerical
experiment involving flexural scattering by a circle-shaped cavity, for which
we obtained an analytical solution. The results of this experiment confirm
that both the IP-FEM and the BP-FEM successfully suppress the oscillations of
the bending moment on the cavity boundary, leading to a significant
improvement in the bending moment, while maintaining the accuracy of the
displacement compared to the regular linear FEM. Furthermore, we extended the
model to handle flexural scattering problems with cavities of different shapes
and compared the results with corresponding reference solutions. The numerical
results demonstrated that the convergence rates for the displacement and
bending moment achieved by the IP-FEM and the BP-FEM approach optimal
convergence orders.

As part of our future work, we aim to conduct further investigations into the
existence of decomposed problems using the variational method and delve deeper
into the related mathematical theory of the IP-FEM and the BP-FEM and
numerical calculations with complex-valued penalty parameters. These research
endeavors are expected to significantly contribute to enhancing the
understanding and applicability of our proposed numerical methods in the field
of flexural scattering and other related problems.

## References

  * [1] O. Aklouche, A. Pelat, S. Maugeais, and F. Gautier, Scattering of flexural waves by a pit of quadratic profile inserted in an infinite thin plate, J. Sound Vib., 375 (2016), 38–52.
  * [2] M. Amara and F. Dabaghi, An optimal $C^{0}$ finite element algorithm for the 2D biharmonic problem: theoretical analysis and numerical results, Numer. Math., 90 (2001), 19–46.
  * [3] E. Behrens and J. Guzmán, A mixed method for the biharmonic problem based on a system of first-order equations, SIAM J. Numer. Anal., 49 (2011), 789–817.
  * [4] J. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys., 114 (1994), 185–200.
  * [5] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, New York, 2012.
  * [6] L. Cai and S. Hambric, Multiple scattering of flexural waves on thin plates, J. Vib. Acoust., 138 (2016), 011009.
  * [7] H. Chen, Z. Zhang, and Q. Zou, A recovery based linear finite element method for 1D bi-harmonic problems, J. Sci. Comput., 68 (2016), 375–394.
  * [8] P. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM, Philadelphia, 2002.
  * [9] A. Climente, P. Gao, L. Wu, and J. Sánchez-Dehesa, Scattering of flexural waves from an N-beam resonator in a thin plate, J. Acoust. Soc. Am., 142 (2017), 3205–3215.
  * [10] A. Climente, A. Norris, and J. Sánchez-Dehesa, Scattering of flexural waves from a hole in a thin plate with an internal beam, J. Acoust. Soc. Am., 137 (2015), 293–302.
  * [11] D. Colton and R. Kress, Integral Equation Methods in Scattering Theory, SIAM, Philadelphia, 2013.
  * [12] D. Colton, P. Monk, and J. Sun, Analytical and computational methods for transmission eigenvalues, Inverse Probl., 26 (2010), 045011.
  * [13] A. Darabi, A. Zareei, M. Alam, and M. Leamy, Experimental demonstration of an ultrabroadband nonlinear cloak for flexural waves, Phys. Rev. Lett., 121 (2018), 174301.
  * [14] H. Dong and P. Li, A novel boundary integral formulation for the biharmonic wave scattering problem, arXiv:2301.10142.
  * [15] Y. Du and H. Wu, Preasymptotic error analysis of higher order FEM and CIP-FEM for Helmholtz equation with high wave number, SIAM J. Numer. Anal., 53 (2015), 782–804.
  * [16] B. Engquist and A. Majda, Absorbing boundary conditions for numerical simulation of waves, Proc. Natl. Acad. Sci. USA, 74 (1977), 1765–1766.
  * [17] D. Evans and R. Porter, Penetration of flexural waves through a periodically constrained thin elastic plate floating on water, J. Eng. Math., 58 (2007), 317–337.
  * [18] M. Farhat, S. Guenneau, and S. Enoch, Ultrabroadband elastic cloaking in thin plates, Phys. Rev. Lett., 103 (2009), 024301.
  * [19] R. Glowinski and O. Pironneau, Numerical methods for the first biharmonic equation and for the two-dimensional stokes problem, SIAM Rev., 21 (1978), 167–212.
  * [20] T. Gudi, N. Nataraj, and A. Pani, Mixed discontinuous Galerkin finite element method for the biharmonic equation, J. Sci. Comput., 37 (2008), 139–161.
  * [21] H. Guo, Z. Zhang, and Q. Zou, A $C^{0}$ linear finite element method for biharmonic problems, J. Sci. Comput., 74 (2018), 1397–1422.
  * [22] S. Haslinger, Mathematical Modelling of Flexural Waves in Structured Elastic Plates, University of Liverpool, Liverpool, 2014.
  * [23] S. Haslinger, R. Craster, A. Movchan, N. Movchan, and I. Jones, Dynamic interfacial trapping of flexural waves in structured plates, Proc. R. Soc. A, 472 (2016), 20150658.
  * [24] S. Haslinger, N. Movchan, A. Movchan, I. Jones, and R. Craster, Controlling flexural waves in semi-infinite platonic crystals with resonator-type scatterers, Q. J. Mech. Appl. Math., 70 (2017), 216–247.
  * [25] J. Hu and Z. Shi, A lower bound of the $L^{2}$ norm error estimate for the Adini element of the biharmonic equation, SIAM J. Numer. Anal., 51 (2013), 2651–2659.
  * [26] Y. Huang, H. Wei, W. Yang, and N. Yi, Recovery based finite element method for biharmonic equation in two dimensional, J. Comput. Math., 38 (2020), 84–102.
  * [27] X. Jiang, P. Li, and W. Zheng, Numerical solution of acoustic scattering by an adaptive DtN finite element method, Commun. Comput. Phys., 13 (2013), 1277–1244.
  * [28] X. Jiang, P. Li, J. Lv, and W. Zheng, An adaptive finite element method for the wave scattering with transparent boundary condition, J. Sci. Comput., 72 (2017), 936–956.
  * [29] B. Lamichhane, A finite element method for a biharmonic equation based on gradient recovery operators, BIT, 54 (2014), 469–484.
  * [30] W. Lee and J. Chen, Scattering of flexural wave in a thin plate with multiple circular inclusions by using the null-field integral equation approach, J. Sound Vib., 329 (2010), 1042–1061.
  * [31] W. Lee and J. Chen, Scattering of flexural wave in a thin plate with multiple circular holes by using the multipole Trefftz method, Int. J. Solids Struct., 47 (2010), 1118–1129.
  * [32] W. Lee and J. Chen, Scattering of flexural wave in a thin plate with multiple circular inclusions by using the multipole method, Int. J. Mech. Sci., 53 (2011), 617–627.
  * [33] F. Li, C. Hu, and W. Huang, Elastic wave scattering and dynamic stress concentrations in cylindrical shells with a circular cutout, J. Sound Vib., 259 (2003), 1209–1223.
  * [34] P. Li and X. Yuan, An adaptive finite element DtN method for the elastic wave scattering problem, Numer. Math., 150 (2022), 993–1033.
  * [35] Y. Li, W. Zheng, and X. Zhu, A CIP-FEM for high-frequency scattering problem with the truncated DtN boundary condition, CSIAM Trans. Appl. Math., 1 (2020), 530–560.
  * [36] M. Liu and W. Zhu, Nonlinear transformation-based broadband cloaking for flexural waves in elastic thin plates, J. Sound Vib., 445 (2019), 270–287.
  * [37] V. Matus and V. Emets, T-matrix method formulation applied to the study of flexural waves scattering from a through obstacle in a plate, J. Sound Vib., 329 (2010), 2843–2850.
  * [38] A. Norris and C. Vemula, Scattering of flexural waves on thin plates, J. Sound Vib., 181 (1995), 115–125.
  * [39] F. Olver, D. Lozier, R. Boisvert, and C. Clark, The NIST Handbook of Mathematical Functions, Cambridge University Press, New York, 2010.
  * [40] Z. Tang, Z. Fu, and C. Chen, A localized MAPS using polynomial basis functions for the fourth-order complex-shape plate bending problems, Arch. Appl. Mech., 90 (2020), 2241–2253.
  * [41] G. Victor, Structural Health Monitoring with Piezoelectric Wafer Active Sensors, 2nd edition, Academic Press, Oxford, 2014.
  * [42] C. Wang and F. Chang, Scattering of plate waves by a cylindrical inhomogeneity, J. Sound Vib., 282 (2005), 429–451.
  * [43] Z. Wang, T. Li, and H. Dong, Flexural wave scattering by varying-thickness annular inclusions on infinite thin plates, Int. J. Mech. Sci., 159 (2019), 406–416.
  * [44] Z. Wang and T. Li, A semi-analytical model for energy harvesting of flexural wave propagation on thin plates by piezoelectric composite beam resonators, Mech. Syst. Signal Pr., 147 (2021), 107137.
  * [45] M. Wang and J. Xu, The Morley element for fourth order elliptic equations in any dimensions, Numer. Math., 103 (2006), 155–169.
  * [46] H. Wu, Pre-asymptotic error analysis of CIP-FEM and FEM for the Helmholtz equation with high wave number. Part I: linear version, IMA J. Numer. Anal., 34 (2014), 1266–1288.