<figure><img src="image_1.png" /></figure>

HAL
open science

# Numerical approximation of the Poisson problem with small holes, using augmented finite elements and defective boundary conditions.

Muriel Boulakia, Céline Grandmont, Fabien Lespagnol, Paolo Zunino

► To cite this version:

Muriel Boulakia, Céline Grandmont, Fabien Lespagnol, Paolo Zunino. Numerical approximation of the Poisson problem with small holes, using augmented finite elements and defective boundary conditions.. 2023. hal-03501521v2

HAL Id: hal-03501521

https://inria.hal.science/hal-03501521v2

Preprint submitted on 30 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Numerical approximation of the Poisson problem with small
holes, using augmented finite elements and defective boundary
conditions.

MURIEL BOULAKIA

Université Paris-Saclay, UVSQ, CNRS,
Laboratoire de Mathématiques de Versailles, 78000, Versailles, France

CÉLINE GRANDMONT

Laboratoire Jacques-Louis Lions, UMR 7597, Sorbonne Université
Centre de Recherche Inria de Paris, Paris, France

Département de Mathématique, ULB, Bruxelles, Belgique

FABIEN LESPAGNOL

MOX, Department of Mathematics, Politecnico di Milano, Milan, Italy
Laboratoire Jacques-Louis Lions, UMR 7597, Sorbonne Université

Centre de Recherche Inria de Paris, Paris, France

AND

PAOLO ZUNINO

MOX, Department of Mathematics, Politecnico di Milano, Milan, Italy

January 30, 2023

Abstract

We consider the Poisson problem in a domain with small holes, as a template for developing efficient and accurate numerical approximation schemes for partial differential equations defined on domains with low-dimensional inclusions, such as embedded fibers. We propose a reduced model based on the projection of Dirichlet boundary constraints on a finite dimensional approximation space, obtaining in this way a Poisson problem with defective interface conditions. We analyze the existence of the solution of the reduced problem and for arbitrarily small holes we prove its convergence towards the original problem, the rate of which depends on the size of the inclusion and on the number of modes of the finite dimensional space. The numerical discretization of the reduced problem is addressed by the finite element method, using a computational mesh that does not fit to the holes in the framework of a fictitious domain approach. We propose a stable, optimally convergent and robust formulation with respect to the hole size that exploits an augmented Galerkin formulation based on the addition of suitable non-polynomial functions to the finite element approximation space. The properties of the discretization method are supported by numerical experiments.

**Keywords** : elliptic PDEs, small inclusions, asymptotic analysis, fictitious domain approach, finite element approximation

# 1 Introduction

Many engineering problems involve the interaction of slender bodies immersed into three-dimensional media. The interaction of thin structures with fluids or the modeling of fiber-reinforced materials are important examples of this wide class of applications. To solve these problems computationally, it is common to model slender structures as one-dimensional (1D) bodies, taking advantage of the slenderness to develop reduced order models and correspondingly decrease the computational cost. In particular, slender structures can be modelled as beams. In the case where one chooses not to reconstruct the coupling conditions on the 2D interface generated from the central line, the reduction of the model dimension creates difficulties at the level of the mathematical formulation. First, if the volume occupied by the three-dimensional (3D) structure is filled by the surrounding continuum, an overlap between the two regions is generated. This requires resorting to some specific computational modeling techniques,
such as the immersed boundary method or the fictitious domain approach. Second, the dimensional reduction of the slender body creates a mismatch in the interface conditions between the continuous medium and the fibers embedded into it. This is an additional difficulty that must be addressed in the mathematical formulation.

These issues have been only partially explored in the literature, but they affect solid-beam interaction or fluid-beam interaction, recently addressed in [40] and [20] respectively. Similar difficulties also appear in widely used simplified models for fluid-fiber interaction such as the slender body theory [34, 35]. Slender body theory is a one-way approach that accounts for the influence of the fluid on the fibers but not vice versa. To model the mismatch of kinematic conditions between the 3D fluid and the 1D fiber, a constant velocity profile is assumed at the interface between the two regions. For this reason, this approach does not account for the torsional stress on the fibers. In the future perspective, when it will be applied to more complex scenarios, our approach aims to overcome the previous limitation of fluid-fiber interaction.

Indeed, the main motivations of this work are rooted in the context of fluid-structure interaction involving slender bodies. In particular, we aim to develop the theoretical foundations for the analysis and the numerical approximation of a 3D-1D fluid-beam interaction model. Our approach relies on three main approximation steps: for the model, for the interface conditions and for the numerical approximation. First, we approximate the 3D slender structure (denoted as $\Omega_{\epsilon}^B$ in Figure 1) with a beam (namely $B_{\epsilon}$ in the figure), and correspondingly we fill the 3D volume of the original structure with fluid (being $\Omega_{\epsilon}^f$ the real fluid domain), resorting to a fictitious domain approach. Second, we approximate the physical interface as a generalized cylinder and, exploiting dimensionality reduction, we approximate the coupling interface conditions by means of projection on a space with finite dimension (on the transversal sections). This technique resolves the dimensional mismatch between the fluid and the structure. Finally, we propose a numerical approximation based on the finite element method. Since the classical discretization would lead to sub-optimal convergence rates, we present an augmented Galerkin approach that restores the desired approximation properties.

Before addressing all these issues at the level of the 3D-1D formulation, where the mathematical difficulties would be overwhelming, we take here an intermediate step, performing the rigorous analysis of the 2D-0D case that arises by studying the problem on a cross section of the 3D-1D one. The transition from 3D-1D to 2D-0D is illustrated in Figure 1. In the 2D case, instead of a genuine two-ways interaction, we study the obstacle problem where the interface condition is replaced by a nonhomogeneous Dirichlet condition on an internal boundary. When this problem is formulated and discretized in the usual way, for example resorting to a Lagrange multiplier defined on the whole internal boundary, we talk about the 2D-1D formulation. In our setting, we will maintain the high dimensional gap at the interface, by projecting the interface constraint on a finite-dimensional space. This dimensional reduction will be introduced thanks to a projection on the $2N + 1$ first Fourier modes of the interface condition, for an arbitrary $N \in \mathbb{N}^*$. Due to this 0D representation of the obstacle, we call this problem the 2D-0D formulation.

We also make a drastic simplification of the fluid model, stepping back to the Poisson equation. However, this first study already contains the main difficulties of more general mathematical models. The generalization of the 2D-0D approach to the Stokes problem is currently in progress with promising results. Keeping the general picture of fluid-beam interaction in mind, from now on we will address the simplified 2D-0D Poisson problem, that is basically a Poisson problem with one small hole of small diameter $\epsilon$, named $\omega_{\epsilon} = \epsilon\omega$ with $\omega = B(0, 1)$ being the unit ball. The generalization to the case of many holes is straightforward in light of the results that we will obtain.

After introducing a suitable mathematical formulation for the 2D-0D problem, we consider the convergence of the solution to the one of the full 2D-1D problem in the limit $\epsilon \to 0$. For standard boundary conditions, such asymptotic analysis can be found in many works, we refer the reader to [33, 32, 30, 36, 37, 6, 7] for a non-exhaustive review. The available results will be adapted here to the case of the dimensional mismatch at the interface. Our approximation results highlight that the reduced model proposed here is a consistent and reliable approximation of the full obstacle problem and, for this reason, it seems to be a good starting point to tackle the fluid-beam interaction problem. To be more specific, we will state that the accuracy of the approximation is of order $\epsilon^N$ (we refer to Theorem 3.5 for the exact statement), so that, for a given $\epsilon$, the number of Fourier modes can be chosen in order to reach some wanted precision.

At the discrete level, several methods such as penalty methods, Nitsche's method, Lagrange multiplier methods have been proposed for imposing boundary conditions, see for example [41] for a review and [28, 29, 3, 4, 16] for a non exhaustive list of specific examples. In our particular problem, to be consistent with the principle of a 2D-0D coupling, we have been inspired in particular by the *fictitious domain* methods [17], where a Lagrange multiplier is applied to impose the Dirichlet boundary conditions on the hole, while using a regular unfitted grid for the finite element discretization of the problem. A
<figure><img src="image_2.png" /><figcaption>Figure 1: Transition from the 3D-1D model to the 2D-0D model.</figcaption></figure>

consequence of this choice is that the mesh will not resolve the internal boundary. In this framework, the classical finite element method presents decreasing performances when $\varepsilon \to 0$, unless severe restrictions on the computational mesh step-size are satisfied. We therefore expect to find similar difficulties for our simplified 2D-0D problem. Such drawbacks of the computational method are generally called locking effects, see for example [11]. To avoid this effect on general unfitted computational meshes, it is possible to use methods that deviate from finite elements like the boundary element method [1]. Another possible strategy is to take advantage of the asymptotic development of the solution combined with an augmented Galerkin approach to obtain a convergence of arbitrary order with respect to the mesh size (see [11, 6, 10]). Here, after showing the limitations of the standard finite element method applied to our case, we propose an augmented Galerkin strategy that restores optimal convergence. More precisely, our approximation result (which is given by Theorem 4.6) states that the discretization error is of order $h^k$ where $h$ is the mesh size and $k$ is the polynomial order of the finite element space.

The case when the number of Fourier modes is given by $N = 0$ (which consists in only enforcing the continuity of the mean value on the interface) with homogeneous Dirichlet boundary conditions on the hole (namely $\varphi = 0$ in (2.1)) has been treated in [11] for the Helmholtz equations and obstacles of arbitrary shape. In addition to the generalization of the method to an arbitrary $N \in \mathbb{N}$ and to non-homogeneous boundary condition in the case of circular inclusion, an important step forward of our approach is that it allows to suitably choose the accuracy in the approximation of the Lagrange multiplier. Accordingly, besides the error estimates of the solution, we provide error estimates of the Lagrange multiplier that depend explicitly on the parameters $N$ and $h$. This study is particularly important in view of the forthcoming application of this coupling technique to 3D-1D problems modeling solid-beam or fluid-beam interaction, since the Lagrange multiplier is related to the forces that are exchanged at the interface. Furthermore, our setting has the advantage to enable any balance between the numerical approximation error and the model reduction error, giving rise to a computational framework that is extremely flexible. These features represent a significant improvement with respect to the previously published works [26, 27, 25] in the context of approximation of 3D-1D or 2D-0D interface conditions.

Overall, this work covers two main original topics, the asymptotic analysis of the reduced 2D-0D model toward the full 2D-1D one and the numerical approximation of it with a thorough analysis of the stability and convergence properties of different variants of the method. The body of the work is organized as follows. In Section 2, we address the asymptotic properties of the problem when the hole becomes arbitrarily small and we use them to motivate the derivation of the reduced 2D-0D model. In Section 3, we perform the well posedness and asymptotic analysis of the reduced model showing that it enjoys optimal properties. The analysis of the fully discrete version of the reduced 2D-0D model is addressed in Section 4, the results of which are supported by the numerical experiments of Section 5.

# 2 The Poisson problem in a domain with small holes

## 2.1 Problem set up

Throughout the paper, for a Lipschitz domain $\mathcal{D}$ in $\mathbb{R}^n$, we will use the classical notation $(\cdot, \cdot)_\mathcal{D}$ for the inner product on $L^2(\mathcal{D})$, more generally, for an Hilbert space $X$ defined on $\mathcal{D}$, $(\cdot, \cdot)_X$ denotes the inner product on $X$. For a vector space $V$ defined on $\mathcal{D}$, we denote by $V'$ its dual space and $\langle \cdot, \cdot \rangle_{V'}$ is the
pairing between $V'$ and $V$. For $1 \le p \le \infty$ and $m \in \mathbb{N}$, the standard notation $W^{p,m}(D)$ is used to denote the Sobolev space of functions on $D$ with all derivatives up to the order $m$ in $L^p(D)$. In the specific case $p=2$, we denote $W^{2,m}(D)$ by $H^m(D)$.

In all the paper, the notation $a \lesssim b$ is equivalent to $a \le Cb$ where $C$ will denote the constant of a generic upper bound assumed to be independent of the variables of the inequality and of the mesh size $h$, the size of the hole $\varepsilon$ and $N$ which characterizes the number of scalar constraints considered to approximate the Dirichlet condition on the hole. When it is necessary, the dependency on some parameters will be made precise, for example, if the constant depends on a domain $D$, we will write $\lesssim_D$ or $C(D)$.

Let $\Omega$ be a convex, polygonal domain of $\mathbb{R}^2$ containing the origin $0$. We introduce $\omega_\varepsilon$ an inclusion of size $\varepsilon$ defined for $\varepsilon > 0$ by $\omega_\varepsilon = \varepsilon\omega$ where $\omega = B(0, 1)$ is the open ball of center $0$ and radius $1$. In all the paper, we will assume from now on that $\varepsilon$ belongs to an admissible set $(0, \varepsilon_{max})$ where $\varepsilon_{max} = \sup\{\varepsilon > 0 | \bar{\omega}_\varepsilon \subset \Omega\}$. We denote the complementary of $\omega_\varepsilon$ in $\Omega$ by $\Omega_\varepsilon = \Omega \setminus \bar{\omega}_\varepsilon$.

In our study, we consider the following Poisson problem

$$ \left\{ \begin{array}{ll} -\Delta u_{\varepsilon} = f & \text{in } \Omega_{\varepsilon}, \\ u_{\varepsilon} = 0 & \text{on } \partial\Omega, \\ u_{\varepsilon} = \varphi & \text{on } \partial\omega_{\varepsilon} \end{array} \right. \qquad (2.1) $$

where $\varphi \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ and $f \in L^2(\Omega)$ are given. We assume that $f$ is such that $0 \notin \text{supp}(f)$.

**Remark 2.1.** In a fluid-solid interaction context, we would have a system of equations and the function $\varphi$ would be equal to the velocity of the structure, typically we would have $\varphi = \dot{\mathbf{X}} + \varepsilon\dot{\chi}\mathbf{e}_\theta$ where $\dot{\mathbf{X}}$ and $\dot{\chi}$ are constant on $\partial\omega_\varepsilon$ and $\mathbf{e}_\theta = (-\sin(\theta), \cos(\theta))$ in polar coordinates.

**Remark 2.2.** We consider in this paper a circular obstacle centered in $0$ but all the results can be generalized to an obstacle centered in $\mathbf{z}$ for all $\mathbf{z} \in \Omega$. This construction is also easily generalized to multiple obstacles, provided that they do not intersect each other and do not intersect with the support of $f$.

**Remark 2.3.** In our study, we consider Dirichlet homogeneous boundary conditions on the external boundary $\partial\Omega$. However, our work can be extended without difficulty to nonhomogeneous Dirichlet boundary conditions.

Let us now consider $u_0 \in H_0^1(\Omega)$ the unique solution of

$$ \left\{ \begin{array}{ll} -\Delta u_0 = f & \text{in } \Omega, \\ u_0 = 0 & \text{on } \partial\Omega. \end{array} \right. \qquad (2.2) $$

The function $u_0$ satisfies the following standard energy bound

$$ \|u_0\|_{1,\Omega} \lesssim \|f\|_{0,\Omega}. \qquad (2.3) $$

If $\varphi = 0$, then problem (2.2) represents the limit case of problem (2.1) when $\varepsilon \to 0$. The following result gives the rate of convergence of $u_\varepsilon$ towards $u_0$, namely an estimate of the difference $u_\varepsilon - u_0$ with respect to $\varepsilon$. If we assume that $\omega_\varepsilon \cap \text{supp} f = \emptyset$, then the solution $u_\varepsilon$ of the problem (2.1) satisfies

$$ \|u_\varepsilon - u_0\|_{1,\Omega_\varepsilon} \lesssim (-\ln(\varepsilon))^{-\frac{1}{2}} \|f\|_{0,\Omega}. \qquad (2.4) $$

This estimate is a classical result which can be seen as a consequence of the fact that the $H^1$-capacity of $\omega_\varepsilon$ tends to $0$ as $(-\ln(\varepsilon))^{-1}$, see for example [22] for a discussion about capacity. Thus, in the particular case where $\varphi = 0$, the solution of the problem (2.1) tends to the solution of the problem (2.2) when $\varepsilon \to 0$ and, to approximate the problem (2.1), one may consider the limit problem without inclusion. In other words, one could just ignore the presence of the inclusion. However, inequality (2.4) shows that the convergence is very slow with respect to the size of the inclusion. To give an idea, $(-\ln(\varepsilon))^{-\frac{1}{2}} \approx 0.201$ for $\varepsilon = 10^{-10}$. For this reason, in what follows, we will introduce and analyse a family of problems whose solutions can approximate $u_\varepsilon$ better than $u_0$, with an arbitrarily high accuracy when $\varepsilon \to 0$.
## 2.2 A reduced model based on defective interface conditions

To derive a family of reduced 2D-0D problems, the first idea consists in identifying the domain $\Omega_{\varepsilon}$ with the whole domain $\Omega$ and fill the small inclusion with the external medium. To do so, the solution of problem (2.1) is extended to $\Omega$ and the constraint on the boundary $\partial\omega_{\varepsilon}$ is imposed by an infinite dimensional Lagrange multiplier defined on the boundary $\partial\omega_{\varepsilon}$. By this way, problem (2.1) becomes: find $u_{\varepsilon} \in H_0^1(\Omega)$ and $\lambda_{\varepsilon} \in H^{-\frac{1}{2}}(\partial\omega_{\varepsilon})$ such that

$$ \left\{ \begin{array}{ll} (\nabla u_{\varepsilon}, \nabla v)_{\Omega} + \langle \lambda_{\varepsilon}, v \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}} = (f, v)_{\Omega}, & \forall v \in H_0^1(\Omega), \\ \langle \mu, u_{\varepsilon} \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}} = \langle \mu, \varphi \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}}, & \forall \mu \in H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}) \end{array} \right. \qquad (2.5) $$

where $\langle \cdot, \cdot \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}}$ denotes the pairing between $H^{-\frac{1}{2}}(\partial\omega_{\varepsilon})$ and $H^{\frac{1}{2}}(\partial\omega_{\varepsilon})$.

Let us remark that the Lagrangian multiplier $\lambda_{\varepsilon}$ corresponds to the jump of the normal derivatives of the solution through the interface. Indeed, following similar computations as in the proof of the Lemma 3.8 which is given in Subsection 3.2, we have:

$$ \lambda_{\varepsilon} = \nabla u_{\varepsilon,ext} \cdot \mathbf{n} - \nabla u_{\varepsilon,int} \cdot \mathbf{n} \text{ in } H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}) \qquad (2.6) $$

where $\mathbf{n}$ is the exterior normal to $\partial\omega_{\varepsilon}$, $u_{\varepsilon,ext}$ is the restriction of $u_{\varepsilon}$ to $\Omega_{\varepsilon}$ and $u_{\varepsilon,int}$ is the restriction of $u_{\varepsilon}$ to $\omega_{\varepsilon}$. In the setting of a fluid-structure interaction problem, the quantity of interest is the force exerted by the fluid on the structure which corresponds to $\nabla u_{\varepsilon,ext} \cdot \mathbf{n}$ and therefore, except in the case when $\varphi$ is a constant function, it does not match with $\lambda_{\varepsilon}$. In particular, to implement iterative methods, we need a numerical method to compute $\nabla u_{\varepsilon,ext} \cdot \mathbf{n}$ from $\lambda_{\varepsilon}$. This point will be discussed in more detail in Remark 3.7.

The second main idea of our reduced approach is to introduce an approximation of the constraint on $\partial\omega_{\varepsilon}$ thanks to a projection operator on a finite dimensional space that will approximate the trace space $H^{\frac{1}{2}}(\partial\omega_{\varepsilon})$. Since this approximation will strongly rely on Fourier analysis, before presenting the reduced problem in subsection 2.2.2, we first introduce some useful notations and results.

### 2.2.1 Notations and main results on Fourier analysis

In this subsection, we introduce notations that will be important in the following and some classical properties on Fourier analysis. We refer for example to [42] for more details on this topic. Since this theory will be used in the rescaled domain $\omega_{\varepsilon}$, we also introduce the notations and results related to this rescaling. The properties presented here are the basis of the approximation properties that will be satisfied by the reduced model.

First, taking advantage of the circular geometry of the hole boundary, we switch to polar coordinates, that is for $r \in \mathbb{R}^+$ and $\theta \in [0, 2\pi[$, if $u$ is a function defined in $\Omega$ and $\mathbf{x} = (\cos(\theta), \sin(\theta))$, we write that $u(r\mathbf{x}) = u(r, \theta)$. Here and in what follows, to simplify the writing, we have used the same notations for the function expressed in cartesian or polar coordinates. In particular, for all $\mathbf{x} \in \partial\omega$, we have $u(\mathbf{x}) = u(1, \theta)$ where $\theta \in [0, 2\pi[$ is such that $\mathbf{x} = (\cos(\theta), \sin(\theta))$.

We then consider the space $\mathcal{M}^N$ of trigonometric polynomials of degree less than or equal to $N$ given by

$$ \mathcal{M}^N = \{v \in L^2(\partial\omega) \mid v(\mathbf{x}) = a_0 + \sum_{n=1}^{N} (a_n \cos(n\theta) + b_n \sin(n\theta)), \text{ for a. e. } \mathbf{x} \in \partial\omega, \\ \text{ with } (a_n)_{0 \le n \le N} \in \mathbb{R}^{N+1}, (b_n)_{1 \le n \le N} \in \mathbb{R}^N \}. $$

The space $\mathcal{M}^N$ is introduced to approximate the space of $L^2(\partial\omega)$ and by a rescaling in $\varepsilon$ it will also be used to obtain an approximation space of $L^2(\partial\omega_{\varepsilon})$. The projection operator that we will consider in our reduced model approach is the $L^2$ projection on $\mathcal{M}^N$. Denoting it by $\Pi^N: L^2(\partial\omega) \to \mathcal{M}^N$, it is given by: for all $u \in L^2(\partial\omega)$,

$$ (\Pi^N u, v)_{\partial\omega} = (u, v)_{\partial\omega}, \forall v \in \mathcal{M}^N. \qquad (2.7) $$

We have for $u \in L^2(\partial\omega)$ and $\theta \in [0, 2\pi[$, if $\mathbf{x} = (\cos(\theta), \sin(\theta))$,

$$ (\Pi^N u)(\mathbf{x}) = \mathcal{A}^0(u) + \sum_{n=1}^{N} (\mathcal{A}^n(u) \cos(n\theta) + \mathcal{B}^n(u) \sin(n\theta)), $$
where the operators $\mathcal{A}^n : L^2(\partial\omega) \to \mathbb{R}$ and $\mathcal{B}^n : L^2(\partial\omega) \to \mathbb{R}$ associate to a function of $L^2(\partial\omega)$ its $n^{th}$ Fourier coefficients and are defined by: for all $n \in \mathbb{N}^*$, for all $u \in L^2(\partial\omega)$

$$
\begin{align}
\mathcal{A}^0(u) &= \frac{1}{2\pi} \int_0^{2\pi} u(1, \theta) \mathrm{d}\theta \nonumber \\
\mathcal{A}^n(u) &= \frac{1}{\pi} \int_0^{2\pi} u(1, \theta) \cos(n\theta) \mathrm{d}\theta, \tag{2.8} \\
\mathcal{B}^n(u) &= \frac{1}{\pi} \int_0^{2\pi} u(1, \theta) \sin(n\theta) \mathrm{d}\theta. \nonumber
\end{align}
$$

Moreover, a function $u$ can be written: for $\mathbf{x} = (\cos(\theta), \sin(\theta))$ with $\theta \in [0, 2\pi[$,

$$ u(\mathbf{x}) = a_0 + \sum_{n=1}^{\infty} (a_n \cos(n\theta) + b_n \sin(\theta)). $$

where $a_0 = \mathcal{A}^0(u)$ and, for all $n \ge 1$, $a_n = \mathcal{A}^n(u)$, $b_n = \mathcal{B}^n(u)$. This decomposition is called the Fourier series decomposition of $u$ and $(a_n \cos(n\theta), b_n \sin(n\theta))$ are the $n^{th}$ Fourier modes of $u$.

The operator $\Pi^N$ satisfies the following properties: for all $u \in L^2(\partial\omega)$

$$ \|u - \Pi^N u\|_{0, \partial\omega} \to 0 \text{ when } N \to \infty $$

and

$$ \|\Pi^N u\|_{0, \partial\omega}^2 = 2\pi \left( a_0^2 + \sum_{n=1}^{N} \left( \frac{a_n^2}{2} + \frac{b_n^2}{2} \right) \right) \le 2\pi \left( a_0^2 + \sum_{n=1}^{\infty} \left( \frac{a_n^2}{2} + \frac{b_n^2}{2} \right) \right) = \|u\|_{0, \partial\omega}^2. $$

So in particular $\Pi^N$ is continuous.

In our case, we want to approximate the space $H^{\frac{1}{2}}(\partial\omega)$. Since $H^{\frac{1}{2}}(\partial\omega) \subset L^2(\partial\omega)$, the elements of $H^{\frac{1}{2}}(\partial\omega)$ also admit a Fourier series decomposition that we exploit to introduce a more suitable auxiliary norm $\|\cdot\|_{\frac{1}{2}, \partial\omega}$ on $H^{\frac{1}{2}}(\partial\omega)$ depending on the Fourier coefficients and defined by: for all $u \in H^{\frac{1}{2}}(\partial\omega)$

$$ \overline{\|u\|}_{\frac{1}{2}, \partial\omega} = \left( a_0^2 + \sum_{n=1}^{\infty} (1+n)(a_n^2 + b_n^2) \right)^{\frac{1}{2}}, $$

where $a_0 = \mathcal{A}^0(u)$ and, for all $n \ge 1$, $a_n = \mathcal{A}^n(u)$, $b_n = \mathcal{B}^n(u)$.

This norm $\|\cdot\|_{\frac{1}{2}, \partial\omega}$ is well defined in $H^{\frac{1}{2}}(\partial\omega)$ and both norms $\|\cdot\|_{\frac{1}{2}, \partial\omega}$ and $\overline{\|\cdot\|_{\frac{1}{2}, \partial\omega}}$ are equivalent in $H^{\frac{1}{2}}(\partial\omega)$ (see [[23], Lemma 2.4.5] for example). This implies in particular that the following stability properties hold for $\Pi^N$:

$$ \overline{\|\Pi^N u\|}_{\frac{1}{2}, \partial\omega} \le \overline{\|u\|}_{\frac{1}{2}, \partial\omega} \quad \text{and} \quad \|\Pi^N u\|_{\frac{1}{2}, \partial\omega} \lesssim \|u\|_{\frac{1}{2}, \partial\omega}, \quad \forall u \in H^{\frac{1}{2}}(\partial\omega). $$

Now, since our objective is to work in the small hole $\omega_\epsilon$, we need to rescale all the notations previously introduced. To do so, $\epsilon > 0$ being fixed, we first introduce the rescaling operator $\Psi_\epsilon : L^2(\partial\omega_\epsilon) \to L^2(\partial\omega)$ defined, for all $v \in L^2(\partial\omega_\epsilon)$, by

$$ \Psi_{\epsilon}(v)(\mathbf{x}) = v(\epsilon\mathbf{x}). \tag{2.9} $$

In particular, $\Psi_\epsilon$ satisfies $\Psi_\epsilon(L^2(\partial\omega_\epsilon)) = L^2(\partial\omega)$ and $\Psi_\epsilon(H^{\frac{1}{2}}(\partial\omega_\epsilon)) = H^{\frac{1}{2}}(\partial\omega)$. Then, we set

$$ \mathcal{M}_{\epsilon}^{N} = \{v \in L^{2}(\partial\omega_{\epsilon}) \mid \Psi_{\epsilon}(v) \in \mathcal{M}^{N}\}, $$

and we define the operator $\Pi_\epsilon^N : L^2(\partial\omega_\epsilon) \to \mathcal{M}_\epsilon^N$ (resp. $\mathcal{A}_\epsilon^N : L^2(\partial\omega_\epsilon) \to \mathbb{R}$ and $\mathcal{B}_\epsilon^N : L^2(\partial\omega_\epsilon) \to \mathbb{R}$) by

$$ \Pi_{\epsilon}^{N} = \Psi_{\epsilon}^{-1} \circ \Pi^{N} \circ \Psi_{\epsilon} \tag{2.10} $$

(resp. $\mathcal{A}_\epsilon^n = \mathcal{A}^n \circ \Psi_\epsilon$ and $\mathcal{B}_\epsilon^n = \mathcal{B}^n \circ \Psi_\epsilon$). Note that $\mathcal{M}_\epsilon^N$ is equal to

$$ \mathcal{M}_{\epsilon}^{N} = \{v \in L^2(\partial\omega_{\epsilon}) \mid \Psi_{\epsilon}(v)(\mathbf{x}) = a_{\epsilon,0} + \sum_{n=1}^{N} (a_{\epsilon,n} \cos(n\theta) + b_{\epsilon,n} \sin(n\theta)), \text{ for a. e. } \mathbf{x} \in \partial\omega, \\ \text{with } (a_{\epsilon,n})_{0 \le n \le N} \in \mathbb{R}^{N+1}, (b_{\epsilon,n})_{1 \le n \le N} \in \mathbb{R}^{N}\}, $$
and for $v \in L^2(\partial\omega_\varepsilon)$, for $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial\omega$ with $\theta \in [0, 2\pi[$, we have

$$ (\Pi^N \circ \Psi_\varepsilon)(v)(\mathbf{x}) = a_{\varepsilon,0} + \sum_{n=1}^{N} (a_{\varepsilon,n} \cos(n\theta) + b_{\varepsilon,n} \sin(n\theta)), $$

with $a_{\varepsilon,0} = \mathcal{A}_\varepsilon^0(v)$ and, for all $n \in \mathbb{N}^*$, $a_{\varepsilon,n} = \mathcal{A}_\varepsilon^n(v)$, $b_{\varepsilon,n} = \mathcal{B}_\varepsilon^n(v)$.

We then consider the following norm on $L^2(\partial\omega_\varepsilon)$: for all $v \in L^2(\partial\omega_\varepsilon)$

$$ \|v\|_{0,\varepsilon} = \|\Psi_\varepsilon(v)\|_{0,\partial\omega} $$

and the associated scalar product $(\cdot, \cdot)_\varepsilon$ defined by: for all $u$ and $v \in L^2(\partial\omega_\varepsilon)$

$$ (u, v)_\varepsilon = (\Psi_\varepsilon(u), \Psi_\varepsilon(v))_{\partial\omega}. $$

In a similar way, on the space $H^{\frac{1}{2}}(\partial\omega_\varepsilon)$, we consider the rescaled norm $\|v\|_{\frac{1}{2},\varepsilon} = \|\Psi_\varepsilon(v)\|_{\frac{1}{2},\partial\omega}$. We also set $\|\cdot\|_{\frac{1}{2},\varepsilon}$ the norm on $H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ defined by $\|v\|_{\frac{1}{2},\varepsilon} = \|\Psi_\varepsilon(v)\|_{\frac{1}{2},\partial\omega}$. Let us note that for all $v \in L^2(\partial\omega_\varepsilon)$, we have

$$ \|\Pi_{\varepsilon}^{N} v\|_{0,\varepsilon}^{2} = 2\pi \left( a_{\varepsilon,0}^{2} + \sum_{n=1}^{N} \left( \frac{a_{\varepsilon,n}^{2}}{2} + \frac{b_{\varepsilon,n}^{2}}{2} \right) \right) \quad \text{and} \quad \|v\|_{0,\varepsilon}^{2} = 2\pi \left( a_{\varepsilon,0}^{2} + \sum_{n=1}^{\infty} \left( \frac{a_{\varepsilon,n}^{2}}{2} + \frac{b_{\varepsilon,n}^{2}}{2} \right) \right) \quad (2.11) $$

and if $v \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$,

$$ \overline{\|\Pi_{\varepsilon}^{N} v\|_{\frac{1}{2},\varepsilon}^{2}} = a_{\varepsilon,0}^{2} + \sum_{n=1}^{N} (1+n) (a_{\varepsilon,n}^{2} + b_{\varepsilon,n}^{2}), \quad \text{and} \quad \overline{\|v\|_{\frac{1}{2},\varepsilon}^{2}} = a_{\varepsilon,0}^{2} + \sum_{n=1}^{\infty} (1+n) (a_{\varepsilon,n}^{2} + b_{\varepsilon,n}^{2}) \quad (2.12) $$

Since the norms $\|\cdot\|_{\frac{1}{2},\partial\omega}$ and $\|\cdot\|_{\frac{1}{2},\partial\omega}$ are equivalent, the norms $\|\cdot\|_{\frac{1}{2},\varepsilon}$ and $\|\cdot\|_{\frac{1}{2},\varepsilon}$ are also equivalent. Moreover, since the constants appearing in the norm equivalence are the same as the ones appearing for $\|\cdot\|_{\frac{1}{2},\partial\omega}$ and $\|\cdot\|_{\frac{1}{2},\partial\omega}$, they are independent of $\varepsilon$. Thanks to the last remark, equations (2.11) and equations (2.12), we can deduce some stability properties on $\Pi_\varepsilon^N$ for the $L^2$ and $H^{\frac{1}{2}}$ rescaled norm:

$$ \|\Pi_{\varepsilon}^{N} u\|_{0,\varepsilon} \le \|u\|_{0,\varepsilon}, \quad \forall u \in L^2(\partial\omega_{\varepsilon}) \\ \overline{\|\Pi_{\varepsilon}^{N} u\|_{\frac{1}{2},\varepsilon}} \le \overline{\|u\|_{\frac{1}{2},\varepsilon}} \quad \text{and} \quad \|\Pi_{\varepsilon}^{N} u\|_{\frac{1}{2},\varepsilon} \lesssim \|u\|_{\frac{1}{2},\varepsilon}, \forall u \in H^{\frac{1}{2}}(\partial\omega_{\varepsilon}). $$

To conclude this part, we notice that $H^{-\frac{1}{2}}(\partial\omega_\varepsilon)$ can be equipped with the duality product $\langle \cdot, \cdot \rangle_{-\frac{1}{2},\varepsilon}$ and the norm $\|\cdot\|_{-\frac{1}{2},\varepsilon}$ defined for $\lambda \in H^{-\frac{1}{2}}(\partial\omega_\varepsilon)$ and $\mu \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ by

$$ \langle \lambda, \mu \rangle_{-\frac{1}{2},\varepsilon} = \varepsilon^{-1} \langle \lambda, \mu \rangle_{-\frac{1}{2},\partial\omega_\varepsilon} \quad \text{and} \quad \|\lambda\|_{-\frac{1}{2},\varepsilon} = \sup_{\mu \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)} \frac{\langle \lambda, \mu \rangle_{-\frac{1}{2},\varepsilon}}{\|\mu\|_{\frac{1}{2},\varepsilon}}. \quad (2.13) $$

Let us note that the duality product $\langle \cdot, \cdot \rangle_{-\frac{1}{2},\varepsilon}$ verifies in particular for $\lambda \in L^2(\partial\omega_\varepsilon)$ and $\mu \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$,

$$ \langle \lambda, \mu \rangle_{-\frac{1}{2},\varepsilon} = \varepsilon^{-1} (\lambda, \mu)_{\partial\omega_\varepsilon} = (\lambda, \mu)_\varepsilon. \quad (2.14) $$

### 2.2.2 Formulation of the 2D-0D reduced problem

We denote by $(u_\varepsilon^N, \lambda_\varepsilon^N) \in H_0^1(\Omega) \times H^{-\frac{1}{2}}(\partial\omega_\varepsilon)$ the solution of the 2D-0D reduced problem. Taking into account the hypothesis that $\varepsilon$ is small, our reduced model will be introduced following the idea to replace the boundary condition $u_\varepsilon = \varphi$ on $\partial\omega_\varepsilon$ with the finite-dimensional condition

$$ \Pi_{\varepsilon}^{N} \circ \mathcal{T}_{\partial\omega_{\varepsilon}}(u_{\varepsilon}^{N}) = \Pi_{\varepsilon}^{N}(\varphi), \quad (2.15) $$

for any $N \ge 0$, where $\mathcal{T}_{\partial\omega_\varepsilon}$ is the trace operator on $\partial\omega_\varepsilon$. In other words, instead of imposing that $u = \varphi$ on $\partial\omega_\varepsilon$, we just impose the equality of the first $(2N+1)$ Fourier coefficients of $u$ and $\varphi$. In particular, for $N=0$, the space $\mathcal{M}_\varepsilon^N$ is the space of constant functions on $\partial\omega_\varepsilon$ and this constraint is equivalent to substituting the trace constraint by a constraint on the average of $u$ on $\partial\omega_\varepsilon$.
Associated to the finite-dimensional constraint (2.15), we will look for the reduced Lagrange multiplier $\lambda_{\varepsilon}^N$ in the space $\mathcal{M}_{\varepsilon}^N$ viewed as an approximation space of $H^{-\frac{1}{2}}(\partial\omega_{\varepsilon})$.

Since $\mathcal{M}_{\varepsilon}^N \subset L^2(\partial\omega_{\varepsilon})$ and using (2.14), system (2.5) becomes: for all $v \in H_0^1(\Omega)$

$$ (\nabla u_{\varepsilon}^N, \nabla v)_{\Omega} + (\lambda_{\varepsilon}^N, \Pi_{\varepsilon}^N \circ \mathcal{T}_{\partial\omega_{\varepsilon}}(v))_{\varepsilon} = (f, v)_{\Omega} $$

and for all $\mu \in \mathcal{M}_{\varepsilon}^N$

$$ (\mu, \Pi_{\varepsilon}^N \circ \mathcal{T}_{\partial\omega_{\varepsilon}}(u_{\varepsilon}^N))_{\varepsilon} = (\mu, \varphi)_{\varepsilon}, $$

Moreover, these expressions can be simplified since, for all $v \in H_0^1(\Omega)$ and $\mu \in \mathcal{M}_{\varepsilon}^N$,

$$ (\lambda_{\varepsilon}^N, \Pi_{\varepsilon}^N \circ \mathcal{T}_{\partial\omega_{\varepsilon}}(v))_{\varepsilon} = (\lambda_{\varepsilon}^N, v)_{\varepsilon} \quad \text{and} \quad (\mu, \Pi_{\varepsilon}^N \circ \mathcal{T}_{\partial\omega_{\varepsilon}}(u_{\varepsilon}^N))_{\varepsilon} = (\mu, u_{\varepsilon}^N)_{\varepsilon}. $$

As a consequence, the 2D-0D fictitious domain formulation reads as follows: find $u_{\varepsilon}^N \in H_0^1(\Omega)$ and $\lambda_{\varepsilon}^N \in \mathcal{M}_{\varepsilon}^N$ such that

$$ \left\{ \begin{array}{ll} (\nabla u_{\varepsilon}^N, \nabla v)_{\Omega} + (\lambda_{\varepsilon}^N, v)_{\varepsilon} = (f, v)_{\Omega}, & \forall v \in H_0^1(\Omega), \\ (\mu, u_{\varepsilon}^N)_{\varepsilon} = (\mu, \varphi)_{\varepsilon}, & \forall \mu \in \mathcal{M}_{\varepsilon}^N. \end{array} \right. \qquad (2.16) $$

Eventually, notice that problem (2.16) can be formally written in the strong form: find $u_{\varepsilon}^N \in H_0^1(\Omega)$ and $\lambda_{\varepsilon}^N \in \mathcal{M}_{\varepsilon}^N$ such that

$$ \left\{ \begin{array}{ll} -\Delta u_{\varepsilon}^N + \varepsilon^{-1} \lambda_{\varepsilon}^N \delta_{\partial\omega_{\varepsilon}} = f & \text{in } \Omega, \\ \Pi_{\varepsilon}^N \circ \mathcal{T}_{\partial\omega_{\varepsilon}}(u_{\varepsilon}^N) = \varphi & \text{on } \partial\omega_{\varepsilon}. \end{array} \right. \qquad (2.17) $$

where the distribution $v\delta_{\partial\omega_{\varepsilon}}$ is such that for all $\psi \in C_0^{\infty}(\Omega)$, $\langle v\delta_{\partial\omega_{\varepsilon}}, \psi \rangle_{\Omega} = \langle v, \psi \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}}$.

Let us notice that the operator $\Pi_{\varepsilon}^N \circ \mathcal{T}_{\partial\omega_{\varepsilon}}$ can be interpreted as a restriction operator of co-dimension 2 and it resolves the mismatch about the dimension of variables at the interface, as mentioned in the introduction. The next two sections will prove that this operator allows us to define a well-posed 2D-0D problem equivalent to the 2D-1D problem at order $N$ when $\varepsilon \to 0$.

# 3 Analysis of the 2D-0D reduced Poisson problem

## 3.1 Well-posedness of the 2D-0D reduced problem

In this section, we will study the well-posedness of the problem (2.16) and establish the following theorem.

**Theorem 3.1.** Problem (2.16) admits a unique solution $(u_{\varepsilon}^N, \lambda_{\varepsilon}^N)$ in $H_0^1(\Omega) \times \mathcal{M}_{\varepsilon}^N$. Moreover, we have the following inequality:

$$ \|u_{\varepsilon}^N\|_{1,\Omega} + \|\lambda_{\varepsilon}^N\|_{-\frac{1}{2},\varepsilon} \lesssim \|\varphi\|_{\frac{1}{2},\varepsilon} + \|f\|_{0,\Omega}. \qquad (3.1) $$

We remind that $\varepsilon \in (0, \varepsilon_{max})$ where $\varepsilon_{max}$ is the fixed parameter defined by $\varepsilon_{max} = \sup\{\varepsilon > 0 | \bar{\omega}_{\varepsilon} \subset \Omega\}$.

This result states that the solution of the reduced problem continuously depends on the data and the continuity constant does not depend on $\varepsilon$ and $N$. Before proving this result, we will introduce several useful lemmas. The first lemma gives the classical Poincaré inequality.

**Lemma 3.2 (Poincaré inequality).** Let $\mathcal{D}$ be a domain in $\mathbb{R}^2$. For all $v \in H_0^1(\mathcal{D})$, we have

$$ \|v\|_{0,\mathcal{D}} \lesssim_{\mathcal{D}} \|\nabla v\|_{0,\mathcal{D}}. $$

The second lemma gives an estimate of the trace on the boundary of the hole $\partial\omega_{\varepsilon}$ of a $H^1$ function defined in $\Omega_{\varepsilon}$.

**Lemma 3.3.** For all $v \in H^1(\Omega)$,

$$ \|v\|_{\frac{1}{2},\varepsilon} \lesssim \|v\|_{1,\Omega_{\varepsilon}}. $$
Proof. According to [[31], Section 4.1.3], if we consider the norm defined for $u \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ by

$$ \langle u \rangle_{\frac{1}{2}, \partial\omega_{\varepsilon}} = \inf_{v \in H^1(\Omega_{\varepsilon}), T_{\partial\omega_{\varepsilon}} v = u} \|v\|_{1, \Omega_{\varepsilon}}, \quad (3.2) $$

then the norm $\langle \cdot \rangle_{\frac{1}{2}, \partial\omega_\varepsilon}$ and the norm $\| \cdot \|_{\frac{1}{2}, \varepsilon}$ are equivalent independently of $\varepsilon$. By definition of the norm $\langle \cdot \rangle_{\frac{1}{2}, \partial\omega_\varepsilon}$ given in (3.2), we then have, for all $v \in H^1(\Omega)$,

$$ \|v\|_{\frac{1}{2}, \varepsilon} \lesssim \langle v \rangle_{\frac{1}{2}, \partial\omega_\varepsilon} \lesssim \|v\|_{1, \Omega_\varepsilon}. $$

□

In the next lemma, we introduce a continuous lifting in the whole domain $\Omega$ of functions given in $H^{\frac{1}{2}}(\partial\omega_\varepsilon)$.

**Lemma 3.4.** Let $z \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ be given. We define $v_\varepsilon^z \in H^1(\Omega)$ by

$$ \left\{ \begin{array}{ll} -\Delta v_\varepsilon^z = 0 & \text{in } \Omega_\varepsilon, \\ -\Delta v_\varepsilon^z = 0 & \text{in } \omega_\varepsilon, \\ v_\varepsilon^z = z & \text{on } \partial\omega_\varepsilon, \\ v_\varepsilon^z = 0 & \text{on } \partial\Omega. \end{array} \right. \qquad (3.3) $$

Then,

$$ \|v_\varepsilon^z\|_{1, \Omega} \lesssim \|z\|_{\frac{1}{2}, \varepsilon}. \qquad (3.4) $$

**Proof.** According to Lemma A.3,

$$ \|v_\varepsilon^z\|_{1, \Omega_\varepsilon} \lesssim \|z\|_{\frac{1}{2}, \varepsilon}. $$

Next, for the estimate in $\omega_\varepsilon$, we will give an explicit formula of $v_\varepsilon^z$. Since $z \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$, $z$ can be written in the following way: for $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial\omega$ with $\theta \in [0, 2\pi[$, we have

$$ z(\varepsilon\mathbf{x}) = a_{\varepsilon,0} + \sum_{n=1}^{\infty} (a_{\varepsilon,n} \cos(n\theta) + b_{\varepsilon,n} \sin(n\theta)). $$

Since $v_\varepsilon$ is harmonic in $\omega_\varepsilon$, by the method of separation of variables, we obtain that, for $0 < r < \varepsilon$ and $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial\omega$ with $\theta \in [0, 2\pi[$,

$$ v_{\varepsilon}^{z}(r\mathbf{x}) = a_{\varepsilon,0} + \sum_{n=1}^{\infty} \left(\frac{r}{\varepsilon}\right)^{n} (a_{\varepsilon,n} \cos(n\theta) + b_{\varepsilon,n} \sin(n\theta)). $$

The derivation of this expression can be found for example in [24]. We then deduce by orthogonality of the basis functions $\sin(n\theta)$ and $\cos(n\theta)$ that

$$ \|v_{\varepsilon}^{z}\|_{0, \omega_{\varepsilon}}^{2} = \pi \varepsilon^{2} a_{\varepsilon,0}^{2} + \pi \sum_{n=1}^{\infty} \varepsilon^{-2n} (a_{\varepsilon,n}^{2} + b_{\varepsilon,n}^{2}) \int_{0}^{\varepsilon} r^{2n+1} dr = \pi \varepsilon^{2} a_{\varepsilon,0}^{2} + \pi \sum_{n=1}^{\infty} \frac{1}{2n+2} (a_{\varepsilon,n}^{2} + b_{\varepsilon,n}^{2}), $$

and

$$ \left\| \frac{1}{r} \frac{\partial v_{\varepsilon}^{z}}{\partial \theta} \right\|_{0, \omega_{\varepsilon}}^{2} = \left\| \frac{\partial v_{\varepsilon}^{z}}{\partial r} \right\|_{0, \omega_{\varepsilon}}^{2} = \pi \sum_{n=1}^{\infty} \varepsilon^{-2n} n^2 (a_{\varepsilon,n}^2 + b_{\varepsilon,n}^2) \int_{0}^{\varepsilon} r^{2n-1} dr = \pi \sum_{n=1}^{\infty} \frac{n}{2} (a_{\varepsilon,n}^2 + b_{\varepsilon,n}^2). $$

So, we get that

$$ \|v_{\varepsilon}^{z}\|_{1, \omega_{\varepsilon}}^{2} = \|v_{\varepsilon}^{z}\|_{0, \omega_{\varepsilon}}^{2} + \|\nabla v_{\varepsilon}\|_{0, \omega_{\varepsilon}}^{2} \lesssim \left( a_{\varepsilon,0}^{2} + \sum_{n=1}^{\infty} (1+n) (a_{\varepsilon,n}^{2} + b_{\varepsilon,n}^{2}) \right) = C \|\overline{z}\|_{\frac{1}{2}, \varepsilon}^{2}. $$

To conclude, we use that, as stated in subsection 2.2.1, the norms $\| \cdot \|_{\frac{1}{2}, \varepsilon}$ and $\| \cdot \|_{\frac{1}{2}, \varepsilon}$ are equivalent independently of $\varepsilon$. □

We are now in position to prove the main result of this subsection.
Proof of Theorem 3.1. First, we notice that problem (2.16) can be written under the form: look for $(u_{\varepsilon}^N, \lambda_{\varepsilon}^N) \in H_0^1(\Omega) \times \mathcal{M}_{\varepsilon}^N$ such that

$$ \begin{cases} a(u_{\varepsilon}^N, v) + b(\lambda_{\varepsilon}^N, v) = (f, v)_{\Omega} & \forall v \in H_0^1(\Omega), \\ b(\mu, u_{\varepsilon}^N) = c(\mu) & \forall \mu \in \mathcal{M}_{\varepsilon}^N, \end{cases} \quad (3.5) $$

where the bilinear forms $a: H_0^1(\Omega) \times H_0^1(\Omega) \to \mathbb{R}$, $b: \mathcal{M}_{\varepsilon}^N \times H_0^1(\Omega) \to \mathbb{R}$ and the linear form $c: \mathcal{M}_{\varepsilon}^N \to \mathbb{R}$ are given by

$$ a(u, v) = (\nabla u, \nabla v)_{\Omega}, \quad b(\mu, u) = (\mu, u)_{\varepsilon}, \quad c(\mu) = (\mu, \varphi)_{\varepsilon}. \quad (3.6) $$

Let us mention that, to lighten the writing, we have omitted in the notations of $b$ and $c$ their dependence in $N$ and $\varepsilon$.

To prove our result, we will apply to problem (3.5) the version of the BNB theorem for saddle point problems (we refer to [8] or to [[13], Theorem 2.34]).

We notice that, thanks to the Poincaré inequality recalled in Lemma 3.2 and to Cauchy-Schwarz inequality, the bilinear form $a$ is continuous and coercive with constants independent of $\varepsilon$ and $N$.

Next, we have, for all $\mu \in \mathcal{M}_{\varepsilon}^N$

$$ |c(\mu)| \lesssim \|\mu\|_{-\frac{1}{2}, \varepsilon} \|\varphi\|_{\frac{1}{2}, \varepsilon}. $$

So $c: \mathcal{M}_{\varepsilon}^N \to \mathbb{R}$ is continuous and the continuity constant is bounded by $C\|\varphi\|_{\frac{1}{2}, \varepsilon}$. Moreover, for all $u \in H_0^1(\Omega)$ and $\mu \in \mathcal{M}_{\varepsilon}^N$,

$$ |b(\mu, u)| = |(\mu, u)_{\varepsilon}| \le \|\mu\|_{-\frac{1}{2}, \varepsilon} \|u\|_{\frac{1}{2}, \varepsilon}. $$

According to Lemma 3.3, we have

$$ \|u\|_{\frac{1}{2}, \varepsilon} \lesssim \|u\|_{1, \Omega}. $$

So we get that $b$ is continuous with a continuity constant independent of $N$ and $\varepsilon$. Let us now prove the inf-sup condition which consists in showing that there exists $\beta > 0$ such that, for all $\mu \in \mathcal{M}_{\varepsilon}^N$,

$$ \sup_{v \in H_0^1(\Omega)} \frac{(\mu, v)_{\varepsilon}}{\|v\|_{1, \Omega}} \ge \beta \|\mu\|_{-\frac{1}{2}, \varepsilon}. \quad (3.7) $$

According to Lemma 3.4, for all $z \in H_{\frac{1}{2}}(\partial\omega_{\varepsilon})$, there exists $v_{\varepsilon}^z \in H_0^1(\Omega)$ such that $v_{\varepsilon}^z = z$ on $\partial\omega_{\varepsilon}$ and

$$ \|v_{\varepsilon}^z\|_{1, \Omega} \le C \|z\|_{\frac{1}{2}, \varepsilon}. $$

We deduce that, for all $\mu \in \mathcal{M}_{\varepsilon}^N$,

$$ \|\mu\|_{-\frac{1}{2}, \varepsilon} = \sup_{z \in H_{\frac{1}{2}}(\partial\omega_{\varepsilon})} \frac{(\mu, z)_{\varepsilon}}{\|z\|_{\frac{1}{2}, \varepsilon}} \le C \sup_{v \in H_0^1(\Omega)} \frac{(\mu, v)_{\varepsilon}}{\|v\|_{1, \Omega}} $$

with $C$ independent of $N$ and $\varepsilon$. Thus we get (3.7) with $\beta$ independent of $\varepsilon$ and $N$. By applying BNB theorem, we then deduce that (3.5) admits a unique solution $(u_{\varepsilon}^N, \lambda_{\varepsilon}^N) \in H_0^1(\Omega) \times \mathcal{M}_{\varepsilon}^N$ that satisfies estimate (3.1). $\square$

## 3.2 Convergence of the reduced problem towards the full problem

In the previous results, the domain of validity for the $\varepsilon$ parameter only depended on the $\Omega$ domain. From now on, some results will also depend on an additional parameter $\Upsilon$ defined by $\Upsilon = \min(d(\{0\}, \text{supp}(f)), d(\{0\}, \partial\Omega))$ where $\{0\}$ is the center of $\omega_{\varepsilon}$ and $d(\cdot, \cdot)$ is the Euclidean distance of two sets. Let us note that we have in particular $\Upsilon \le \varepsilon_{\text{max}}$. In this subsection, we will study the convergence when $\varepsilon$ tends to 0 of the solution of the reduced problem (2.16) to the solution of the full problem (2.5) and prove the following theorem

**Theorem 3.5.** We assume that $\varphi$ belongs to $\mathcal{M}_{\varepsilon}^N$. Let $(u_{\varepsilon}, \lambda_{\varepsilon}) \in H_0^1(\Omega) \times H_{-\frac{1}{2}}(\partial\omega_{\varepsilon})$ be the solution of (2.5) and let $(u_{\varepsilon}^N, \lambda_{\varepsilon}^N) \in H_0^1(\Omega) \times \mathcal{M}_{\varepsilon}^N$ be the solution of (2.16). Then, for $\varepsilon \in (0, \Upsilon)$, we have

$$ \|u_{\varepsilon} - u_{\varepsilon}^N\|_{1, \Omega} \lesssim \left(\frac{\varepsilon}{\Upsilon}\right)^{N+1} \left(\|\varphi\|_{\frac{1}{2}, \varepsilon} + \|f\|_{0, \Omega}\right) \quad (3.8) $$

and

$$ \|\lambda_{\varepsilon} - \varepsilon^{-1}\lambda_{\varepsilon}^N\|_{-\frac{1}{2}, \varepsilon} \lesssim \left(\frac{\varepsilon}{\Upsilon}\right)^N \left(\|\varphi\|_{\frac{1}{2}, \varepsilon} + \|f\|_{0, \Omega}\right). \quad (3.9) $$
Thanks to this result, we see that the convergence in $\varepsilon$ of $u_\varepsilon - u_\varepsilon^N$ behaves like $\varepsilon^{N+1}$ where $N$ is the degree of the trigonometric polynomials of the approximation space $\mathcal{M}_\varepsilon^N$. In particular, this result can be used to adjust $N$ in order to reach a certain accuracy, $\varepsilon$ being fixed. Notice that, for $N = 0$, inequality (3.9) does not give the convergence of $\lambda_\varepsilon - \varepsilon^{-1}\lambda_\varepsilon^N$ in the $H^{-\frac{1}{2}}$ rescaled norm $\|\cdot\|_{-\frac{1}{2},\varepsilon}$. However, if we use the definition (2.13) of the $H^{-\frac{1}{2}}$ rescaled norm, we can prove that

$$ \|\mu\|_{-\frac{1}{2},\partial\omega_\varepsilon} \lesssim \varepsilon^{\frac{1}{2}} \|\mu\|_{-\frac{1}{2},\varepsilon}, $$

so, if we consider the standard norm $\|\cdot\|_{-\frac{1}{2},\partial\omega_\varepsilon}$, we obtain

$$ \|\lambda_\varepsilon - \varepsilon^{-1}\lambda_\varepsilon^N\|_{-\frac{1}{2},\partial\omega_\varepsilon} \lesssim \left(\frac{\varepsilon}{\Upsilon}\right)^{N+\frac{1}{2}} \left(\|\varphi\|_{\frac{1}{2},\varepsilon} + \|f\|_{0,\Omega}\right). $$

This parameter $\Upsilon$ plays an important role in the extension of the method to several obstacles. In that case, $\Upsilon$ also depends on the distances between inclusions and so we see that, when the inclusions get closer together, the estimate of $u_\varepsilon - u_\varepsilon^N$ deteriorates. The advantage of our approach is that, in that case, the loss of precision can be compensated by increasing $N$ the number of moments.

**Remark 3.6.** To simplify the presentation, in this subsection, we assume that $\varphi$ belongs to $\mathcal{M}_\varepsilon^N$. This assumption is satisfied for example for $N \ge 1$ if $\varphi$ corresponds to the coordinates of the velocity of a rigid solid in 2D described in Remark 2.1. However, the convergence properties given by Theorem 3.5 remain valid if

$$ \|\varphi - \Pi_\varepsilon^N \varphi\|_{\frac{1}{2},\varepsilon} \lesssim \varepsilon^{N+1}. $$

**Remark 3.7.** As pointed out in the beginning of Section 2.2, for the numerical implementation of iterative methods, the quantity of interest which corresponds to the force exerted by the fluid is $\nabla u_{\varepsilon,ext}^N \cdot \mathbf{n}$ which does not coincide with the Lagrange multiplier $\lambda_\varepsilon^N$. Nevertheless we can notice that this quantity can be computed from $\lambda_\varepsilon^N$ using the next Lemma 3.8 and the explicit expression of $\nabla u_{\varepsilon,int}^N \cdot \mathbf{n}$ provided by formula (3.18). More precisely, we have in $H^{-\frac{1}{2}}(\partial\omega_\varepsilon)$

$$ \nabla u_{\varepsilon,ext}^N \cdot \mathbf{n} = \frac{1}{\varepsilon}\lambda_\varepsilon^N + \frac{1}{\varepsilon}\sum_{n=1}^N n(\mathcal{A}_\varepsilon^n(\varphi)\cos(n\theta) + \mathcal{B}_\varepsilon^n(\varphi)\sin(n\theta)) + R_N \quad (3.10) $$

where $R_N$ is defined by $R_N = \sum_{n \ge N+1} n\varepsilon^{n-1}(\cos(n\theta) + B_{\varepsilon,n}\sin(n\theta))$, with the coefficients $A_{\varepsilon,n}$ and $B_{\varepsilon,n}$ such that

$$ \sum_{n \ge N+1} (1+n)(|A_{\varepsilon,n}|^2 + |B_{\varepsilon,n}|^2) \lesssim \|\varphi\|_{\frac{1}{2},\partial\omega_\varepsilon}^2 + \|f\|_{0,\Omega}^2. $$

This inequality implies that

$$ \|R_N\|_{-\frac{1}{2},\varepsilon} \lesssim \varepsilon^N (\|\varphi\|_{\frac{1}{2},\partial\omega_\varepsilon} + \|f\|_{0,\Omega}). $$

Thus, to compute $\nabla u_{\varepsilon,ext}^N \cdot \mathbf{n}$ with an accuracy of $\varepsilon^N$ it is enough to compute the first two terms of (3.10). Moreover, according to (2.6) and the definition of $u_{\varepsilon,int}$, we have

$$ \nabla u_{\varepsilon,ext}^N \cdot \mathbf{n} - \nabla u_{\varepsilon,ext} \cdot \mathbf{n} = \frac{1}{\varepsilon}\lambda_\varepsilon^N - \lambda_\varepsilon + R_N $$

Therefore, $\nabla u_{\varepsilon,ext}^N \cdot \mathbf{n}$ gives an approximation of $\nabla u_{\varepsilon,ext} \cdot \mathbf{n}$ with an accuracy of $\varepsilon^N$ in the $\|\cdot\|_{-\frac{1}{2},\varepsilon}$ norm.

Before proving Theorem 3.5, we will first state and prove two lemmas which will be useful in the proof. The first lemma gives an expression of $\lambda_\varepsilon^N$ with respect to the jump of the gradient of $u_\varepsilon^N$ at the interface $\partial\omega_\varepsilon$.

**Lemma 3.8.** Let $(u_\varepsilon^N, \lambda_\varepsilon^N) \in H_0^1(\Omega) \times \mathcal{M}_\varepsilon^N$ be the solution of problem (2.16). Then, we have

$$ \lambda_\varepsilon^N = \varepsilon(\nabla u_{\varepsilon,ext}^N \cdot \mathbf{n} - \nabla u_{\varepsilon,int}^N \cdot \mathbf{n}) \text{ in } H^{-\frac{1}{2}}(\partial\omega_\varepsilon). \quad (3.11) $$

where $\mathbf{n}$ is the exterior normal to $\partial\omega_\varepsilon$, $u_{\varepsilon,ext}^N$ is the restriction of $u_\varepsilon^N$ to $\Omega_\varepsilon$ and $u_{\varepsilon,int}^N$ is the restriction of $u_\varepsilon^N$ to $\omega_\varepsilon$.
Proof. For all $\phi \in H_0^1(\Omega)$, we have

$$ 
\begin{aligned}
\int_{\Omega} \nabla u_{\varepsilon}^{N} \nabla \phi \, dx &= \int_{\Omega_{\varepsilon}} \nabla u_{\varepsilon}^{N} \nabla \phi \, dx + \int_{\omega_{\varepsilon}} \nabla u_{\varepsilon}^{N} \nabla \phi \, dx \\
&= - \int_{\Omega_{\varepsilon}} \Delta u_{\varepsilon}^{N} \phi \, dx - \int_{\omega_{\varepsilon}} \Delta u_{\varepsilon}^{N} \phi \, dx + \langle \nabla u_{\varepsilon, \text{int}}^{N} \cdot \mathbf{n} - \nabla u_{\varepsilon, \text{ext}}^{N} \cdot \mathbf{n}, \phi \rangle_{-\frac{1}{2}, \partial \omega_{\varepsilon}}
\end{aligned}
$$

and since $-\Delta u_{\varepsilon}^{N} = f$ in $\Omega_{\varepsilon}$ and $-\Delta u_{\varepsilon}^{N} = 0$ in $\omega_{\varepsilon}$ in a strong sense, we obtain the equality

$$ \int_{\Omega} \nabla u_{\varepsilon}^{N} \nabla \phi \, dx = \langle \nabla u_{\varepsilon, \text{int}}^{N} \cdot \mathbf{n} - \nabla u_{\varepsilon, \text{ext}}^{N} \cdot \mathbf{n}, \phi \rangle_{-\frac{1}{2}, \partial \omega_{\varepsilon}} + \int_{\Omega} f \phi \, dx. $$

On the other hand, we have

$$ \int_{\Omega} \nabla u_{\varepsilon}^{N} \nabla \phi \, dx = -(\lambda_{\varepsilon}^{N}, \phi)_{\varepsilon} + \int_{\Omega} f \phi \, dx. $$

Identifying these two formulations, we get

$$ \langle \nabla u_{\varepsilon, \text{ext}}^{N} \cdot \mathbf{n} - \nabla u_{\varepsilon, \text{int}}^{N} \cdot \mathbf{n}, \phi \rangle_{-\frac{1}{2}, \partial \omega_{\varepsilon}} = (\lambda_{\varepsilon}^{N}, \phi)_{\varepsilon} = \varepsilon^{-1} (\lambda_{\varepsilon}^{N}, \phi)_{\partial \omega_{\varepsilon}} $$

for all $\phi \in H_0^1(\Omega)$. Thus, since the trace operator from $H_0^1(\Omega)$ on $H^{\frac{1}{2}}(\partial \omega_{\varepsilon})$ is surjective, we get (3.11). $\square$

The second lemma describes the behaviour of $u_{\varepsilon}^{N}$ on $\partial \omega_{\varepsilon}$.

**Lemma 3.9.** We assume that $\varphi$ belongs to $\mathcal{M}_{\varepsilon}^{N}$. There exist $(a_{\varepsilon,n})_{n \ge N+1}$ and $(b_{\varepsilon,n})_{n \ge N+1}$ such that for $\varepsilon \in (0, \Upsilon)$, if $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial \omega$, with $\theta \in [0, 2\pi[$, we have

$$ u_{\varepsilon}^{N}(\varepsilon \mathbf{x}) = \varphi(\varepsilon \mathbf{x}) + \sum_{n \ge N+1} \left(\frac{\varepsilon}{\Upsilon}\right)^{n} (a_{\varepsilon,n} \cos(n\theta) + b_{\varepsilon,n} \sin(n\theta)). \quad (3.12) $$

Moreover,

$$ \sum_{n \ge N+1} (1+n) (|a_{\varepsilon,n}|^2 + |b_{\varepsilon,n}|^2) \lesssim \|\varphi\|_{\frac{2}{5}, \partial \omega_{\varepsilon}}^2 + \|f\|_{0, \Omega}^2. $$

Proof. Let $0 < \varepsilon < \Upsilon$ be fixed. The solution $u_{\varepsilon}^{N}$ belongs to $H^1(\Omega)$ and is harmonic in $\omega_{\varepsilon}$, so for $0 < r < \varepsilon$ and $\theta \in [0, 2\pi[$, if $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial \omega$, we have

$$ u_{\varepsilon}^{N}(r\mathbf{x}) = A_{\varepsilon,0} + \sum_{n \in \mathbb{N}^{*}} r^{n} (A_{\varepsilon,n} \cos(n\theta) + B_{\varepsilon,n} \sin(n\theta)), \quad (3.13) $$

with for $1 \le n \le N$,

$$ \left\{ 
\begin{array}{l}
A_{\varepsilon,0} = \mathcal{A}_{\varepsilon}^{0}(\varphi), \\
A_{\varepsilon,n} = \varepsilon^{-n} \mathcal{A}_{\varepsilon}^{n}(\varphi), \\
B_{\varepsilon,n} = \varepsilon^{-n} \mathcal{B}_{\varepsilon}^{n}(\varphi),
\end{array} 
\right. 
\quad (3.14) $$

where the operators $\mathcal{A}_{\varepsilon}^{n}$ and $\mathcal{B}_{\varepsilon}^{n}$ are defined by (2.8). Moreover, the solution $u_{\varepsilon}^{N}$ is also harmonic in the annulus $\omega_{\Upsilon} \setminus \bar{\omega}_{\varepsilon}$, so for $\varepsilon < r < \Upsilon$ and $\theta \in [0, 2\pi[$, if $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial \omega$, we have

$$ u_{\varepsilon}^{N}(r\mathbf{x}) = C_{\varepsilon,0} + D_{\varepsilon,0} \ln(r) + \sum_{n \in \mathbb{N}^{*}} (C_{\varepsilon,n} r^{n} + D_{\varepsilon,n} r^{-n}) \cos(n\theta) + (E_{\varepsilon,n} r^{n} + F_{\varepsilon,n} r^{-n}) \sin(n\theta). \quad (3.15) $$

with for all $1 \le n \le N$,

$$ \left\{ 
\begin{array}{l}
C_{\varepsilon,0} + D_{\varepsilon,0} \ln(\varepsilon) = \mathcal{A}_{\varepsilon}^{0}(\varphi), \\
C_{\varepsilon,n} \varepsilon^{n} + D_{\varepsilon,n} \varepsilon^{-n} = \mathcal{A}_{\varepsilon}^{n}(\varphi), \\
E_{\varepsilon,n} \varepsilon^{n} + F_{\varepsilon,n} \varepsilon^{-n} = \mathcal{B}_{\varepsilon}^{n}(\varphi).
\end{array} 
\right. $$

We refer to [24] for the derivation of formulas (3.13) and (3.15). By applying for all $n \ge N + 1$ the operators $\mathcal{A}_{\varepsilon}^{n}$ and $\mathcal{B}_{\varepsilon}^{n}$ to equations (3.13) and (3.15), we obtain that, for all $n \ge N + 1$

$$ \left\{ 
\begin{array}{l}
A_{\varepsilon,n} \varepsilon^n = C_{\varepsilon,n} \varepsilon^n + D_{\varepsilon,n} \varepsilon^{-n}, \\
B_{\varepsilon,n} \varepsilon^n = E_{\varepsilon,n} \varepsilon^n + F_{\varepsilon,n} \varepsilon^{-n}.
\end{array} 
\right. 
\quad (3.16) $$
Moreover, according to Lemma 3.8, $\lambda_{\varepsilon}^N$ satisfies (3.11). Thanks to formula (3.13), we have in $H^{-\frac{1}{2}}(\partial\omega_{\varepsilon})$

$$ \nabla u_{\varepsilon, int}^N \cdot \mathbf{n} = \frac{\partial u_{\varepsilon, int}^N}{\partial r} = \sum_{n \in \mathbb{N}^*} n \varepsilon^{n-1} (A_{\varepsilon, n} \cos(n\theta) + B_{\varepsilon, n} \sin(n\theta)), $$

and, thanks to formula (3.15), we have in $H^{-\frac{1}{2}}(\partial\omega_{\varepsilon})$

$$ \nabla u_{\varepsilon, ext}^N \cdot \mathbf{n} = \frac{\partial u_{\varepsilon, ext}^N}{\partial r} = \frac{D_{\varepsilon, 0}}{r} + \sum_{n \in \mathbb{N}^*} n(C_{\varepsilon, n}\varepsilon^{n-1} - D_{\varepsilon, n}\varepsilon^{-n-1}) \cos(n\theta) \\ + n(E_{\varepsilon, n}\varepsilon^{n-1} - F_{\varepsilon, n}\varepsilon^{-n-1}) \sin(n\theta). $$

Thus, by applying for all $n \ge N+1$ the operators $A_{\Upsilon}^n$ and $B_{\Upsilon}^n$ to equation (3.11), we obtain the following equations: for all $n \ge N+1$,

$$ \left\{ \begin{array}{l} A_{\varepsilon, n}\varepsilon^n = C_{\varepsilon, n}\varepsilon^n - D_{\varepsilon, n}\varepsilon^{-n} \\ B_{\varepsilon, n}\varepsilon^n = E_{\varepsilon, n}\varepsilon^n - F_{\varepsilon, n}\varepsilon^{-n}. \end{array} \right. \qquad (3.17) $$

From equations (3.16) and (3.17), we thus deduce that, for all $n \ge N+1$,

$$ \left\{ \begin{array}{l} D_{\varepsilon, n} = F_{\varepsilon, n} = 0 \\ A_{\varepsilon, n} = C_{\varepsilon, n} \\ B_{\varepsilon, n} = E_{\varepsilon, n}. \end{array} \right. $$

By applying now the operators $A_{\Upsilon}^n$ and $B_{\Upsilon}^n$ on equation (3.15), we obtain that, for all $n \ge N+1$,

$$ \left\{ \begin{array}{l} C_{\varepsilon, n}\Upsilon^n + D_{\varepsilon, n}\Upsilon^{-n} = A_{\Upsilon}^n(u_{\varepsilon}^N), \\ E_{\varepsilon, n}\Upsilon^n + F_{\varepsilon, n}\Upsilon^{-n} = B_{\Upsilon}^n(u_{\varepsilon}^N). \end{array} \right. $$

From these last two systems, we get that, for all $n \ge N+1$,

$$ \left\{ \begin{array}{l} D_{\varepsilon, n} = F_{\varepsilon, n} = 0 \\ A_{\varepsilon, n} = C_{\varepsilon, n} = \frac{1}{\Upsilon^n} A_{\Upsilon}^n(u_{\varepsilon}^N), \\ B_{\varepsilon, n} = E_{\varepsilon, n} = \frac{1}{\Upsilon^n} B_{\Upsilon}^n(u_{\varepsilon}^N). \end{array} \right. $$

Using (3.13) and (3.14), if we define $a_{\varepsilon, n} = A_{\Upsilon}^n(u_{\varepsilon}^N)$ and $b_{\varepsilon, n} = B_{\Upsilon}^n(u_{\varepsilon}^N)$ for all $n \ge N+1$, we get (3.12). Using now the equivalence of the norms $\|\cdot\|_{\frac{1}{2}, \partial\omega_{\Upsilon}}$ and $\|\cdot\|_{\frac{1}{2}, \partial\omega_{\Upsilon}}$ as well as Theorem 3.1, we get that

$$ \sum_{n \ge N+1} (1+n)(|a_{\varepsilon, n}|^2 + |b_{\varepsilon, n}|^2) = \sum_{n \ge N+1} (1+n)(|A_{\Upsilon}^n(u_{\varepsilon}^N)|^2 + |B_{\Upsilon}^n(u_{\varepsilon}^N)|^2) \lesssim \|u_{\varepsilon}^N\|_{\frac{1}{2}, \partial\omega_{\Upsilon}}^2 \\ \lesssim \|u_{\varepsilon}^N\|_{1, \Omega}^2 \lesssim \|\varphi\|_{\frac{1}{2}, \varepsilon}^2 + \|f\|_{0, \Omega}^2. $$

□

Let us notice that, following the computations made in this proof, we can obtain an expression of $u_{\varepsilon}^N$ in $\omega_{\Upsilon}$ with respect to its Fourier coefficients $A_{\Upsilon}^n(u_{\varepsilon}^N)$ and $B_{\Upsilon}^n(u_{\varepsilon}^N)$ on $\partial\omega_{\Upsilon}$ and the first $N+1$ Fourier coefficients of $\varphi$ on $\partial\omega_{\varepsilon}$. This expression will be useful in the next section for the analysis of the numerical approximation but we state it at this stage. We have for $0 < r < \Upsilon$ and $\theta \in [0, 2\pi[$, if $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial\omega$,

$$ \begin{aligned} u_{\varepsilon}^N(r\mathbf{x}) = & \left( A_{\Upsilon}^0(u_{\varepsilon}^N) \frac{\ln(r/\varepsilon)}{\ln(\Upsilon/\varepsilon)} + \sum_{n=1}^N A_{\Upsilon}^n(u_{\varepsilon}^N) \frac{\left(\frac{r}{\varepsilon}\right)^n - \left(\frac{\varepsilon}{r}\right)^n}{\left(\frac{\Upsilon}{\varepsilon}\right)^n - \left(\frac{\varepsilon}{\Upsilon}\right)^n} \cos(n\theta) + B_{\Upsilon}^n(u_{\varepsilon}^N) \frac{\left(\frac{r}{\varepsilon}\right)^n - \left(\frac{\varepsilon}{r}\right)^n}{\left(\frac{\Upsilon}{\varepsilon}\right)^n - \left(\frac{\varepsilon}{\Upsilon}\right)^n} \sin(n\theta) \right. \\ & + A_{\varepsilon}^0(\varphi) \frac{\ln(r/\Upsilon)}{\ln(\varepsilon/\Upsilon)} + \sum_{n=1}^N A_{\varepsilon}^n(\varphi) \frac{\left(\frac{r}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{r}\right)^n}{\left(\frac{\varepsilon}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{\varepsilon}\right)^n} \cos(n\theta) + B_{\varepsilon}^n(\varphi) \frac{\left(\frac{r}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{r}\right)^n}{\left(\frac{\varepsilon}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{\varepsilon}\right)^n} \sin(n\theta) \left. \right) \mathbf{1}_{r \ge \varepsilon} \\ & + \left( A_{\varepsilon}^0(\varphi) + \sum_{n=1}^N \left(\frac{r}{\varepsilon}\right)^n \left( A_{\varepsilon}^n(\varphi) \cos(n\theta) + B_{\varepsilon}^n(\varphi) \sin(n\theta) \right) \right) \mathbf{1}_{r < \varepsilon} \\ & + \sum_{n \ge N+1} \left(\frac{r}{\Upsilon}\right)^n \left( A_{\Upsilon}^n(u_{\varepsilon}^N) \cos(n\theta) + B_{\Upsilon}^n(u_{\varepsilon}^N) \sin(n\theta) \right). \end{aligned} \qquad (3.18) $$
Let us note that expression (3.18) is valid for any $\varphi \in H^{\frac{1}{2}}(\partial\omega_{\varepsilon})$.

We are now ready to prove Theorem 3.5.

**Proof of Theorem 3.5.** Let $0 < \varepsilon < \Upsilon$ be fixed. Let us first note that $u_{\varepsilon}^N - u_{\varepsilon}$ satisfies (3.3) with $z = u_{\varepsilon}^N - \varphi$ and so, according to Lemma 3.4, we have

$$ \|u_{\varepsilon}^N - u_{\varepsilon}\|_{1,\Omega} \lesssim \|u_{\varepsilon}^N - \varphi\|_{\frac{1}{2},\varepsilon} \lesssim \|u_{\varepsilon}^N - \varphi\|_{\frac{3}{2},\varepsilon}. \quad (3.19) $$

Then, according to Lemma 3.9, if $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial\omega$, with $\theta \in [0, 2\pi]$, we have

$$ u_{\varepsilon}^{N}(\varepsilon\mathbf{x}) = \varphi(\varepsilon\mathbf{x}) + \sum_{n=N+1}^{+\infty} \left(\frac{\varepsilon}{\Upsilon}\right)^{n} \left( a_{\varepsilon,n} \cos(n\theta) + b_{\varepsilon,n} \sin(n\theta) \right) $$

with

$$ \sum_{n=N+1}^{+\infty} (1+n)(|a_{\varepsilon,n}|^2 + |b_{\varepsilon,n}|^2) \lesssim \|\varphi\|_{\frac{1}{2},\varepsilon}^2 + \|f\|_{0,\Omega}^2. $$

Moreover, thanks to (2.12), we have

$$ \begin{aligned} \overline{\|u_{\varepsilon}^N - \varphi\|}_{\frac{1}{2},\varepsilon} & \lesssim \left(\frac{\varepsilon}{\Upsilon}\right)^{N+1} \left( \sum_{n=N+1}^{+\infty} (1+n) \left(\frac{\varepsilon}{\Upsilon}\right)^{2(n-N-1)} (|a_{\varepsilon,n}|^2 + |b_{\varepsilon,n}|^2) \right)^{\frac{1}{2}} \\ & \lesssim \left(\frac{\varepsilon}{\Upsilon}\right)^{N+1} \left( \sum_{n=N+1}^{+\infty} (1+n)(|a_{\varepsilon,n}|^2 + |b_{\varepsilon,n}|^2) \right)^{\frac{1}{2}} \lesssim \left(\frac{\varepsilon}{\Upsilon}\right)^{N+1} \left( \|\varphi\|_{\frac{1}{2},\varepsilon} + \|f\|_{0,\Omega} \right). \end{aligned} $$

Using this inequality in (3.19) allows to get (3.8).

Let us now prove (3.9). Using the weak formulations (2.5) and (2.16) and the formulas (2.13) and (2.14) for the relationship between $\langle \cdot, \cdot \rangle_{-\frac{1}{2},\partial\omega_{\varepsilon}}$, $\langle \cdot, \cdot \rangle_{-\frac{1}{2},\varepsilon}$ and $(\cdot, \cdot)_{\varepsilon}$, we notice that $(u_{\varepsilon} - u_{\varepsilon}^N, \lambda_{\varepsilon} - \varepsilon^{-1}\lambda_{\varepsilon}^N)$ verifies

$$ (\nabla(u_{\varepsilon} - u_{\varepsilon}^N), \nabla v)_{\Omega} + \varepsilon \langle \lambda_{\varepsilon} - \varepsilon^{-1} \lambda_{\varepsilon}^N, v \rangle_{-\frac{1}{2},\varepsilon} = 0, \forall v \in H_0^1(\Omega). $$

Next, we have that

$$ \|\lambda_{\varepsilon} - \varepsilon^{-1}\lambda_{\varepsilon}^N\|_{-\frac{1}{2},\varepsilon} = \sup_{\mu \in H^{\frac{1}{2}}(\partial\omega_{\varepsilon})} \frac{\langle \lambda_{\varepsilon} - \varepsilon^{-1}\lambda_{\varepsilon}^N, \mu \rangle_{-\frac{1}{2},\varepsilon}}{\|\mu\|_{\frac{1}{2},\varepsilon}}. $$

For all $\mu \in H^{\frac{1}{2}}(\partial\omega_{\varepsilon})$, if we take $v \in H_0^1(\Omega)$ defined in Lemma 3.4 such that $v = \mu$ on $\partial\omega_{\varepsilon}$, we have

$$ \|\lambda_{\varepsilon} - \varepsilon^{-1}\lambda_{\varepsilon}^N\|_{-\frac{1}{2},\varepsilon} \lesssim \sup_{v \in H_0^1(\Omega)} \frac{\langle \lambda_{\varepsilon} - \varepsilon^{-1}\lambda_{\varepsilon}^N, v \rangle_{-\frac{1}{2},\varepsilon}}{\|v\|_{1,\Omega}}. $$

We deduce that

$$ \|\lambda_{\varepsilon} - \varepsilon^{-1}\lambda_{\varepsilon}^N\|_{-\frac{1}{2},\varepsilon} \lesssim \frac{1}{\varepsilon} \sup_{v \in H_0^1(\Omega)} \frac{(\nabla(u_{\varepsilon} - u_{\varepsilon}^N), \nabla v)}{\|v\|_{1,\Omega}} \lesssim \frac{1}{\varepsilon} \|u_{\varepsilon} - u_{\varepsilon}^N\|_{1,\Omega}, $$

and conclude, thanks to (3.8), that (3.9) holds. $\square$

### 3.3 Regularity of the solution

Due to the presence of the Dirac source $\lambda_{\varepsilon}^N \delta_{\partial\omega_{\varepsilon}}$ in problem (2.17), the global $H^2$ regularity for $u_{\varepsilon}^N$ cannot be obtained. However, we will prove that the solution $u_{\varepsilon}^N$ is more regular than $H^1$. To do so, we will rewrite the strong formulation (2.17) in the following form: $u_{\varepsilon}^N \in H_0^1(\Omega)$ is solution of

$$ \begin{cases} -\Delta u_{\varepsilon}^N = f - \varepsilon^{-1} \lambda_{\varepsilon}^N \delta_{\partial\omega_{\varepsilon}} & \text{in } \Omega, \\ u_{\varepsilon}^N = 0 & \text{on } \partial\Omega. \end{cases} $$

**Theorem 3.10.** Let $0 < \eta < \frac{1}{2}$ be given. Then, if $\varphi \in H^1(\partial\omega_{\varepsilon})$, the solution of problem (2.16) satisfies the additional regularity $u_{\varepsilon}^N \in H^{\frac{3}{2}-\eta}(\Omega)$ and the following estimate holds:

$$ \|u_{\varepsilon}^N\|_{\frac{3}{2}-\eta,\Omega} \lesssim_{\varepsilon} \|\varphi\|_{1,\partial\omega_{\varepsilon}} + \|f\|_{0,\Omega}. \quad (3.20) $$
To prove this theorem, we will need an auxiliary lemma presented in [18].

**Lemma 3.11.** Let $\mathcal{D}$ be a bounded, convex domain in $\mathbb{R}^2$. Let $\gamma \subset \mathcal{D}$ be a $\mathcal{C}^2$-surface such that the distance between $\gamma$ and $\partial\mathcal{D}$ is positive. We consider the following problem

$$ \begin{cases} -\Delta y = \zeta \delta_{\gamma} & \text{in } \mathcal{D}, \\ y = 0 & \text{on } \partial\mathcal{D}, \end{cases} \qquad (3.21) $$

where $\zeta \in L^2(\gamma)$. Then problem (3.21) admits a unique solution $y$ which belongs to $H^{\frac{3}{2}-\eta}(\mathcal{D})$ for any $\eta > 0$. Furthermore there exists a constant $C$ such that

$$ \|y\|_{\frac{3}{2}-\eta, \mathcal{D}} \leq C \|\zeta \delta_{\gamma}\|_{-\frac{1}{2}-\eta, \mathcal{D}} $$

where $C$ is independent of $\zeta$ and $\gamma$.

**Proof of Theorem 3.10.** Since $\varepsilon^{-1}\lambda_{\varepsilon}^N \in L^2(\partial\omega_{\varepsilon})$, we can apply Lemma 3.11 and get that $u_{\varepsilon}^N \in H^{\frac{3}{2}-\eta}(\Omega)$. Then the main objective of the proof is to show the estimate (3.20). We start observing that thanks to Lemma 3.11 we also have,

$$ \|u_{\varepsilon}^N\|_{\frac{3}{2}-\eta, \Omega} \lesssim \varepsilon^{-1} \|\lambda_{\varepsilon}^N \delta_{\partial\omega_{\varepsilon}}\|_{-\frac{1}{2}-\eta, \Omega}, \qquad (3.22) $$

and that

$$ \|\lambda_{\varepsilon}^N \delta_{\partial\omega_{\varepsilon}}\|_{-\frac{1}{2}-\eta, \Omega} = \sup_{v \in H_0^{\frac{1}{2}+\eta}(\Omega)} \frac{(\lambda_{\varepsilon}^N, v)_{\partial\omega_{\varepsilon}}}{\|v\|_{\frac{1}{2}+\eta, \Omega}}. $$

Note that here $(\cdot, \cdot)_{\partial\omega_{\varepsilon}}$ is the standard scalar product on $L^2(\partial\omega_{\varepsilon})$. Using Cauchy-Schwarz inequality, we have, for all $v \in H_0^{\frac{1}{2}+\eta}(\Omega)$,

$$ (\lambda_{\varepsilon}^N, v)_{\partial\omega_{\varepsilon}} \leq \|\lambda_{\varepsilon}^N\|_{0, \partial\omega_{\varepsilon}} \|v\|_{0, \partial\omega_{\varepsilon}}. $$

Thanks to the classical trace inequality, we have, for all $v \in H_0^{\frac{1}{2}+\eta}(\Omega)$,

$$ \|v\|_{0, \partial\omega_{\varepsilon}} \lesssim_{\varepsilon} \|v\|_{\frac{1}{2}+\eta, \Omega}. \qquad (3.23) $$

Thus, we get that

$$ \|\lambda_{\varepsilon}^N \delta_{\partial\omega_{\varepsilon}}\|_{-\frac{1}{2}-\eta, \Omega} \lesssim_{\varepsilon} \|\lambda_{\varepsilon}^N\|_{0, \partial\omega_{\varepsilon}}. $$

In addition, thanks to equality (3.11) as well as expression (3.18), we can have an explicit expression of $\lambda_{\varepsilon}^N$ on $\partial\omega_{\varepsilon}$ according to the first $N+1$ Fourier moments of $u_{\varepsilon}^N$ on $\partial\omega_{\Upsilon}$ and $\varphi$ on $\partial\omega_{\varepsilon}$. Indeed, after some manipulations, we can prove that

$$ \|\lambda_{\varepsilon}^N\|_{0, \partial\omega_{\varepsilon}} \lesssim_{\varepsilon} \|u_{\varepsilon}^N\|_{1, \partial\omega_{\Upsilon}} + \|\varphi\|_{1, \partial\omega_{\varepsilon}}. $$

Using successively trace inequality and elliptic regularity on the domain $\Omega_{\Upsilon} = \Omega \setminus \omega_{\Upsilon}$, we get

$$ \|u_{\varepsilon}^N\|_{1, \partial\omega_{\Upsilon}} \lesssim \|u_{\varepsilon}^N\|_{2, \Omega_{\Upsilon}} \lesssim \|f\|_{0, \Omega_{\Upsilon}} \leq \|f\|_{0, \Omega} $$

Then, if $\varphi \in H^1(\partial\omega_{\varepsilon})$,

$$ \|\lambda_{\varepsilon}^N\|_{0, \partial\omega_{\varepsilon}} \lesssim_{\varepsilon} \|\varphi\|_{1, \partial\omega_{\varepsilon}} + \|f\|_{0, \Omega}, $$

where the constant of the inequality depends a-priori on $\varepsilon$. Using these last two inequalities and coming back to (3.22), we deduce that

$$ \|u_{\varepsilon}^N\|_{\frac{3}{2}-\eta, \Omega} \lesssim_{\varepsilon} \|\varphi\|_{1, \partial\omega_{\varepsilon}} + \|f\|_{0, \Omega}. $$

□

Note that the dependence with respect to $\varepsilon$ in the estimate (3.20) is not explicit. This comes directly from the fact that the dependence with respect to $\varepsilon$ of the constant in the trace inequality (3.23) is not known itself.
# 4 Numerical approximation

We aim at approximating problem (2.16) using the finite element method. Interestingly, we just need to discretize the space $H_0^1(\Omega)$ as the space $\mathcal{M}_\varepsilon^N$ is already of finite dimension. Let us introduce a family of shape regular, quasi-uniform triangulations $\mathcal{T}_h^\Omega$ of $\Omega$, where $h$ is the characteristic size of the mesh. We set for $k \ge 1$,

$$X_h^k(\Omega) = \{v_h \in \mathcal{C}(\Omega) \cap H_0^1(\Omega) \mid v_h|_K \in \mathbb{P}^k, \forall K \in \mathcal{T}_h^\Omega\},$$

where $\mathbb{P}^k$ is the vector space of polynomials of degree less than or equal to $k$. Then, the discrete version of problem (2.16) writes, find $(u_{\varepsilon,h}^N, \lambda_{\varepsilon,h}^N) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N$ such that

$$\left\{ \begin{array}{ll} \left(\nabla u_{\varepsilon,h}^N, \nabla v_h\right)_\Omega + \left(\lambda_{\varepsilon,h}^N, v_h\right)_\varepsilon = (f, v_h)_\Omega, & \forall v_h \in X_h^k(\Omega), \\ (\mu, u_{\varepsilon,h}^N)_\varepsilon = (\mu, \varphi)_\varepsilon, & \forall \mu \in \mathcal{M}_\varepsilon^N. \end{array} \right. \quad (4.1)$$

Problem (4.1) is a saddle point problem that will be analyzed in the framework of the available general theory, see for example [5, 13]. In particular the central property is the inf-sup stability of the form $(\lambda_{\varepsilon,h}^N, v_h)_\varepsilon$ that will be discussed in the next section.

For the approximation of problem (4.1) we introduce the Scott-Zhang operator $R_h^k : H_0^1(\Omega) \to X_h^k(\Omega)$, defined in [38]. $R_h^k$ is continuous in $H_0^1(\Omega)$, namely, for all $v \in H_0^1(\Omega)$, $\|R_h^k v\|_{1,\Omega} \lesssim \|v\|_{1,\Omega}$ and strongly consistent, that is $R_h^k v = v$ for all $v \in X_h^k(\Omega)$. The Scott-Zhang operator satisfies the following approximation properties. Let $k \in \mathbb{N}^*$ and $s \ge 1$, for $K \in \mathcal{T}_h^\Omega$ and $v \in H^s(S_K)$, then

$$h^{-\frac{1}{2}} \|R_h^k v - v\|_{0,K} + h^{\frac{1}{2}} \|\nabla(R_h^k v - v)\|_{0,K} \lesssim h^{l+\frac{1}{2}} |v|_{s,S_K}, \quad l = \min\{k, s-1\}, \quad (4.2)$$

where $S_K$ is a domain made of the elements neighboring $K$. Moreover, we have for $v \in H_0^s(\Omega)$,

$$\|v - R_h^k v\|_{1,\Omega} \lesssim h^l |v|_{s,\Omega}, \quad l = \min\{k, s-1\}. \quad (4.3)$$

## 4.1 Inf-Sup Stability

The aim of this section is to discuss the discrete counterpart of the inf-sup condition for problem (4.1). We show below that it is satisfied under the assumption that $k \ge N$ which corresponds to the case where the polynomial degree in the discrete space $X_h^k(\Omega)$ is larger than or equal to the degree of the trigonometric polynomials in the reduced space $\mathcal{M}_\varepsilon^N$.

**Lemma 4.1.** We assume that $k \ge N$. There exists a constant $\beta > 0$ uniformly independent of $h, N$ and $\varepsilon$ such that for all $\lambda \in \mathcal{M}_\varepsilon^N$,

$$\sup_{v_h \in X_h^k(\Omega)} \frac{(\lambda, v_h)_\varepsilon}{\|v_h\|_{1,\Omega}} \ge \frac{\beta}{(1 + \frac{h}{\varepsilon})^N} \|\lambda\|_{-\frac{1}{2},\varepsilon}$$

*Proof.* Let $\mu \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ be given. We define $a_n$, for all $n \ge 0$, and $b_n$, for all $n \ge 1$, by

$$a_n = \mathcal{A}_\varepsilon^n(\mu) \text{ and } b_n = \mathcal{B}_\varepsilon^n(\mu).$$

Then we consider $v \in H_0^1(\Omega)$ such that in $\omega_{\varepsilon+h} = (\varepsilon + h)\omega$,

$$v = a_0 + \sum_{n=1}^{N} a_n \left(\frac{r}{\varepsilon}\right)^n \cos(n\theta) + b_n \left(\frac{r}{\varepsilon}\right)^n \sin(n\theta)$$

and outside $\omega_{\varepsilon+h}$, $v$ is defined as the harmonic lifting in $H_0^1(\Omega)$. We notice that $v$ satisfies equation (3.3) where $\omega_\varepsilon$ has been replaced by $\omega_{\varepsilon+h}$ and $z$ is given by $z = \Pi_\varepsilon^N(\mu)$. Thus, according to Lemma 3.4, we have

$$\|v\|_{1,\Omega} \lesssim \|v\|_{\frac{1}{2},\varepsilon+h}.$$

Moreover on $\partial\omega_{\varepsilon+h}$,

$$v = a_0 + \sum_{n=1}^{N} a_n \left(1 + \frac{h}{\varepsilon}\right)^n \cos(n\theta) + b_n \left(1 + \frac{h}{\varepsilon}\right)^n \sin(n\theta).$$
This implies that

$$
\begin{align*}
\|v\|_{\frac{1}{2}, \varepsilon+h} &\lesssim \left( a_0^2 + \sum_{n=1}^N (1+n) \left(1+\frac{h}{\varepsilon}\right)^{2n} (a_n^2+b_n^2) \right)^{\frac{1}{2}} \\
&\le \left(1+\frac{h}{\varepsilon}\right)^N \left( a_0^2 + \sum_{n \in \mathbb{N}^*} (1+n)(a_n^2+b_n^2) \right)^{\frac{1}{2}} \lesssim \left(1+\frac{h}{\varepsilon}\right)^N \|\mu\|_{\frac{1}{2}, \varepsilon}.
\end{align*}
$$

Furthermore, if we consider the trigonometric polynomials $c_n(r, \theta) = (\frac{r}{\epsilon})^n \cos(n\theta)$ and $s_n(r, \theta) = (\frac{r}{\epsilon})^n \sin(n\theta)$, for $1 \le n \le N$, we can write them as polynomials of $x$ and $y$ of degree smaller than or equal to $n$, where $(x, y)$ satisfies $(x, y) = (r \cos(\theta), r \sin(\theta))$. Indeed, using Chebyshev polynomials, we have

$$
c_n(r, \theta) = \frac{1}{\epsilon^n} \sum_{0 \le 2j \le n} \binom{n}{2j} (-1)^j y^{2j} x^{n-2j},
$$

and

$$
s_n(r, \theta) = \frac{1}{\epsilon^n} \sum_{0 \le 2j+1 \le n} \binom{n}{2j+1} (-1)^j y^{2j+1} x^{n-2j-1}.
$$

Since, by assumption, $N \le k$, we deduce that $v|_{\omega_{\epsilon+h}}$ belongs to $\mathbb{P}^k$.

Let us now set $v_h = R_h^k v$. By definition of $h$, the domain $\omega_{\epsilon+h}$ contains all the elements $K \in \mathcal{T}_h^\Omega$ that intersect $\partial \omega_\epsilon$. Thus, thanks to the properties of stability and strong consistency satisfied by $R_h^k$, we have $v_h = v = \Pi_\epsilon^N \mu$ on $\partial \omega_\epsilon$, and $\|v_h\|_{1,\Omega} \le \beta^{-1}(1+\frac{h}{\epsilon})^N \|\mu\|_{\frac{1}{2},\epsilon}$, where $\beta$ is a constant independent of $h$ and $\epsilon$ and also uniformly independent of $N$. From these two properties on $v_h$, it comes that, for any $\lambda \in \mathcal{M}_\epsilon^N$, we have

$$
\frac{(v_h, \lambda)_\epsilon}{\|v_h\|_{1,\Omega}} \ge \frac{\beta}{(1+\frac{h}{\epsilon})^N} \frac{(\Pi_\epsilon^N \mu, \lambda)_\epsilon}{\|\mu\|_{\frac{1}{2},\epsilon}} = \frac{\beta}{(1+\frac{h}{\epsilon})^N} \frac{(\mu, \lambda)_\epsilon}{\|\mu\|_{\frac{1}{2},\epsilon}}
$$

Thus, it comes that

$$
\sup_{v_h \in X_h^k(\Omega)} \frac{(v_h, \lambda)_\epsilon}{\|v_h\|_{1,\Omega}} \ge \frac{\beta}{(1+\frac{h}{\epsilon})^N} \sup_{\mu \in H^{\frac{1}{2}}(\partial \omega_\epsilon)} \frac{(\mu, \lambda)_\epsilon}{\|\mu\|_{\frac{1}{2},\epsilon}} = \frac{\beta}{(1+\frac{h}{\epsilon})^N} \|\lambda\|_{-\frac{1}{2},\epsilon}.
$$

□

The same observation applies, in what follows, to Lemmas 4.2, 4.5.

Thanks to Lemma 4.1, noting that the operators $a$ and $b$ defined in (3.6) are bounded by constants uniformly independent of $\epsilon$, we obtain the following general a-priori estimates for the error,

$$
\|u_{\epsilon}^{N} - u_{\epsilon,h}^{N}\|_{1,\Omega} \lesssim (1+\frac{h}{\epsilon})^{N} \inf_{v_h \in X_h^k(\Omega)} \|u_{\epsilon}^{N} - v_h\|_{1,\Omega}, \quad (4.4)
$$

$$
\|\lambda_{\epsilon}^{N} - \lambda_{\epsilon,h}^{N}\|_{-\frac{1}{2},\epsilon} \lesssim (1+\frac{h}{\epsilon})^{2N} \inf_{v_h \in X_h^k(\Omega)} \|u_{\epsilon}^{N} - v_h\|_{1,\Omega}. \quad (4.5)
$$

Note that we are in a particular case of the general results for saddle-point problems, because here the discretization space for the Lagrange multipliers coincides with the continuous one. As a result in this case $\lambda_{\epsilon}^N, \lambda_{\epsilon,h}^N \in \mathcal{M}_{\epsilon}^N$ and consequently the approximation error for the Lagrange multipliers vanishes, namely $\inf_{\mu \in \mathcal{M}_{\epsilon}^N} \|\lambda_{\epsilon}^N - \mu\|_{-\frac{1}{2}, \epsilon} = 0$.

Estimates (4.4)-(4.5) are not fully satisfactory due to the presence of the factor $(1 + \frac{h}{\epsilon})^N$ that comes from the inf-sup condition of Lemma 4.1. It shows that the scheme (4.1) is affected by a loss of stability when $\epsilon \to 0$ for a fixed computational mesh. This is an example of the well known locking phenomenon (see for example [5, 13]) that appears because the inclusion becomes too small with respect to the mesh size and by consequence the solution $u_{\epsilon,h}^N$ becomes over-constrained. To overcome the limitations of the scheme (4.1) with respect to stability, we propose in the next section a stabilized variant of it.
## 4.2 Stabilized finite element method

In order to remove the stability condition $k \ge N$ that appears in Lemma 4.1, we introduce in this section a stabilized formulation based on a penalty method. We look for $u_{\varepsilon,h}^N \in X_h^k(\Omega)$ and $\lambda_{\varepsilon,h}^N \in \mathcal{M}_{\varepsilon}^N$ such that

$$ \left\{ \begin{array}{ll} \left(\nabla u_{\varepsilon,h}^N, \nabla v_h\right)_\Omega + \left(\lambda_{\varepsilon,h}^N, v_h\right)_\varepsilon = (f, v_h)_\Omega, & \forall v_h \in X_h^k(\Omega), \\ (\mu, u_{\varepsilon,h}^N)_\varepsilon = (\mu, \varphi)_\varepsilon + s_h(\mu, \lambda_{\varepsilon,h}) & \forall \mu \in \mathcal{M}_\varepsilon^N, \end{array} \right. \quad (4.6) $$

where $s_h(\mu, \lambda) : \mathcal{M}_\varepsilon^N \times \mathcal{M}_\varepsilon^N \to \mathbb{R}$ is a stabilization term defined for all $(\mu, \lambda) \in \mathcal{M}_\varepsilon^N \times \mathcal{M}_\varepsilon^N$ by

$$ s_h(\mu, \lambda) = \gamma \left(\frac{h}{\varepsilon}\right) \varepsilon^{-1} (\mu, \lambda)_\varepsilon, \quad (4.7) $$

where $\gamma$ is a positive constant independent of $h$ and $\varepsilon$. We also introduce the following (augmented) bilinear form $B_h : (X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N) \times (X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N)$ defined for all $(u_h, \lambda), (v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N$ by

$$ B_h((u_h, \lambda), (v_h, \mu)) = a(u_h, v_h) + b(\lambda, v_h) + b(\mu, u_h) - s_h(\mu, \lambda). $$

Problem (4.6) is equivalent to find $(u_{\varepsilon,h}^N, \lambda_{\varepsilon,h}^N) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N$ such that

$$ B_h((u_{\varepsilon,h}^N, \lambda_{\varepsilon,h}^N), (v_h, \mu)) = (f, v_h)_\Omega + (\mu, \varphi)_\varepsilon, \quad \forall (v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N. \quad (4.8) $$

For the analysis of problem (4.6) we introduce the following discrete trace lemma (see [21]). Let $\Gamma$ be an internal smooth boundary dividing $\Omega$ into two open sets and $K \in \mathcal{T}_h^\Omega$. We set $\Gamma_K := \Gamma \cap K$. There exists a constant $C > 0$ depending only on the maximum curvature of $\Gamma$ such that for all $v \in H_0^1(\Omega)$,

$$ \|v\|_{0,\Gamma_K}^2 \le C (h^{-1} \|v\|_{0,K}^2 + h \|v\|_{1,K}^2). \quad (4.9) $$

For the analysis of the augmented bilinear form, we define the following norm:

$$ \|(v_h, \mu)\|_{\varepsilon} = \left( \|v_h\|_{1,\Omega}^2 + \|\mu\|_{-\frac{1}{2},\varepsilon}^2 + \left(\frac{h}{\varepsilon}\right) \varepsilon^{-1} \|\mu\|_{0,\varepsilon}^2 \right)^{\frac{1}{2}}, \quad \forall (v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N. $$

Using such norm, it is straightforward to show that $B_h(\cdot, \cdot)$ is continuous, precisely

$$ |B_h((u_h, \lambda), (v_h, \mu))| \lesssim \|(u_h, \lambda)\|_{\varepsilon} \|(v_h, \mu)\|_{\varepsilon}, \quad \forall (u_h, \lambda), (v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N. $$

The stability of the stabilized problem is addressed in the following lemma.

**Lemma 4.2.** There exists a constant $\beta > 0$ uniformly independent of $h$, $N$ and $\varepsilon$ such that

$$ \inf_{(u_h, \lambda) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N} \sup_{(v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N} \frac{B_h((u_h, \lambda), (v_h, \mu))}{\|(u_h, \lambda)\|_{\varepsilon} \|(v_h, \mu)\|_{\varepsilon}} \ge \beta. $$

*Proof.* Let $(u_h, \lambda)$ be given in $X_h^k(\Omega) \times \mathcal{M}_\varepsilon^N$. First, taking $(v_h, \mu) = (u_h, -\lambda)$ in the definition of $B_h$, we observe that

$$ B_h((u_h, \lambda), (u_h, -\lambda)) = \|\nabla u_h\|_{0,\Omega}^2 + \gamma \left(\frac{h}{\varepsilon}\right) \varepsilon^{-1} \|\lambda\|_{0,\varepsilon}^2. \quad (4.10) $$

Next, using definition (2.13) of the norm $\|\cdot\|_{-\frac{1}{2},\varepsilon}$, since $\mathcal{M}_\varepsilon^N$ is a closed subspace of $L^2(\partial\omega_\varepsilon)$, there exists $\mu \in \mathcal{M}_\varepsilon^N$ such that $\|\mu\|_{\frac{1}{2},\varepsilon} = 1$ and

$$ \|\lambda\|_{-\frac{1}{2},\varepsilon} = (\lambda, \mu)_\varepsilon. $$

Taking now $\underline{\mu} = \mu\|\lambda\|_{-\frac{1}{2},\varepsilon}$, we have $\|\underline{\mu}\|_{\frac{1}{2},\varepsilon} = \|\lambda\|_{-\frac{1}{2},\varepsilon}$ and

$$ \|\lambda\|_{-\frac{1}{2},\varepsilon}^2 = (\lambda, \underline{\mu})_\varepsilon. $$

According to Lemma 3.4, there exists $\underline{v}_\varepsilon \in H_0^1(\Omega)$ such that $\underline{v}_\varepsilon = \underline{\mu}$ on $\partial\omega_\varepsilon$ and

$$ \|\underline{v}_\varepsilon\|_{1,\Omega} \lesssim \|\underline{\mu}\|_{\frac{1}{2},\varepsilon} = \|\lambda\|_{-\frac{1}{2},\varepsilon}. \quad (4.11) $$
We will prove that, for $\alpha > 0$ small enough, there exists $C_m > 0$ such that

$$ B_h((u_h, \lambda), (u_h - \alpha R_h^k \underline{v}_\epsilon, -\lambda)) \geq C_m \|(u_h, \lambda)\|_{\epsilon}^2. \quad (4.12) $$

Property (4.9) implies that for the rescaled function $\Psi_{\epsilon}(\underline{v}_{\epsilon})$, we have, for all $K \in \mathcal{T}_h$ such that $K \cap \partial\omega_{\epsilon} \neq \emptyset$,

$$ \|R_h^k \Psi_{\epsilon}(\underline{v}_{\epsilon}) - \Psi_{\epsilon}(\underline{v}_{\epsilon})\|_{0, \frac{K}{\epsilon} \cap \partial \omega_{\epsilon}}^2 \lesssim \left(\frac{h}{\epsilon}\right)^{-1} \|R_h^k \Psi_{\epsilon}(\underline{v}_{\epsilon}) - \Psi_{\epsilon}(\underline{v}_{\epsilon})\|_{0, \frac{K}{\epsilon} + \left(\frac{h}{\epsilon}\right)}^2 + \left(\frac{h}{\epsilon}\right) \|\nabla(R_h^k \Psi_{\epsilon}(\underline{v}_{\epsilon}) - \Psi_{\epsilon}(\underline{v}_{\epsilon}))\|_{0, \frac{K}{\epsilon}}^2. $$

Using the approximation property of $R_h^k$ given in (4.2), we obtain that

$$ \|R_h^k \Psi_{\epsilon}(\underline{v}_{\epsilon}) - \Psi_{\epsilon}(\underline{v}_{\epsilon})\|_{0, \frac{K}{\epsilon} \cap \partial \omega_{\epsilon}}^2 \lesssim \left(\frac{h}{\epsilon}\right) \|\nabla \Psi_{\epsilon}(\underline{v}_{\epsilon})\|_{0, \frac{S_K}{\epsilon}}^2. $$

A scaling argument then leads to

$$ \|R_h^k \underline{v}_{\epsilon} - \underline{v}_{\epsilon}\|_{0, K \cap \partial \omega_{\epsilon}}^2 \lesssim \left(\frac{h}{\epsilon}\right) \|\nabla \underline{v}_{\epsilon}\|_{0, S_K}^2. $$

Summing over all tiles $K$ intersecting $\partial\omega_{\epsilon}$, we get, according to (4.11)

$$ \|R_h^k \underline{v}_{\epsilon} - \underline{v}_{\epsilon}\|_{0, \partial \omega_{\epsilon}} \lesssim \left(\frac{h}{\epsilon}\right)^{\frac{1}{2}} \|\underline{v}_{\epsilon}\|_{1, \Omega} \lesssim \left(\frac{h}{\epsilon}\right)^{\frac{1}{2}} \|\lambda\|_{-\frac{1}{2}, \epsilon}. $$

Then using Cauchy-Schwarz inequality and the scaling between the norms $\| \cdot \|_{0, \epsilon}$ and $\| \cdot \|_{0, \partial \omega_{\epsilon}}$ deduced from (2.14), we get

$$ \begin{aligned} |(\lambda, R_h^k \underline{v}_{\epsilon} - \underline{v}_{\epsilon})_{\epsilon}| &\leq C \left(\frac{h}{\epsilon}\right)^{\frac{1}{2}} \epsilon^{-\frac{1}{2}} \|\lambda\|_{0, \epsilon} \|\lambda\|_{-\frac{1}{2}, \epsilon}, \\ &\leq \frac{1}{4} \|\lambda\|_{-\frac{1}{2}, \epsilon}^2 + C \left(\frac{h}{\epsilon}\right) \epsilon^{-1} \|\lambda\|_{0, \epsilon}^2. \end{aligned} $$

Moreover, using the continuity of $R_h^k$ and (4.11), we have

$$ |(\nabla u_h, \nabla R_h^k \underline{v}_{\epsilon})_{\Omega}| \lesssim \|\nabla u_h\|_{0, \Omega} \|\lambda\|_{-\frac{1}{2}, \epsilon} \leq C \|\nabla u_h\|_{0, \Omega}^2 + \frac{1}{4} \|\lambda\|_{-\frac{1}{2}, \epsilon}^2. $$

Combining these last two results and taking $(v_h, \mu) = (R_h^k \underline{v}_{\epsilon}, 0)$ in the definition of $B_h$, we have

$$ \begin{aligned} B_h((u_h, \lambda), (R_h^k \underline{v}_{\epsilon}, 0)) &= (\nabla u_h, \nabla R_h^k \underline{v}_{\epsilon})_{\Omega} + (\lambda, R_h^k \underline{v}_{\epsilon})_{\epsilon}, \\ &= (\nabla u_h, \nabla R_h^k \underline{v}_{\epsilon})_{\Omega} + (\lambda, R_h^k \underline{v}_{\epsilon} - \underline{v}_{\epsilon})_{\epsilon} + (\lambda, \underline{v}_{\epsilon})_{\epsilon}, \\ &\geq -C \|\nabla u_h\|_{0, \Omega}^2 + \frac{1}{2} \|\lambda\|_{-\frac{1}{2}, \epsilon}^2 - C \left(\frac{h}{\epsilon}\right) \epsilon^{-1} \|\lambda\|_{0, \epsilon}^2. \end{aligned} $$

Gathering finally this inequality and (4.10), we get

$$ \begin{aligned} B_h((u_h, \lambda), (u_h + \alpha R_h^k \underline{v}_{\epsilon}, -\lambda)) \geq & \|\nabla u_h\|_{0, \Omega}^2 + \gamma \epsilon^{-1} \left(\frac{h}{\epsilon}\right) \|\lambda\|_{0, \epsilon}^2 \\ & + \alpha(-C \|\nabla u_h\|_{0, \Omega}^2 + \frac{1}{2} \|\lambda\|_{-\frac{1}{2}, \epsilon}^2 - C \left(\frac{h}{\epsilon}\right) \epsilon^{-1} \|\lambda\|_{0, \epsilon}^2), \end{aligned} $$

which leads to (4.12) if we take $\alpha > 0$ small enough with

$$ C_m = \min \left( 1 - \alpha C, \frac{1}{2}\alpha, \gamma - \alpha C \right). $$

To end the proof, we notice that, from the triangle inequality, the continuity of $R_h^k$ and inequality (4.11), we have

$$ \|(u_h - \alpha R_h^k \underline{v}_{\epsilon}, -\lambda)\|_{\epsilon} \leq \|(u_h, \lambda)\|_{\epsilon} + \alpha C \|\lambda\|_{-\frac{1}{2}, \epsilon} \leq C_M \|(u_h, \lambda)\|_{\epsilon}, \quad (4.13) $$

where

$$ C_M = 1 + \alpha C. $$
Thus, dividing (4.12) by (4.13), we get, with $v_h = u_h - \alpha R_h^k v_\epsilon$ and $\mu = -\lambda$

$$ \frac{B_h((u_h, \lambda), (v_h, \mu))}{\|(v_h, \mu)\|_{\epsilon}} \geq \frac{C_m}{C_M} \|(u_h, \lambda)\|_{\epsilon}, $$

which implies that, for all $(u_h, \lambda) \in X_h^k(\Omega) \times \mathcal{M}_{\epsilon}^N$,

$$ \sup_{(v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_{\epsilon}^N} \frac{B_h((u_h, \lambda), (v_h, \mu))}{\|(v_h, \mu)\|_{\epsilon}} \geq \frac{C_m}{C_M} \|(u_h, \lambda)\|_{\epsilon}, $$

leading to the desired result.

<figure><img src="image_3.png" /></figure>

To conclude the analysis of the stabilized formulation, we address the consistency of the method. Since $\lambda_{\epsilon,h}^N - \lambda_{\epsilon}^N \in \mathcal{M}_{\epsilon}^N$, we can apply Lemma 4.2 to have an estimate of $(u_{\epsilon,h}^N - R_h^k u_{\epsilon}^N, \lambda_{\epsilon,h}^N - \lambda_{\epsilon}^N)$. Precisely, we have,

$$ \|(u_{\epsilon,h}^N - R_h^k u_{\epsilon}^N, \lambda_{\epsilon,h}^N - \lambda_{\epsilon}^N)\|_{\epsilon} \lesssim \sup_{(v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_{\epsilon}^N} \frac{B_h((u_{\epsilon,h}^N - R_h^k u_{\epsilon}^N, \lambda_{\epsilon,h}^N - \lambda_{\epsilon}^N), (v_h, \mu))}{\|(v_h, \mu)\|_{\epsilon}}. $$

Then, by combining the equations of the continuous problem (3.5) and the discrete problem (4.8), we obtain that for all $(v_h, \mu) \in X_h^k(\Omega) \times \mathcal{M}_{\epsilon}^N$

$$ B_h((u_{\epsilon,h}^N - R_h^k u_{\epsilon}^N, \lambda_{\epsilon,h}^N - \lambda_{\epsilon}^N), (v_h, \mu)) = a(u_{\epsilon}^N - R_h^k u_{\epsilon}^N, v_h) + b(\mu, u_{\epsilon}^N - R_h^k u_{\epsilon}^N) + s_h(\lambda_{\epsilon}^N, \mu). $$

Now, using Lemma 3.3, we deduce that

$$ \|(u_{\epsilon,h}^N - R_h^k u_{\epsilon}^N, \lambda_{\epsilon,h}^N - \lambda_{\epsilon}^N)\|_{\epsilon} \lesssim \|(u_{\epsilon}^N - R_h^k u_{\epsilon}^N, 0)\|_{\epsilon} + \frac{h^{\frac{1}{2}}}{\epsilon} \|\lambda_{\epsilon}^N\|_{0,\epsilon}. $$

In conclusion, according to the triangle inequality, we obtain the following a-priori estimates in terms of the approximation error,

$$ \begin{aligned} \|(u_{\epsilon}^N - u_{\epsilon,h}^N, \lambda_{\epsilon}^N - \lambda_{\epsilon,h}^N)\|_{\epsilon} &\lesssim \|(u_{\epsilon}^N - R_h^k u_{\epsilon}^N, 0)\|_{\epsilon} + \|(R_h^k u_{\epsilon}^N - u_{\epsilon,h}^N, \lambda_{\epsilon}^N - \lambda_{\epsilon,h}^N)\|_{\epsilon}, \\ &\lesssim \|(u_{\epsilon}^N - R_h^k u_{\epsilon}^N, 0)\|_{\epsilon} + \frac{h^{\frac{1}{2}}}{\epsilon} \|\lambda_{\epsilon}^N\|_{0,\epsilon}. \end{aligned} \quad (4.14) $$

The a-priori error estimate (4.14) must be compared with the corresponding results (4.4)-(4.5) obtained for the non stabilized scheme. We remark that (4.14) is not affected by any restriction on $N$ and it does not suffer from locking effects. However, it shows that the stabilized scheme necessarily converges at a sub-optimal rate $h^{\frac{1}{2}}$ for any polynomial order $k > 0$. Then, if low order convergence is acceptable, one should look at the stabilized scheme with low order finite element approximation ($k = 1$). Otherwise, if a sufficiently refined resolution of the inclusions is affordable, namely $h \simeq \epsilon$ and if the constraint $k \ge N$ can be satisfied, the estimate results (4.4)-(4.5) seem more attractive, provided that the finite element space suitably approximates the solution of the problem. This is another difficulty of this problem (see for example Remark 4.3) that will be discussed and resolved in the next two sections.

## 4.3 Approximation properties and convergence

Because of the limited regularity of the solution, the approximation properties of the finite element method are not optimal. In particular, due to the approximation properties (4.3) satisfied by the Scott-Zhang interpolation operator, since $u_{\epsilon}^N \in H_0^s(\Omega)$ with $s = \frac{3}{2} - \eta$, we have the following interpolation error estimate, for any $k > 0$:

$$ \|u_{\epsilon}^N - R_h^k u_{\epsilon}^N\|_{0,\Omega} + h \|\nabla(u_{\epsilon}^N - R_h^k u_{\epsilon}^N)\|_{0,\Omega} \lesssim h^{\frac{1}{2}-\eta} \|u_{\epsilon}^N\|_{\frac{3}{2}-\eta,\Omega}. \quad (4.15) $$

Thus, for the discrete problem without stabilization studied in Section 4.1, if we combine (4.15) with the a-priori error estimates (4.4) and (4.5) we obtain that, under the assumption that $k \ge N$,

$$ \|u_{\epsilon}^N - u_{\epsilon,h}^N\|_{1,\Omega} \lesssim (1 + \frac{h}{\epsilon})^N h^{\frac{1}{2}-\eta} \|u_{\epsilon}^N\|_{\frac{3}{2}-\eta,\Omega}, \quad (4.16) $$

$$ \|\lambda_{\epsilon}^N - \lambda_{\epsilon,h}^N\|_{-\frac{1}{2},\epsilon} \lesssim (1 + \frac{h}{\epsilon})^{2N} h^{\frac{1}{2}-\eta} \|u_{\epsilon}^N\|_{\frac{3}{2}-\eta,\Omega}. \quad (4.17) $$
On the other hand, for the stabilized method studied in Section 4.2, adopting a similar approach, we have for any $k > 0$ the following result

$$ \|(u_{\varepsilon}^N - u_{\varepsilon,h}^N, \lambda_{\varepsilon}^N - \lambda_{\varepsilon,h}^N)\|_{\varepsilon} \le h^{\frac{1}{2}-\eta} (\|u_{\varepsilon}^N\|_{\frac{3}{2}-\eta,\Omega} + \varepsilon^{-1} \|\lambda_{\varepsilon}^N\|_{0,\varepsilon}). \quad (4.18) $$

The previous inequalities show that with respect to the approximation properties, the numerical schemes with or without stabilization are equivalent and both are sub-optimal. We also observe that there is no advantage in using high-order finite elements. As a consequence, since we are limited to the case $k = 1$, the restriction $k \ge N$ for the scheme without stabilization entails that we can not exploit the additional accuracy provided by high-order Fourier modes. In this respect, the stabilized scheme guarantees more flexibility.

**Remark 4.3.** In the limit when $\varepsilon \to 0$, the regularity of the continuous problem further decreases. In particular, the $H^s$-norm of $u_{\varepsilon}^N$ for $s \ge 1$, which pop out in (4.16)-(4.17) and (4.18) may tend to infinity when $\varepsilon \to 0$. This can be illustrated by the very simple exterior problem

$$ \begin{cases} -\Delta \hat{v}_{\varepsilon} = 0 & \text{in } \omega_{\Upsilon} \setminus \omega_{\varepsilon}, \\ \hat{v}_{\varepsilon} = 0 & \text{on } \partial \omega_{\Upsilon}, \\ \hat{v}_{\varepsilon} = L & \text{on } \partial \omega_{\varepsilon}, \end{cases} $$

that features the following analytical solution

$$ \hat{v}_{\varepsilon}(|x|) = L \frac{\ln(|x|/\Upsilon)}{\ln(\varepsilon/\Upsilon)}. $$

Thus, we obtain, for $s \ge 1$,

$$ \|\hat{v}_{\varepsilon}\|_{s,\Omega} \sim \frac{1}{\varepsilon^{s-1}} \quad \text{when } \varepsilon \to 0. $$

These considerations show that for small inclusions it is particularly difficult to achieve a desired accuracy with a standard approximation method.

**Remark 4.4.** Under some additional assumptions on the mesh, we can prove convergence results better than those stated in (4.16)-(4.17) and (4.18). Let $\mathcal{T}_h^\Omega$ be a $\delta$-resolving mesh with respect to the interface, that is the edges of $\mathcal{T}_h^\Omega$ have a maximum distance of $\delta$ to $\partial \omega_\varepsilon$. This condition is in particular fulfilled with $\delta = O(h^2)$ when the nodes of the mesh $\mathcal{T}_h^\Omega$ fall on the interface $\partial \omega_\varepsilon$. With a little abuse of notation, we call this case the conforming mesh configuration. In this case, following the approach of Theorem 5.2 in [26], under the restriction that $k = 1 \ge N$, we obtain the following quasi-optimal convergence results for the numerical scheme without stabilization,

$$ \begin{aligned} \|u_{\varepsilon}^N - u_{\varepsilon,h}^N\|_{1,\Omega} &\lesssim (1 + \frac{h}{\varepsilon})^N h (\|u_{\varepsilon}^N\|_{2,\Omega_{\varepsilon}} + \|u_{\varepsilon}^N\|_{2,\omega_{\varepsilon}}), \\ \|\lambda_{\varepsilon}^N - \lambda_{\varepsilon,h}^N\|_{-\frac{1}{2},\varepsilon} &\lesssim (1 + \frac{h}{\varepsilon})^{2N} h (\|u_{\varepsilon}^N\|_{2,\Omega_{\varepsilon}} + \|u_{\varepsilon}^N\|_{2,\omega_{\varepsilon}}). \end{aligned} $$

We note that for $\delta$-meshes, the stabilized scheme remains sub-optimal because the consistency error, of order $h^{\frac{1}{2}}$, dominates over the approximation error.

Overall, the approximation properties of the proposed (original and stabilized) schemes are not satisfactory, in particular because quadratic convergence rates or higher are prevented. To overcome this limitation, we discuss an improvement of the approximation method in the next section.

## 4.4 Augmented finite element method

As pointed out in the previous section, the classical finite element method may suffer from numerical locking when $\varepsilon/\hbar \to 0$. The other limitation of the classical finite element method is the sub-optimality of the convergence in $h$ for continuous finite elements when the mesh does not conform to the $\partial \omega_\varepsilon$ interface. This is a classical limitation of fictitious domain problems and can be attributed to the poor estimation of the gradient jump by $C^1$ functions across the interface. In this section, we will introduce an augmented finite element method which allows us to obtain a numerical approximation error independent of $\varepsilon$ and of arbitrary order in $h$.
### 4.4.1 Definition of the augmented finite element space

In order to overcome the limitations described above, we propose to enrich the finite element space by adding well-chosen functions that allow to better approximate the singular behaviour of $u_{\varepsilon}^N$ when $\varepsilon$ tends to 0 and the jump at the interface $\partial\omega_{\varepsilon}$. The design of these functions directly comes from the expression of the solution $u_{\varepsilon}^N$ in $\omega_{\Upsilon}$ with respect to its Fourier modes on $\partial\omega_{\Upsilon}$ and $\partial\omega_{\varepsilon}$ given by equation (3.18).

In what follows, we assume that $\varepsilon > 0$ is small enough to have $2\varepsilon < \Upsilon$ where the parameter $\Upsilon$ is defined at section 3.2. We first define $\chi$ a radial cut-off $C^{\infty}$ function satisfying

$$ \chi(r) = \begin{cases} 1 & \text{if } r \le \frac{1}{2}\Upsilon, \\ 0 & \text{if } r > \frac{3}{4}\Upsilon. \end{cases} $$

Then, we set, for $r > 0$ and $\theta \in [0, 2\pi[$,

$$ \phi_{\varepsilon}^{0}(r, \theta) = \chi(r) \frac{\ln(r/\Upsilon)}{\ln(\varepsilon/\Upsilon)} \mathbb{1}_{r \ge \varepsilon} + \mathbb{1}_{r < \varepsilon} $$

and, for $1 \le n \le N$,

$$ \left\{ \begin{aligned} \phi_{\varepsilon}^{n,0}(r, \theta) &:= \chi(r) \left( \frac{\left(\frac{\Upsilon}{r}\right)^n - \left(\frac{r}{\Upsilon}\right)^n}{\left(\frac{\Upsilon}{\varepsilon}\right)^n - \left(\frac{r}{\Upsilon}\right)^n} \right) \cos(n\theta) \mathbb{1}_{r \ge \varepsilon} + \left(\frac{r}{\varepsilon}\right)^n \cos(n\theta) \mathbb{1}_{r < \varepsilon}, \\ \phi_{\varepsilon}^{n,1}(r, \theta) &:= \chi(r) \left( \frac{\left(\frac{\Upsilon}{r}\right)^n - \left(\frac{r}{\Upsilon}\right)^n}{\left(\frac{\Upsilon}{\varepsilon}\right)^n - \left(\frac{r}{\Upsilon}\right)^n} \right) \sin(n\theta) \mathbb{1}_{r \ge \varepsilon} + \left(\frac{r}{\varepsilon}\right)^n \sin(n\theta) \mathbb{1}_{r < \varepsilon}. \end{aligned} \right. $$

We consider the augmented finite element space $X_h^{k,N}(\Omega)$ equal to

$$ X_h^{k,N}(\Omega) = X_h^k(\Omega) \cup \operatorname{span}\{\phi_{\varepsilon}^0, \phi_{\varepsilon}^{n,j}, 1 \le n \le N, 0 \le j \le 1\}. $$

In particular, we notice that the number of additional functions is equal to $2N + 1$, that is the dimension of $\mathcal{M}_{\varepsilon}^N$. Indeed, as we will see in the proof of Theorem 4.6, adding a Lagrange multiplier allows to better solve the condition at the interface but also adds a term in the singular part when $\varepsilon \to 0$ and in the gradient jump at the $\partial\omega_{\varepsilon}$ interface.

We can now introduce the discrete formulation: find $(\tilde{u}_{\varepsilon,h}^N, \tilde{\lambda}_{\varepsilon,h}^N) \in X_h^{k,N}(\Omega) \times \mathcal{M}_{\varepsilon}^N$ such that

$$ \left\{ \begin{aligned} & (\nabla \tilde{u}_{\varepsilon,h}^N, \nabla v_h)_\Omega + (\tilde{\lambda}_{\varepsilon,h}^N, v_h)_\varepsilon = (f, v_h)_\Omega, && \forall v_h \in X_h^{k,N}(\Omega), \\ & (\mu, \tilde{u}_{\varepsilon,h}^N)_\varepsilon = (\mu, \varphi)_\varepsilon, && \forall \mu \in \mathcal{M}_{\varepsilon}^N. \end{aligned} \right. \quad (4.19) $$

### 4.4.2 Inf-sup stability

In the augmented finite element space, we are able to derive an inf-sup condition uniformly independent of $h$, $N$ and $\varepsilon$.

**Lemma 4.5.** There exists $\mathbb{R} \ni \beta > 0$ uniformly independent of $h$, $N$ and $\varepsilon$ such that for all $\lambda \in \mathcal{M}_{\varepsilon}^N$,

$$ \sup_{v_h \in X_h^{k,N}(\Omega)} \frac{(\lambda, v_h)_\varepsilon}{\|v_h\|_{1,\Omega}} \ge \beta \|\lambda\|_{-\frac{1}{2},\varepsilon}. $$

*Proof.* For all $1 \le n \le N$ and $j \in \{0, 1\}$, we define $\phi_{\varepsilon,z}^{n,j}$ as the solution of the equation (3.3) with

$$ z = \cos(n\theta) \text{ if } j = 0 \text{ or } z = \sin(n\theta) \text{ if } j = 1. $$

Then, we notice that $\phi_{\varepsilon}^{n,j} = \chi(r)\phi_{\varepsilon,z}^{n,j}$. Similarly, we define $\phi_{\varepsilon,z}^0$ as the solution of the equation (3.3) with $z = 1$. Then, we notice that $\phi_{\varepsilon}^0 = \chi(r)\phi_{\varepsilon,z}^0$.

For all $\lambda \in \mathcal{M}_{\varepsilon}^N$, we have

$$ \|\lambda\|_{-\frac{1}{2},\varepsilon} = \sup_{\zeta \in H^{\frac{1}{2}}(\partial\omega_{\varepsilon})} \frac{(\lambda, \zeta)_{\varepsilon}}{\|\zeta\|_{\frac{1}{2},\varepsilon}}. $$
For a given $\zeta \in H^{\frac{1}{2}}(\partial\omega_\epsilon)$, we define $v_\zeta^z$ and $v_h \in H_0^1(\Omega)$ by

$$v_{\zeta}^{z} = \mathcal{A}_{\epsilon}^{0}(\zeta)\phi_{\epsilon,z}^{0} + \sum_{n=1}^{N} \mathcal{A}_{\epsilon}^{n}(\zeta)\phi_{\epsilon,z}^{n,0} + \mathcal{B}_{\epsilon}^{n}(\zeta)\phi_{\epsilon,z}^{n,1} \quad \text{and} \quad v_h = \chi(r)v_{\zeta}^{z}.$$

We have by linearity that $v_\zeta^z$ satisfies equation (3.3) with $z = \Pi_\epsilon^N(\zeta)$ and according to Lemma 3.4 $\|v_\zeta^z\|_{1,\Omega} \lesssim \|\Pi_\epsilon^N(\zeta)\|_{\frac{1}{2},\epsilon} \lesssim \|\zeta\|_{\frac{1}{2},\epsilon}$. Moreover, since $\chi(r) = 1$ in $\partial\omega_\epsilon$, we have $v_h = v_\zeta^z = \Pi_\epsilon^N(\zeta)$ on $\partial\omega_\epsilon$. Thus $\|v_h\|_{1,\Omega} \lesssim \|v_\zeta^z\|_{1,\Omega} \lesssim \|\zeta\|_{\frac{1}{2},\epsilon}$, that is $\|v_h\|_{1,\Omega} \le C_0\|\zeta\|_{\frac{1}{2},\epsilon}$ where $C_0$ is a generic constant independent of $h, N$ and $\epsilon$. At last, by definition of $X_h^{k,N}(\Omega)$, $v_h \in X_h^{k,N}(\Omega)$. Thus, we obtain that, for all $\zeta \in H^{\frac{1}{2}}(\partial\omega_\epsilon)$

$$\frac{(\lambda, \zeta)_\epsilon}{\|\zeta\|_{\frac{1}{2},\epsilon}} = \frac{(\lambda, \Pi_\epsilon^N(\zeta))_\epsilon}{\|\zeta\|_{\frac{1}{2},\epsilon}} \le C_0 \sup_{v_h \in X_h^k(\Omega)} \frac{(\lambda, v_h)_\epsilon}{\|v_h\|_{1,\Omega}}.$$

By taking the supremum over $\zeta \in H^{\frac{1}{2}}(\partial\omega_\epsilon)$, this gives the result with $\beta = \frac{1}{C_0}$. $\square$

### 4.4.3 Approximation properties and convergence

To conclude this section, we give an estimate of the approximation error for the augmented finite element method.

**Theorem 4.6.** Let $k \in \mathbb{N}^*$. We assume that $\Omega$ is of class $C^{k+1}$ and that $f$ belongs to $H^{k-1}(\Omega)$. Let $(u_\epsilon^N, \lambda_\epsilon^N)$ be the solution of problem (2.16) and $(\tilde{u}_\epsilon^N, \tilde{\lambda}_\epsilon^N)$ be the solution of problem (4.19). Then, we have

$$\|u_{\epsilon}^{N} - \tilde{u}_{\epsilon,h}^{N}\|_{1,\Omega} + \|\lambda_{\epsilon}^{N} - \tilde{\lambda}_{\epsilon,h}^{N}\|_{-\frac{1}{2},\epsilon} \lesssim h^{k} \left( \|f\|_{k-1,\Omega} + \|\varphi\|_{\frac{1}{2},\epsilon} \right). \quad (4.20)$$

Let us mention that, up to now, we have assumed that the domain is polygonal in order to have meshes that match the boundary. In the theorem, we assume that the domain is regular and we neglect the error due to the approximation of this regular domain by a polygonal domain.

If we compare the result given by our theorem to the approximation properties of the original and stabilized schemes presented in Section 4.3, we see that the augmented finite element method overcomes the previous limitations of suboptimal convergence. In addition, if we gather the estimate (3.8) of the error $u_\epsilon - u_\epsilon^N$ coming from the model reduction and the estimate (4.20) of the error $u_\epsilon^N - \tilde{u}_\epsilon^N$ coming from the numerical approximation, we get this estimate for the global error:

$$\|u_{\epsilon} - \tilde{u}_{\epsilon,h}^{N}\|_{1,\Omega} \lesssim (\epsilon^{N+1} + h^k) \left( \|f\|_{k-1,\Omega} + \|\varphi\|_{\frac{1}{2},\epsilon} \right).$$

This estimate shows that, $h$ and $\epsilon$ being fixed, it is possible to reach a given accuracy by adjusting the parameters of the method $k$ and $N$.

*Proof of Theorem 4.6.* Since, for all $n \ge 0$ and $j = 0, 1$, $\phi_\epsilon^{n,j}$ belongs to $H_0^1(\Omega)$, the enrichment of the space $X_h^k(\Omega)$ preserves the coercivity condition. Moreover, according to Lemma 4.5, the inf-sup condition is satisfied for the spaces $X_h^{k,N}(\Omega)$ and $\mathcal{M}_\epsilon^N$. So, we have

$$\|u_{\epsilon}^{N} - \tilde{u}_{\epsilon,h}^{N}\|_{1,\Omega} + \|\lambda_{\epsilon}^{N} - \tilde{\lambda}_{\epsilon,h}^{N}\|_{-\frac{1}{2},\epsilon} \lesssim \inf_{w_h \in X_h^{k,N}(\Omega)} \|u_{\epsilon}^{N} - w_h\|_{1,\Omega}. \quad (4.21)$$

As mentioned previously, the enrichment of the classical finite element space $X_h^k(\Omega)$ by $X_h^{k,N}(\Omega)$ will allow to get a better approximation of $u_\epsilon^N$ by providing a better estimate of the singular part of $u_\epsilon^N$. To highlight this feature, we use (3.18) and write $u_\epsilon^N$ as the sum of a regular part $u_{\epsilon,R}^N$ and a singular part $u_{\epsilon,S}^N$ defined respectively by: for $0 < r < \Upsilon$ and $\theta \in [0, 2\pi[$, if $\mathbf{x} = (\cos(\theta), \sin(\theta)) \in \partial\omega$,

$$u_{\epsilon,R}^{N}(\mathbf{r}\mathbf{x}) := \mathcal{A}_{\Upsilon}^{0}(u_{\epsilon}^{N}) + \sum_{n \in \mathbb{N}^{*}} \left(\frac{r}{\Upsilon}\right)^{n} \left(\mathcal{A}_{\Upsilon}^{n}(u_{\epsilon}^{N}) \cos(n\theta) + \mathcal{B}_{\Upsilon}^{n}(u_{\epsilon}^{N}) \sin(n\theta)\right)$$
and

$$
\begin{aligned}
u_{\varepsilon,S}^N(\mathbf{r}\mathbf{x}) := & \left( -\mathcal{A}_{\Upsilon}^0(u_{\varepsilon}^N) \frac{\ln(r/\Upsilon)}{\ln(\varepsilon/\Upsilon)} - \sum_{1 \le n \le N} \left(\frac{\varepsilon}{\Upsilon}\right)^n \frac{\left(\frac{r}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{r}\right)^n}{\left(\frac{\varepsilon}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{\varepsilon}\right)^n} (\mathcal{A}_{\Upsilon}^n(u_{\varepsilon}^N) \cos(n\theta) + \mathcal{B}_{\Upsilon}^n(u_{\varepsilon}^N) \sin(n\theta)) \right. \\
& \left. + \mathcal{A}_{\varepsilon}^0(\varphi) \frac{\ln(r/\Upsilon)}{\ln(\varepsilon/\Upsilon)} + \sum_{1 \le n \le N} \left(\frac{\frac{r}{\Upsilon}}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{r}\right)^n - \left(\frac{\varepsilon}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{\varepsilon}\right)^n \left(\mathcal{A}_{\varepsilon}^n(\varphi) \cos(n\theta) + \mathcal{B}_{\varepsilon}^n(\varphi) \sin(n\theta)\right) \right) \mathbf{1}_{r \ge \varepsilon} \\
& + \left( -\mathcal{A}_{\Upsilon}^0(u_{\varepsilon}^N) - \sum_{1 \le n \le N} \left(\frac{r}{\Upsilon}\right)^n \left(\mathcal{A}_{\Upsilon}^n(u_{\varepsilon}^N) \cos(n\theta) + \mathcal{B}_{\Upsilon}^n(u_{\varepsilon}^N) \sin(n\theta)\right) \right. \\
& \left. + \mathcal{A}_{\varepsilon}^0(\varphi) + \sum_{1 \le n \le N} \left(\frac{r}{\varepsilon}\right)^n \left(\mathcal{A}_{\varepsilon}^n(\varphi) \cos(n\theta) + \mathcal{B}_{\varepsilon}^n(\varphi) \sin(n\theta)\right) \right) \mathbf{1}_{r < \varepsilon} \right)
\end{aligned}
$$

The functions $u_{\varepsilon,R}^N$ and $u_{\varepsilon,S}^N$ belong to $H^1(\Omega)$ and, according to (3.18), we have $u_{\varepsilon}^N = u_{\varepsilon,R}^N + u_{\varepsilon,S}^N$ in $\omega_{\Upsilon}$, which leads in particular to

$$ \chi u_{\varepsilon}^{N} = \chi u_{\varepsilon,R}^{N} + \chi u_{\varepsilon,S}^{N}. $$

Let us now prove that the truncated singular part belongs to the approximation space, that is that $\chi u_{\varepsilon,S}^N$ belongs to $X_h^{k,N}(\Omega)$. First, we notice that, by definition of $\phi_{\varepsilon}^0$, we have, for all $r > 0$ and $\theta \in [0, 2\pi[$,

$$ \chi(r) (-\mathcal{A}_{\Upsilon}^{0}(u_{\varepsilon}^{N}) + \mathcal{A}_{\varepsilon}^{0}(\varphi)) \left( \frac{\ln(r/\Upsilon)}{\ln(\varepsilon/\Upsilon)} \mathbf{1}_{r \ge \varepsilon} + \mathbf{1}_{r < \varepsilon} \right) = (-\mathcal{A}_{\Upsilon}^{0}(u_{\varepsilon}^{N}) + \mathcal{A}_{\varepsilon}^{0}(\varphi)) \phi_{\varepsilon}^{0}(r, \theta) $$

Moreover, by definition of $\phi_{\varepsilon}^{n,j}$, we have, for all $1 \le n \le N$

$$
\begin{aligned}
& \chi(r) \left( -\mathcal{A}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{\varepsilon}{\Upsilon}\right)^n \frac{\left(\frac{r}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{r}\right)^n}{\left(\frac{\varepsilon}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{\varepsilon}\right)^n} + \mathcal{A}_{\varepsilon}^n(\varphi) \frac{\left(\frac{r}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{r}\right)^n}{\left(\frac{\varepsilon}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{\varepsilon}\right)^n} \right) \cos(n\theta) \mathbf{1}_{r \ge \varepsilon} \right. \\
& \left. + \chi(r) \left( -\mathcal{A}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{r}{\Upsilon}\right)^n + \mathcal{A}_{\varepsilon}^n(\varphi) \left(\frac{r}{\varepsilon}\right)^n \right) \cos(n\theta) \mathbf{1}_{r < \varepsilon} = \left( -\mathcal{A}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{\varepsilon}{\Upsilon}\right)^n + \mathcal{A}_{\varepsilon}^n(\varphi) \right) \phi_{\varepsilon}^{n,0} \right)
\end{aligned}
$$

and

$$
\begin{aligned}
& \chi(r) \left( -\mathcal{B}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{\varepsilon}{\Upsilon}\right)^n \frac{\left(\frac{r}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{r}\right)^n}{\left(\frac{\varepsilon}{\Upsilon}\right)^n - \left(\frac{\Upsilon}{\varepsilon}\right)^n} + \mathcal{B}_{\Upsilon}^n(\varphi) \left(\frac{r}{\Upsilon}\right)^n - \left(\frac{r}{\varepsilon}\right)^n \right) \sin(n\theta) \mathbf{1}_{r \ge \varepsilon} \\
& + \chi(r) \left( -\mathcal{B}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{r}{\Upsilon}\right)^n + \mathcal{B}_{\varepsilon}^n(\varphi) \left(\frac{r}{\varepsilon}\right)^n \right) \sin(n\theta) \mathbf{1}_{r < \varepsilon} = \left( -\mathcal{B}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{\varepsilon}{\Upsilon}\right)^n + \mathcal{B}_{\varepsilon}^n(\varphi) \right) \phi_{\varepsilon}^{n,1}.
\end{aligned}
$$

We deduce that for $r > 0$ and $\theta \in [0, 2\pi[$,

$$
\begin{aligned}
\chi(r) u_{\varepsilon,S}^N(r, \theta) = & (-\mathcal{A}_{\Upsilon}^0(u_{\varepsilon}^N) + \mathcal{A}_{\varepsilon}^0(\varphi)) \phi_{\varepsilon}^0(r, \theta) + \sum_{1 \le n \le N} \left( -\mathcal{A}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{\varepsilon}{\Upsilon}\right)^n + \mathcal{A}_{\varepsilon}^n(\varphi) \right) \phi_{\varepsilon}^{n,0}(r, \theta) \\
& + \sum_{1 \le n \le N} \left( -\mathcal{B}_{\Upsilon}^n(u_{\varepsilon}^N) \left(\frac{\varepsilon}{\Upsilon}\right)^n + \mathcal{B}_{\varepsilon}^n(\varphi) \right) \phi_{\varepsilon}^{n,1}(r, \theta),
\end{aligned}
$$

and thus we get that $\chi u_{\varepsilon,S}^N \in X_h^{k,N}(\Omega)$. Therefore, for the right hand side of (4.21), we have in particular

$$ \inf_{w_h \in X_h^{k,N}(\Omega)} \|u_{\varepsilon}^N - w_h\|_{1,\Omega} \le \inf_{v_h \in X_h^k(\Omega)} \|u_{\varepsilon}^N - \chi u_{\varepsilon,S}^N - v_h\|_{1,\Omega}. $$

Moreover, writing $u_{\varepsilon}^N = \chi u_{\varepsilon}^N + (1 - \chi)u_{\varepsilon}^N$, it comes

$$ u_{\varepsilon}^N - \chi u_{\varepsilon,S}^N = \chi u_{\varepsilon,R}^N + (1 - \chi)u_{\varepsilon}^N, $$

so that estimate (4.21) becomes

$$ \|u_{\varepsilon}^{N} - \tilde{u}_{\varepsilon,h}^{N}\|_{1,\Omega} + \|\lambda_{\varepsilon}^{N} - \tilde{\lambda}_{\varepsilon,h}^{N}\|_{-\frac{1}{2},\varepsilon} \lesssim \inf_{v_h \in X_h^k(\Omega)} \|\chi u_{\varepsilon,R}^{N} + (1-\chi)u_{\varepsilon}^{N} - v_h\|_{1,\Omega}. \quad (4.22) $$

It is now sufficient to prove that the function $\chi u_{\varepsilon,R}^N + (1 - \chi)u_{\varepsilon}^N$ is well estimated by the elements of $X_h^k(\Omega)$. To do this, we will prove that this function belongs to $H^{k+1}(\Omega)$.
First, since $u_{\varepsilon,R}^N$ is an harmonic function, thanks to the interior regularizing property of elliptic operators (we refer to [14]), we have

$$ \|u_{\varepsilon,R}^N\|_{H^{k+1}(\omega_{3\Upsilon/4})}^2 \lesssim \|u_{\varepsilon,R}^N\|_{L^2(\omega_{\Upsilon})}^2 \lesssim |A_{\Upsilon}^0(u_{\varepsilon}^N)|^2 + \sum_{n \in \mathbb{N}^*} |A_{\Upsilon}^n(u_{\varepsilon}^N)|^2 + |B_{\Upsilon}^n(u_{\varepsilon}^N)|^2 \\ \lesssim \|u_{\varepsilon}^N\|_{L^2(\omega_{\Upsilon})} \lesssim \|f\|_{0,\Omega} + \|\varphi\|_{\frac{1}{2},\varepsilon} $$

thanks to Theorem 3.1. Thus, by definition of $\chi$, we get that $\chi u_{\varepsilon,R}^N$ belongs to $H^{k+1}(\Omega)$ and we have

$$ \|\chi u_{\varepsilon,R}^N\|_{k+1,\Omega} \lesssim \|f\|_{0,\Omega} + \|\varphi\|_{\frac{1}{2},\varepsilon} $$

Moreover, since $u_{\varepsilon}^N$ satisfies $-\Delta u_{\varepsilon}^N = f$ in $\Omega_{\varepsilon}$ and $u_{\varepsilon}^N = 0$ on $\partial\Omega$ and since $\Omega \setminus \omega_{\Upsilon/2}$ is relatively compact in $\Omega_{\varepsilon} \cup \partial\Omega$, the regularizing property of elliptic problems (we also refer to [14] for this result) implies that $u_{\varepsilon}^N$ belongs to $H^{k+1}(\Omega \setminus \omega_{\Upsilon/2})$ and we have

$$ \|u_{\varepsilon}^N\|_{H^{k+1}(\Omega \setminus \omega_{\Upsilon/2})} \lesssim \|f\|_{k-1,\Omega}. $$

This implies that $(1 - \chi)u_{\varepsilon}^N$ belongs to $H^{k+1}(\Omega)$ and we have

$$ \|(1 - \chi)u_{\varepsilon}^N\|_{k+1,\Omega} \lesssim \|f\|_{k-1,\Omega}. $$

Gathering these results, we deduce that $\chi u_{\varepsilon,R}^N + (1 - \chi)u_{\varepsilon}^N$ belongs to $H^{k+1}(\Omega)$ and

$$ \|\chi u_{\varepsilon,R}^N + (1 - \chi)u_{\varepsilon}^N\|_{k+1,\Omega} \lesssim \|f\|_{k-1,\Omega} + \|\varphi\|_{\frac{1}{2},\varepsilon}. $$

By this way, using interpolation estimates for the finite element space $X_h^k(\Omega)$ in $H^{k+1}(\Omega)$ [13], (4.22) becomes

$$ \|u_{\varepsilon}^N - \tilde{u}_{\varepsilon,h}^N\|_{1,\Omega} + \|\lambda_{\varepsilon}^N - \tilde{\lambda}_{\varepsilon,h}^N\|_{-\frac{1}{2},\varepsilon} \lesssim Ch^k \|\chi u_{\varepsilon,R}^N + (1 - \chi)u_{\varepsilon}^N\|_{k+1,\Omega} $$

and thus we get (4.20).

□

**Remark 4.7.** Let us notice that the use of an enriched element space allows us to derive an error estimate that is at the same time independent of $\varepsilon$ and can be of any order in $h$. This is first made possible by getting an inf-sup condition independent of $\varepsilon$ and $h$, but also by the fact that the singular part of the solution $u_{\varepsilon}^N$ is only in the first Fourier modes, and is therefore fully captured by the enrichment functions.

Eventually, we see that the estimate of Theorem 3.1 requires only a regularity in $H^{\frac{1}{2}}(\partial\omega_{\varepsilon})$ for $\varphi$. This assumption, which may seem rather weak, is due to the fact that thanks to the enrichment of the finite element space the convergence analysis of the scheme only relies on the regularity of the solution away from the inclusion.

# 5 Numerical experiments

Following the main results of this work, we subdivide this section in two parts. First, we analyze the behavior of the error due to the defective boundary conditions, referring to this part as the convergence in $\varepsilon$. We perform numerical tests to illustrate that the properties stated in Theorem 3.5 are actually observed in practice. Second, we analyze the approximation properties of the reduced problem discretized by the finite element method, using both the standard formulation and the augmented formulation. The numerical tests shown in the second section address the error analysis with respect to the space discretization parameter $h$ and complement the analysis shown in Section 4.

## 5.1 Convergence in $\varepsilon$

In this section we analyze first the case of a single inclusion, for which a rigorous analysis has been presented here. However, as mentioned in the introduction, the proposed approach can be easily extended to several obstacles of different sizes. The behavior of our method in this case is addressed in the second paragraph of this section. Finally, in the third part of this section, for the case of multiple inclusions, we analyze the effect of the distance between inclusions, which is an important factor in many applications.
### 5.1.1 Single inclusion

We consider a rectangular domain $\Omega$ of width $L$ and height $l$ with an inclusion located in $z$ (we have taken $L = l = 1$ and $z = (0.2, 0.1)$) and we make $\varepsilon$ vary in the set $\{0.05, 0.04, 0.03, 0.02, 0.01\}$. On the external boundary $\partial\Omega$, we have considered non homogeneous Dirichlet boundary conditions presented in Figure 2. Although our mathematical study focused on homogeneous Dirichlet boundary conditions, the results of the previous sections remain valid in this more general case. As for the inclusion, we set $\varphi = 0$. To illustrate Theorem 3.5, we have computed the convergence rate with respect to $\varepsilon$ of the reduced problem towards the full problem for increasing values of $N$ in the definition of $\mathcal{M}_{\varepsilon}^N$. In the general case, the continuous solutions are not known a priori. Therefore, we have considered the linear finite element approximation of the reduced and full problems for a conforming mesh (defined as $\delta$-resolving mesh in Remark 4.4) with $h$ small enough ($h \approx 1e-3$) in order to have a negligible discretization error in comparison with the model reduction error. The convergence curves are given in Figure 3. They represent the $H^1$-norms of the difference $e_{\varepsilon,h}^N$ between the discrete solutions of the full problem and the reduced problem for $N = 0, 1$ and $2$ moments. If we compare the slopes of the error curves with the expected slopes (represented by the black lines), we notice that the match is between theory and observations is very satisfactory. The numerical solution of the reduced problem and the errors $e_{\varepsilon,h}^0, e_{\varepsilon,h}^1$ and $e_{\varepsilon,h}^2$ for $\varepsilon = 0.05$ are reported in Figure 4.

<figure><img src="image_4.png" /><figcaption>Figure 2: Boundary conditions for the square domain.</figcaption></figure>

### 5.1.2 Multiple inclusions

We consider a rectangular domain $\Omega = [-0.5, 0.5] \times [-0.5, 0.5]$ with four inclusions positioned respectively in $(0.2, 0.1)$, $(-0.3, 0.1)$, $(-0.2, 0.2)$, $(0.1, -0.2)$ and of respective initial sizes 0.05, 0.04, 0.06 and 0.03. As in the previous tests, the boundary conditions on $\partial\Omega$ are described by Figure 2 and $\varphi = 0$ on the inclusions. The results obtained for $N = 0, 1$ and $2$ are reported in Figure 5. There, the parameter $\varepsilon$ corresponds to the size of the first inclusion and all the other inclusions are scaled proportionally. Again, we can notice that the convergence rates coincide with the expected ones. The numerical solution of the reduced problem and the errors $e_{\varepsilon,h}^0, e_{\varepsilon,h}^1$ and $e_{\varepsilon,h}^2$ are reported in Figure 6.

### 5.1.3 Behavior of the model error for close obstacles

We consider here a test case to illustrate the behavior of the method when obstacles are close to each other. The domain is the same as previously and we consider uniform Dirichlet boundary conditions equal to 1 and two inclusions of radius 0.05 initially centered in $(-0.1, 0)$ and $(0.1, 0)$. The boundary values imposed on each inclusion are respectively 0.5 and 1.5 and are enforced with $N \in \{0, 1, 2\}$ modes. The numerical solution of such problem is displayed in Figure 7 panel (a). Then, the distance between the two inclusions is gradually reduced. The numerical solution with inclusions separated by a distance of 0.01 is shown in panel (b). The comparison of the two top panels shows that the solution in the circle
<figure><img src="image_5.png" /><figcaption>Figure 3: The modelling error behavior with respect to the radius for different values of N in the case of a single inclusion.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Contour plots of the solution $u_{\epsilon}$ and of the discrete model errors $e_{\epsilon,h}^0$, $e_{\epsilon,h}^1$ and $e_{\epsilon,h}^2$ on a log-scale axis for an inclusion of radius $\epsilon = 0.05$.</figcaption></figure>

deviates from the constant as the two inclusions get closer. More quantitative results are given in Figure 7 panel (c). We see that the model error increases as the distance between the two obstacles decreases, as expected with the inverse dependence on $\Upsilon$ in the estimate stated in Theorem 3.5. We also notice
<figure><img src="image_7.png" /><figcaption>Figure 5: The modelling error with respect to the radius of four inclusions for different values of N.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Contour plots of the solution $u_{\epsilon}$ and of the discrete model errors $e_{\epsilon,h}^0$, $e_{\epsilon,h}^1$ and $e_{\epsilon,h}^2$ on a log-scale axis for multiple obstacles of radius $\epsilon = 0.05$.</figcaption></figure>

that this effect becomes more severe as the number of modes increases.
<figure><img src="image_9.png" /><figcaption>Figure 7: Effect of the inter-inclusions distance on the model error. By comparing panels (a) and (b) we see that the approximation of the solution at the interior boundary worsens when the centers of the inclusions approach. In panel (c) we see the variation of the model error when the distance between the inclusions progressively changes.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>(a) $H^1$ rel. err. $u_{\epsilon,h}^N$, $N=0$, non-augmented FEM.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) $H^1$ rel. err. $u_{\epsilon,h}^N$, $N=0$, augmented FEM.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e) $L^2$ rel. err. $\lambda_{\epsilon,h}^N$, $N=0$, non-augmented FEM.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(g) $L^2$ rel. err. $\lambda_{\epsilon,h}^N$, $N=0$, augmented FEM.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) $H^1$ rel. err. $u_{\epsilon,h}^N$, $N=1$, non-augmented FEM.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) $H^1$ rel. err. $u_{\epsilon,h}^N$, $N=1$, augmented FEM.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f) $L^2$ rel. err. $\lambda_{\epsilon,h}^N$, $N=1$, non-augmented FEM.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(h) $L^2$ rel. err. $\lambda_{\epsilon,h}^N$, $N=1$, augmented FEM.</figcaption></figure>

<caption>Figure 8: Comparison of the numerical error for the finite element method without and with enrichment.</caption>
## 5.2 Convergence in h

We now study the behaviour of the discretization error when the mesh size decreases in order to illustrate the convergence results given by (4.16) and (4.17) for the classical finite element method without stabilization and by Theorem 4.6 for the augmented finite element method.

Here, we have considered $\Omega = [-1.5, 1.5] \times [1.5, 1.5]$ and a single inclusion of center $z = (0, 0)$ and radius $\epsilon \in \{0.1, 0.01, 0.001\}$. For the number of modes, we have taken $N = 0$ or $1$. Moreover, we consider non homogeneous boundary conditions on $\partial\omega_{\epsilon}$. More precisely, denoting by $u_{\epsilon}^e$ the continuous solution, we take

$$ (\mu, u_{\epsilon}^{e})_{\epsilon} = (\mu, \varphi)_{\epsilon}, \forall \mu \in \mathcal{M}_{\epsilon}^{N}, \text{ with } \varphi = 1 + \cos(\theta) + \sin(\theta). $$

By setting the appropriate Dirichlet boundary conditions, the continuous solution writes for $N = 0$,

$$ \left\{ \begin{array}{ll} u_{\epsilon}^{e} = & \begin{cases} 1 + \left(\frac{r}{\rho}\right) \cos(\theta) + \left(\frac{r}{\rho}\right) \sin(\theta) & \text{if } r \le \epsilon, \\ \frac{\ln(r/\rho)}{\ln(\epsilon/\rho)} + \left(\frac{r}{\rho}\right) \cos(\theta) + \left(\frac{r}{\rho}\right) \sin(\theta) & \text{if } r > \epsilon, \end{cases} \\ \\ \lambda_{\epsilon}^{e} = & \frac{1}{\epsilon \ln(\epsilon/\rho)}, \end{array} \right. $$

and for $N = 1$,

$$ \left\{ \begin{array}{ll} u_{\epsilon}^{e} = & \begin{cases} 1 + \left(\frac{r}{\epsilon}\right) \cos(\theta) + \left(\frac{r}{\epsilon}\right) \sin(\theta) + \left(\frac{r}{\rho}\right)^2 \cos(2\theta) + \left(\frac{r}{\rho}\right)^2 \sin(2\theta) & \text{if } r \le \epsilon, \\ \frac{\ln(r/\rho)}{\ln(\epsilon/\rho)} + \frac{(r/\rho) - (\rho/r)}{(\epsilon/\rho) - (\rho/\epsilon)} \cos(\theta) + \frac{(r/\rho) - (\rho/r)}{(\epsilon/\rho) - (\rho/\epsilon)} \sin(\theta) + \left(\frac{r}{\rho}\right)^2 \cos(2\theta) + \left(\frac{r}{\rho}\right)^2 \sin(2\theta) & \text{if } r > \epsilon, \end{cases} \\ \\ \lambda_{\epsilon}^{e} = & \frac{1}{\epsilon \ln(\epsilon/\rho)} + \left( \frac{1/\rho - \rho/\epsilon^2}{\epsilon/\rho - \rho/\epsilon} - \frac{1}{\epsilon} \right) \cos(\theta) + \left( \frac{1/\rho - \rho/\epsilon^2}{\epsilon/\rho - \rho/\epsilon} - \frac{1}{\epsilon} \right) \sin(\theta), \end{array} \right. $$

with $\rho = 0.4$. To solve this problem, we have implemented the classical finite element method given by the formulation (4.1) and the augmented finite element method given by the formulation (4.19). In what follows we analyze the relative error for the solution $u_{\epsilon,h}^N$ in the $H^1$ and $L^2$ norms and for the Lagrange multiplier $\lambda_{\epsilon,h}^N$ in the $L^2$ norm, namely

$$ \frac{\|u_{\epsilon}^{N} - u_{\epsilon,h}^{N}\|_{H^{1}(\Omega)}}{\|u_{\epsilon}^{N}\|_{H^{1}(\Omega)}}, \quad \frac{\|u_{\epsilon}^{N} - u_{\epsilon,h}^{N}\|_{L^{2}(\Omega)}}{\|u_{\epsilon}^{N}\|_{L^{2}(\Omega)}}, \quad \frac{\|\lambda_{\epsilon}^{N} - \lambda_{\epsilon,h}^{N}\|_{L^{2}(\partial\omega_{\epsilon})}}{\|\lambda_{\epsilon}^{N}\|_{L^{2}(\partial\omega_{\epsilon})}}. $$

The numerical results for the h-convergence of the method without and with enrichment functions are compared in Figure 8. We focus our attention on two items: the convergence rate with respect to $h$ and the influence of $\epsilon$ on the error, as in some cases the error estimates show the dependence on $(1 + h/\epsilon)^N$. For the first item, we test the method for linear (FEM-P1) and quadratic (FEM-P2) finite elements (corresponding to different Markers in each plot). For the second item, we exploit the fact that the exact solutions are parametrized with respect to $\epsilon$ and we calculate the error for $\epsilon = 10^{-1}, 10^{-2}, 10^{-3}$. Again, each test is identified by a different marker in the plots (as shown in the legend). Also, in both cases we perform the error analysis for $N = 0$ (on the left column) and $N = 1$ (on the right column).

For the standard finite element formulation, we observe that independently of the choice of finite elements order, linear or quadratic, the convergence rate in the $H^1$-norm is limited to 0.5 because of the low regularity of the solution on the whole domain, for any the value of $N$. We also observe that the scatter of the error curves with respect to $\epsilon$ is more pronounced in the case $N = 1$, as expected because of the factor $(1 + h/\epsilon)^N$ of (4.16). According to the estimate, we see that the error increases when decreasing $\epsilon$ for a fixed value of $h$. This behavior, predicted by the estimate (4.17) for the Lagrange multipliers in the $\| \cdot \|_{-1/2, \epsilon}$-norm is only partially observed in the numerical experiments, as the effect of $\epsilon$ increases moving from $\tilde{N} = 0$ to $N = 1$. We justify this discrepancy observing that we are not calculating the same norm of (4.17). Surprisingly, but not in contrast with the theory, the convergence rate for the Lagrange multiplier increases when increasing the polynomial order of the finite element method. This super-convergence behavior may be due to the particular properties of the exact solution.

For the augmented finite element method, the situation changes completely, in fact the theoretical order of convergence predicted by (4.20) is observed in all the numerical tests. For the $H^1$-norm of the solution, the theoretical order of convergence is almost exactly verified in the experiments. For the Lagrange multipliers, super-convergence is clearly observed for linear finite elements, but the theoretical rate is respected by the quadratic ones. Reminding again that we are not testing the same norm of
the error estimate (but we are actually using a stronger norm), this result may be due to the particular structure of the exact solution. Finally, we observe that for the augmented finite element method the convergence error of the solution and of the Lagrange multipliers is not negatively affected by the size of the inclusion, namely $\epsilon$. In fact, we observe a slight influence of $\epsilon$ on the error, but it decreases almost proportionally with $\epsilon$. Indeed, this is a desirable feature of this numerical approach.

## Conclusions

In this work we have considered a second order PDE defined on a domain with small circular holes subject to non homogeneous Dirichlet boundary conditions. This problem can be seen as a template for different families of applications. It may represent solid mechanics problem where the holes play the role of small supports of a mechanical part, but it may also be used as a prototype problem for fluid-structure interaction applications where the small inclusions represent particles or fibers immersed into a fluid. To address these challenging applications using computational models, a thorough mathematical understanding of the fundamental mathematical aspects of the problem is extremely useful. As highlighted in [2], a mathematically-informed approach is a prerequisite for safe and reliable computations.

For these reasons, we focused on the fundamental aspects of the approximation of the problem. On one hand, we addressed the approximation of the Dirichlet boundary conditions on the inclusion by means of a simplified modeling approach based on the projection on Fourier modes. On the other hand we have studied the properties of the finite element method used for the approximation of the simplified model.

A particularly important question to be addressed is the robustness of this approach with respect to the size of the holes, which may become arbitrarily small with respect to the domain. To this purpose, we have studied three relevant problems: i) the full problem, corresponding to the standard enforcement of the non homogeneous Dirichlet boundary condition on the holes by means of Lagrange multipliers; ii) the simplified problem, characterized by the approximate weak enforcement of Dirichlet boundary conditions by projection on Fourier modes; iii) the limit problem obtained when the diameter of the holes vanishes. Understanding the mutual interaction of these problems characterizes what we call the model error in terms of the size of the holes. By means of suitable a-priori estimates of the model error and of the finite element approximation error, we provide guidelines to optimally balance the approximation parameters of the proposed simplified modeling approach. In particular a crucial issue is to provide a good approximation strategy of the Lagrange multiplier so that these theoretical results could be especially useful in view of forthcoming applications of this methodology to fluid-structure interaction problems.

## Acknowledgments

Fabien Lespagnol is supported by the project A new computational approach for the fluid-structure interaction of slender bodies immersed in three-dimensional flows granted by the Università Italo-Francese, in the framework Vinci 2019. Paolo Zunino is a member of the Gruppo Nazionale per il Calcolo Scientifico (GNCS), Istituto Nazionale di Alta Matematica (INdAM).

## References

[1] Xavier Antoine, Karim Ramdani, and Bertrand Thierry. Wide frequency band numerical approaches for multiple scattering problems by disks. *Journal of Algorithms & Computational Technology*, 6(2):241–259, 2012.

[2] Ivo Babuška, Ana Maria Soane, and Manil Suri. The computational modeling of problems on domains with small holes. *Computer Methods in Applied Mechanics and Engineering*, 322:563–589, 2017.

[3] S. Bertoluzza, A. Decoene, L. Lacouture, and S. Martin. Local error estimates of the finite element method for an elliptic problem with a dirac source term. *Numerical Methods for Partial Differential Equations*, 34(1):97–120, 2018.

[4] D. Boffi, N. Cavallini, and L. Gastaldi. The finite element immersed boundary method with distributed lagrange multiplier. *SIAM JOURNAL on Numerical Analysis*, 53(6):2584–2604, 2015.

[5] Daniele Boffi, Franco Brezzi, Michel Fortin, et al. *Mixed finite element methods and applications*, volume 44. Springer, 2013.
[6] Virginie Bonnaillie-Noël and Marc Dambrine. Interactions between moderately close circular inclusions: the dirichlet–laplace equation in the plane. *Asymptotic Analysis*, 84(3-4):197–227, 2013.

[7] Virginie Bonnaillie-Noël, Marc Dambrine, Sébastien Tordeux, and Grégory Vial. On moderately close inclusions for the laplace equation. *Comptes Rendus Mathématique*, 345(11):609–614, 2007.

[8] Franco Brezzi. On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers. *Publications mathématiques et informatique de Rennes, 1*(S4):1–26, 1974.

[9] Fabien Caubet, Carlos Conca, and Matías Godoy. On the detection of several obstacles in 2d stokes flow: topological sensitivity and combination with shape derivatives. *Inverse Problems & Imaging*, 10(2):327, 2016.

[10] Lucas Chesnel and Xavier Claeys. A numerical approach for the poisson equation in a planar domain with a small inclusion. *BIT Numerical Mathematics*, 56(4):1237–1256, 2016.

[11] Xavier Claeys and Francis Collino. Augmented galerkin schemes for the numerical solution of scattering by small obstacles. *Numerische Mathematik*, 116(2):243–268, 2010.

[12] Robert Dautray and Jacques-Louis Lions. Analyse mathématique et calcul numérique pour les sciences et les techniques. *Collection du Commissariat a l'Energie Atomique. Serie Scientifique*, 1985.

[13] Alexandre Ern and Jean-Luc Guermond. *Theory and practice of finite elements*, volume 159. Springer, 2004.

[14] Lawrence C. Evans. *Partial differential equations*, volume 19 of *Graduate Studies in Mathematics*. American Mathematical Society, Providence, RI, 1998.

[15] Jean Giroire. *Etude de quelques problemes aux limites extérieurs et résolution par équations intégrales*. PhD thesis, Paris 6, 1987.

[16] R. Glowinski and Yu. Kuznetsov. Distributed lagrange multipliers based on fictitious domain method for second order elliptic problems. *Computer Methods in Applied Mechanics and Engineering*, 196(8):1498–1506, 2007.

[17] R. Glowinski, T.-W. Pan, and J. Periaux. A fictitious domain method for external incompressible viscous flow modeled by navier-stokes equations. *Computer Methods in Applied Mechanics and Engineering*, 112(1-4):133–148, 1994.

[18] Wei Gong, Gengsheng Wang, and Ningning Yan. Approximations of elliptic optimal control problems with controls acting on a lower dimensional manifold. *SIAM JOURNAL on Control and Optimization*, 52(3):2008–2035, 2014.

[19] Ph Guillaume and K Sid Idris. The topological asymptotic expansion for the dirichlet problem. *SIAM JOURNAL on Control and Optimization*, 41(4):1042–1072, 2002.

[20] N. Hagmeyer, M. Mayr, I. Steinbrecher, and A. Popp. One-way coupled fluid–beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa. *Advanced Modeling and Simulation in Engineering Sciences*, 9(1), 2022.

[21] Anita Hansbo and Peter Hansbo. An unfitted finite element method, based on nitsche's method, for elliptic interface problems. *Computer methods in applied mechanics and engineering*, 191(47-48):5537–5552, 2002.

[22] Antoine Henrot and Michel Pierre. *Variation et optimisation de formes: une analyse géométrique*, volume 48. Springer Science & Business Media, 2006.

[23] George C Hsiao and Wolfgang L Wendland. *Boundary integral equations*. Springer, 2008.

[24] Ognyan Kounchev. *Multivariate polysplines: applications to numerical and wavelet analysis*. Academic Press, 2001.
[25] M. Kuchta, F. Laurino, K.-A. Mardal, and P. Zunino. Analysis and approximation of mixed-dimensional pdes on 3d-1d domains coupled with lagrange multipliers. *SIAM JOURNAL on Numerical Analysis*, 59(1):558–582, 2021.

[26] T. Köppl, E. Vidotto, B. Wohlmuth, and P. Zunino. Mathematical modeling, analysis and numerical approximation of second-order elliptic problems with inclusions. *Mathematical Models and Methods in Applied Sciences*, 28(5):953–978, 2018.

[27] Federica Laurino and Paolo Zunino. Derivation and analysis of coupled pdes on manifolds with high dimensionality gap arising from topological model reduction. *ESAIM: Mathematical Modelling and Numerical Analysis*, 53(6):2047–2080, 2019.

[28] B. Maury. A fat boundary method for the poisson problem in a domain with holes. *JOURNAL of Scientific Computing*, 16(3):319–339, 2001.

[29] B. Maury. Numerical analysis of a finite element/volume penalty method. *SIAM JOURNAL on Numerical Analysis*, 47(2):1126–1148, 2009.

[30] Vladimir Maz'Ya, Serguei Nazarov, and Boris Plamenevskij. *Asymptotic theory of elliptic boundary value problems in singularly perturbed domains*, volume 1. Springer Science & Business Media, 2000.

[31] Vladimir G Maz'ya and Sergei V Poborchi. *Differentiable Functions on Bad Domains*. WORLD SCIENTIFIC, 1998.

[32] Vladimir Gilelevich Maz'ya and Sergei Aleksandrovich Nazarov. Asymptotic behavior of energy integrals under small perturbations of the boundary near corner and conic points. *Trudy Moskovskogo Matematicheskogo Obshchestva*, 50:79–129, 1987.

[33] Vladimir Gilelevich Maz'ya, Sergei Aleksandrovich Nazarov, and Boris Alekseevich Plamenevskii. Asymptotic expansions of the eigenvalues of boundary value problems for the laplace operator in domains with small holes. *Izvestiya Rossiiiskoi Akademii Nauk. Seriya Matematicheskaya*, 48(2):347–371, 1984.

[34] Y. Mori, L. Ohm, and D. Spirn. Theoretical justification and error analysis for slender body theory. *Communications on Pure and Applied Mathematics*, 73(6):1245–1314, 2020.

[35] Y. Mori, L. Ohm, and D. Spirn. Theoretical justification and error analysis for slender body theory with free ends. *Archive for Rational Mechanics and Analysis*, 235(3):1905–1978, 2020.

[36] SA Nazarov. Asymptotic conditions at a point, selfadjoint extensions of operators, and the method of matched asymptotic expansions. *Translations of the American Mathematical Society-Series 2*, 193:77–126, 1999.

[37] Serguei A Nazarov and Jan Sokolowski. Asymptotic analysis of shape functionals. *Journal de Mathématiques pures et appliquées*, 82(2):125–196, 2003.

[38] L Ridgway Scott and Shangyou Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions. *Mathematics of Computation*, 54(190):483–493, 1990.

[39] Kahina Sid Idris. *Sensibilité topologique en optimisation de forme*. PhD thesis, Toulouse, INSA, 2001.

[40] I. Steinbrecher, A. Popp, and C. Meier. Consistent coupling of positions and rotations for embedding 1d cosserat beams into 3d solid volumes. *Computational Mechanics*, 69(3):701–732, 2022.

[41] R. Stenberg. On some techniques for approximating boundary conditions in the finite element method. *JOURNAL of Computational and Applied Mathematics*, 63(1-3):139–148, 1995.

[42] Antoni Zygmund. *Trigonometric series*, volume 1. Cambridge university press, 2002.
# A Appendix

In this appendix, we state several results that have been used in the previous sections. We recall that $\varepsilon \in (0, \varepsilon_{max})$ where $\varepsilon_{max}$ is the fixed parameter defined by $\varepsilon_{max} = \sup\{\varepsilon > 0 | \bar{\omega}_{\varepsilon} \subset \Omega\}$.

**Lemma A.1.** For $\Phi \in H^{\frac{1}{2}}(\partial\Omega)$ and $f \in L^2(\Omega)$ such that $\omega_{\varepsilon} \cap \text{supp} f = \emptyset$, the following problem

$$ \begin{cases} -\Delta v_{\varepsilon} = f & \text{in } \Omega_{\varepsilon}, \\ v_{\varepsilon} = \Phi & \text{on } \partial\Omega, \\ v_{\varepsilon} = 0 & \text{on } \partial\omega_{\varepsilon}, \end{cases} \qquad (A.1) $$

admits a unique weak solution in $H^1(\Omega_{\varepsilon})$. Moreover,

$$ \|v_{\varepsilon}\|_{1,\Omega_{\varepsilon}} \lesssim (\|\Phi\|_{\frac{1}{2},\partial\Omega} + \|f\|_{0,\Omega}). $$

*Proof.* The proof is similar to the proof of [[9], Lemma C.1] for the Stokes equations which is itself based on the results of [[39], Chapter 3]. In this proof, we will consider two cases, the first where $f = 0$ and the second where $\Phi = 0$ and conclude by linearity.

Let us first suppose that $f = 0$. We consider $v_{\varepsilon_{max}}$ the solution of (A.1) for $\varepsilon = \varepsilon_{max}$. It satisfies

$$ |v_{\varepsilon_{max}}|_{1,\Omega_{\varepsilon_{max}}} = \left( \int_{\Omega_{\varepsilon_{max}}} |\nabla v_{\varepsilon_{max}}|^2 dx \right)^{\frac{1}{2}} \lesssim \|\Phi\|_{\frac{1}{2},\partial\Omega}. \qquad (A.2) $$

Now consider $\tilde{v}_{\varepsilon_{max}} \in H^1(\Omega)$ the extension by 0 of $v_{\varepsilon_{max}}$ to all $\Omega$. Notice that since $\varepsilon < \varepsilon_{max}$, $\varepsilon\omega \subset \varepsilon_{max}\omega$ and $\Omega_{\varepsilon_{max}} \subset \Omega_{\varepsilon}$, so, by minimization of energy, we have

$$ |v_{\varepsilon}|_{1,\Omega_{\varepsilon}} \le |\tilde{v}_{\varepsilon_{max}}|_{1,\Omega_{\varepsilon}} = |v_{\varepsilon_{max}}|_{1,\Omega_{\varepsilon_{max}}}, $$

and thanks to equation (A.2),

$$ |v_{\varepsilon}|_{1,\Omega_{\varepsilon}} \lesssim \|\Phi\|_{\frac{1}{2},\partial\Omega}. \qquad (A.3) $$

Let $u_0$ be the solution of the problem (2.2). Since $\tilde{v}_{\varepsilon} - u_0 \in H_0^1(\Omega)$, we can apply the Poincaré inequality given by Lemma 3.2, we get

$$ \|v_{\varepsilon}\|_{0,\Omega_{\varepsilon}} = \|\tilde{v}_{\varepsilon}\|_{0,\Omega} \le \|\tilde{v}_{\varepsilon} - u_0\|_{0,\Omega} + \|u_0\|_{0,\Omega} \lesssim (|v_{\varepsilon}|_{1,\Omega_{\varepsilon}} + |u_0|_{1,\Omega}) + \|u_0\|_{0,\Omega}. $$

The well-posedness of the problem (2.2) also gives

$$ \|u_0\|_{1,\Omega} \lesssim \|\Phi\|_{\frac{1}{2},\partial\Omega}. \qquad (A.4) $$

Combining equations (A.3) and (A.4) finally gives us

$$ \|v_{\varepsilon}\|_{0,\Omega_{\varepsilon}} \lesssim \|\Phi\|_{\frac{1}{2},\partial\Omega}. $$

Now, let us suppose that $\Phi = 0$. Let us note that for all $v \in H_0^1(\Omega_{\varepsilon})$, denoting by $\tilde{v}$ the extension by zero of $v$ to $\Omega$, we have

$$ \|v\|_{0,\Omega_{\varepsilon}} = \|\tilde{v}\|_{0,\Omega} \lesssim \|\nabla \tilde{v}\|_{0,\Omega} = \|\nabla v\|_{0,\Omega_{\varepsilon}}. $$

Using this inequality and Lax-Milgram theorem, we get

$$ \|v_{\varepsilon}\|_{1,\Omega_{\varepsilon}} \lesssim \|f\|_{0,\Omega}. $$

This concludes the proof of the lemma.

□

**Lemma A.2.** For any $L \in \mathbb{R}$, the problem

$$ \begin{cases} -\Delta v_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\ v_{\varepsilon} = 0 & \text{on } \partial\Omega, \\ v_{\varepsilon} = L & \text{on } \partial\omega_{\varepsilon}, \end{cases} $$

admits a unique weak solution in $H^1(\Omega_{\varepsilon})$. Moreover,

$$ \|v_{\varepsilon}\|_{1,\Omega_{\varepsilon}} \lesssim (-\ln(\varepsilon))^{-\frac{1}{2}}|L|. $$
Proof. As for the previous lemma, the proof is adapted from [[9], Lemma C.3] for the Stokes equations itself presented in [[39], Chapter 3]. Let us define $w_\epsilon$ the unique solution of the system

$$ \begin{cases} -\Delta w_\epsilon = 0 & \text{in } \omega_{\epsilon \max / \epsilon} \setminus \bar{\omega}, \\ w_\epsilon = 0 & \text{on } \partial \omega_{\epsilon \max / \epsilon}, \\ w_\epsilon = L & \text{on } \partial \omega. \end{cases} $$

We also consider the function $\hat{v}_\epsilon$ defined on $\Omega_\epsilon \setminus \bar{\omega}$ by $\hat{v}_\epsilon(\mathbf{x}) = v_\epsilon(\epsilon\mathbf{x})$ for all $\mathbf{x} \in \Omega_\epsilon \setminus \bar{\omega}$. The function $\hat{v}_\epsilon$ satisfies

$$ \left\{ \begin{array}{ll} -\Delta \hat{v}_{\epsilon} = 0 & \text{in } \Omega_{\epsilon} \setminus \bar{\omega}, \\ \hat{v}_{\epsilon} = 0 & \text{on } \frac{1}{\epsilon} \partial \Omega, \\ \hat{v}_{\epsilon} = L & \text{on } \partial \omega. \end{array} \right. \qquad (A.5) $$

Notice that we have $\bar{\omega} \subset \omega_{\epsilon \max / \epsilon} \subset \Omega_\epsilon$. Now we consider $\tilde{w}_\epsilon$ the extension of $w_\epsilon$ to $\Omega_\epsilon \setminus \bar{\omega}$ by zero in the outer part of the extended domain. Therefore, by the principle of minimization of energy, we have

$$ |v_{\epsilon}|_{1, \Omega_{\epsilon}} = |\hat{v}_{\epsilon}|_{1, \Omega_{\epsilon} \setminus \bar{\omega}} \le |\tilde{w}_{\epsilon}|_{1, \Omega_{\epsilon} \setminus \bar{\omega}} = |w_{\epsilon}|_{1, \omega_{\epsilon \max / \epsilon} \setminus \bar{\omega}}. \qquad (A.6) $$

A computation provides for all $\mathbf{x} \in \omega_{\epsilon \max / \epsilon} \setminus \bar{\omega}$,

$$ w_{\epsilon}(\mathbf{x}) = L \frac{\ln(\epsilon_{\max} / \epsilon) - \ln(|\mathbf{x}|)}{\ln(\epsilon_{\max} / \epsilon)} $$

and

$$ |w_{\epsilon}|_{1, \omega_{\epsilon \max / \epsilon} \setminus \bar{\omega}} \lesssim (-\ln(\epsilon))^{-\frac{1}{2}} |L|. $$

So we get

$$ |v_{\epsilon}|_{1, \Omega_{\epsilon}} \lesssim (-\ln(\epsilon))^{-\frac{1}{2}} |L|. \qquad (A.7) $$

Finally, we consider $\tilde{v}_\epsilon$ the extension of $v_\epsilon$ to $\Omega$ by $L$. Since this extension is in $H_0^1(\Omega)$, we can use the Poincaré inequality given by Lemma 3.2,

$$ \|v_{\epsilon}\|_{0, \Omega_{\epsilon}} \le \|\tilde{v}_{\epsilon}\|_{0, \Omega} \lesssim |\tilde{v}_{\epsilon}|_{1, \Omega} \lesssim |v_{\epsilon}|_{1, \Omega_{\epsilon}}. $$

Using equation (A.7), we get the result.

<figure><img src="image_18.png" /></figure>

**Lemma A.3.** For any $\varphi \in H^{\frac{1}{2}}(\partial \omega_\epsilon)$, the problem

$$ \begin{cases} -\Delta v_{\epsilon} = 0 & \text{in } \Omega_{\epsilon}, \\ v_{\epsilon} = 0 & \text{on } \partial \Omega, \\ v_{\epsilon} = \varphi & \text{on } \partial \omega_{\epsilon}, \end{cases} $$

admits a unique weak solution in $H^1(\Omega_\epsilon)$. Moreover,

$$ \|v_{\epsilon}\|_{1, \Omega_{\epsilon}} \lesssim \|\varphi(\epsilon\mathbf{x})\|_{\frac{1}{2}, \partial\omega}. $$

Proof. Once again, this proof is adapted from a similar proof conducted in [[9], Lemma B.2 and Lemma 4.2] for the Stokes problem which is itself inspired from a proof described in [[39], Chapter 3]. Lax-Milgram theorem allows to prove that the problem

$$ \left\{ \begin{array}{ll} -\Delta V = 0 & \text{in } \mathbb{R}^2 \setminus \bar{\omega}, \\ V = \varphi(\epsilon \mathbf{x}) & \text{in } \partial \omega, \end{array} \right. \qquad (A.8) $$

is well posed and has a unique solution in

$$ W_0^{1,2}(\mathbb{R}^2 \setminus \bar{\omega}) = \{u \in \mathcal{D}'(\mathbb{R}^2 \setminus \bar{\omega}) \mid \ln(\rho)^{-1} u \in L_{-1}^2(\mathbb{R}^2 \setminus \bar{\omega}), \nabla u \in L^2(\mathbb{R}^2 \setminus \bar{\omega})\} $$

where

$$ \rho(\mathbf{x}) = (1 + |\mathbf{x}|^2)^{\frac{1}{2}} $$
and

$$L_{-1}^{2}(\mathbb{R}^{2}\setminus\bar{\omega}) = \{u \in D'(\mathbb{R}^{2}\setminus\bar{\omega}) \mid \rho^{-1}u \in L^{2}(\mathbb{R}^{2}\setminus\bar{\omega})\},$$

(see [15] for example). We will try to give an explicit representation of $V$. By setting $-\Delta V = 0$ in $\omega$, the problem (A.8) has a unique solution in $\mathbb{R}^2$ and we have that

$$-\Delta V = \nabla V \cdot \mathbf{n} \delta_{\partial\omega}$$

in $D'(\mathbb{R}^2)$ where $\mathbf{n}$ is the exterior normal on $\partial\omega$. Now let us define

$$W = E * (\nabla V \cdot \mathbf{n} \delta_{\partial\omega})$$

where $E$ is the fundamental solution of the Laplace equation given for $\mathbf{x} \in \mathbb{R}^2 \setminus \{0\}$ by

$$E(\mathbf{x}) = -\frac{\ln(|\mathbf{x}|)}{2\pi},$$

and $*$ is the convolution product. We have

$$-\Delta W = \nabla V \cdot \mathbf{n} \delta_{\partial\omega}$$

in $D'(\mathbb{R}^2)$. Then $V-W$ is a harmonic tempered distribution. A classical result of Fourier analysis states that harmonic tempered distribution are polynomials (see [12] for example). Then $V = L+W$ with $L$ a polynomial and for $\mathbf{y} \in \mathbb{R}^2\setminus\omega$,

$$W(\mathbf{y}) = \int_{\partial\omega} t(\mathbf{x})E(\mathbf{y} - \mathbf{x})ds(\mathbf{x}),$$

with $t(\mathbf{x}) = \nabla V \cdot \mathbf{n}$. Using a Taylor development for $E$, we get

$$E(\mathbf{y} - \mathbf{x}) = E(\mathbf{y}) - \nabla E(\mathbf{y} - \alpha\mathbf{x}) \cdot \mathbf{x}$$

for some $\alpha \in (0, 1)$. We then have

$$W(\mathbf{y}) = E(\mathbf{y}) \int_{\partial\omega} t(\mathbf{x})ds(\mathbf{x}) - \int_{\partial\omega} t(\mathbf{x})\nabla E(\mathbf{y} - \alpha\mathbf{x}) \cdot \mathbf{x}ds(\mathbf{x}).$$

Let us denote

$$U(\mathbf{y}) = \int_{\partial\omega} t(\mathbf{x})\nabla E(\mathbf{y} - \alpha\mathbf{x}) \cdot \mathbf{x}ds(\mathbf{x}).$$

By computation, we get that $U(\mathbf{y}) = o(1/|\mathbf{y}|)$ when $|\mathbf{y}| \to \infty$ so $(\ln(\rho))^{-1}U \in L_{-1}^2(\mathbb{R}^2\setminus\bar{\omega})$. As $(\ln(\rho))^{-1}V \in L_{-1}^2(\mathbb{R}^2\setminus\bar{\omega})$ and $(\ln(\rho))^{-1} \notin L_{-1}^2(\mathbb{R}^2\setminus\bar{\omega})$, we necessarily have that

$$\int_{\partial\omega} t(\mathbf{x})ds(\mathbf{x}) = 0$$

and that $L$ is a constant. By computation, we have that for $\|\mathbf{y}\|$ sufficiently large,

$$|W(\mathbf{y})| \lesssim \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2}, \partial\omega} \frac{1}{\|\mathbf{y}\|} \quad \text{and} \quad |\nabla W(\mathbf{y})| \lesssim \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2}, \partial\omega} \frac{1}{\|\mathbf{y}\|^2}.$$

Let $A > 0$ such that the previous inequality is satisfied for $\|\mathbf{y}\| > A$. We have for $\|\mathbf{y}\| > A$,

$$|L| \lesssim |V(\mathbf{y})| + \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2}, \partial\omega} \frac{1}{\|\mathbf{y}\|}.$$

Integrating for $\|\mathbf{y}\| > A$, we get

$$|L| \left( \int_{\|\mathbf{y}\| > A} \frac{1}{\log(\|\mathbf{y}\|)^2 \|\mathbf{y}\|^2} \right)^{\frac{1}{2}} \lesssim \left( \int_{\|\mathbf{y}\| > A} \frac{|V(\mathbf{y})|^2}{\log(\|\mathbf{y}\|)^2 \|\mathbf{y}\|^2} \right)^{\frac{1}{2}} + \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2}, \partial\omega} \left( \int_{\|\mathbf{y}\| > A} \frac{1}{\log(\|\mathbf{y}\|)^2 \|\mathbf{y}\|^4} \right)^{\frac{1}{2}}.$$
The fact that A is independent of $\varphi$ and the well-posedness of the problem (A.8) give

$$|L| \lesssim \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega}.$$

Using similar computations as in [[19], Lemma 7.1], we also have

$$\|W(\frac{\mathbf{x}}{\varepsilon})\|_{1,\Omega_{\varepsilon}} \lesssim \varepsilon \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega}.$$

We then define $z_{\varepsilon} := v_{\varepsilon} - W(\frac{\mathbf{x}}{\varepsilon})$, $z_{\varepsilon}$ satisfies

$$\begin{cases} -\Delta z_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\ z_{\varepsilon} = -W(\frac{\mathbf{x}}{\varepsilon}) & \text{on } \partial\Omega, \\ z_{\varepsilon} = L & \text{on } \partial\omega_{\varepsilon}. \end{cases}$$

Using Lemma A.1 and Lemma A.2 we get that

$$\begin{aligned} \|z_{\varepsilon}\|_{1,\Omega_{\varepsilon}} & \lesssim \|W(\frac{\mathbf{x}}{\varepsilon})\|_{\frac{1}{2},\partial\Omega} + (-\ln(\varepsilon))^{-\frac{1}{2}}|L|, \\ & \lesssim \varepsilon \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega} + (-\ln(\varepsilon))^{-\frac{1}{2}}\|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega}. \end{aligned}$$

So finally, we get

$$\begin{aligned} \|v_{\varepsilon}\|_{1,\Omega_{\varepsilon}} & \lesssim \|z_{\varepsilon}\|_{1,\Omega_{\varepsilon}} + \|W(\frac{\mathbf{x}}{\varepsilon})\|_{1,\Omega_{\varepsilon}}, \\ & \lesssim \varepsilon \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega} + (-\ln(\varepsilon))^{-\frac{1}{2}} \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega} + \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega}, \\ \|v_{\varepsilon}\|_{1,\Omega_{\varepsilon}} & \lesssim \|\varphi(\varepsilon\mathbf{x})\|_{\frac{1}{2},\partial\omega}. \end{aligned}$$

□