# An a Posteriori Error Estimate for a 0D-2D Coupled Model

Hussein Albazzal, Alexei Lozinski, Roberta Tittarelli

► To cite this version:

Hussein Albazzal, Alexei Lozinski, Roberta Tittarelli. An a Posteriori Error Estimate for a 0D-2D Coupled Model. SIAM Journal on Scientific Computing, 2024, 46 (6), pp.B903-B924. 10.1137/22M161651X. http://05221594

HAL Id: hal-05221594

https://hal.science/hal-05221594v1

Submitted on 25 Aug 2025

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# An a posteriori error estimate for a 0D/2D coupled model *

Hussein Albazzal†, Alexei Lozinski‡, and Roberta Tittarelli§

**Abstract.** This work is motivated by the need of efficient numerical simulations of gas flows in the serpentine channels used in proton-exchange membrane fuel cells. In particular, we consider the Poisson problem in a 2D domain composed of several long straight rectangular sections and of several bends corners. In order to speed up the resolution, we propose a 0D model in the rectangular parts of the channel and a Finite Element resolution in the bends. To find a good compromise between precision and time consuming, the challenge is double: how to choose a suitable position of the interface between the 0D and the 2D models and how to control the discretization error in the bends. We shall present an *a posteriori* error estimator based on an equilibrated flux reconstruction in the subdomains where the Finite Element method is applied. The estimates give a global upper bound on the error measured in the energy norm of the difference between the exact and approximate solutions on the whole domain. They are guaranteed, meaning that they feature no undetermined constants. (global) Lower bounds for the error are also derived. An adaptive algorithm is proposed to use smartly the estimator for aforementioned double challenge. A numerical validation of the estimator and the algorithm completes the work.

**Key words.** *A posteriori* error estimate, mixed dimensional coupling, adaptive algorithm, Finite Element Method, FreeFEM.

**MSC codes.** 65N15; 65N30; 65N50

**1. Introduction.** The present work is motivated by models of the serpentine cathode-anode flow channels in Proton Exchange Membrane Fuel Cells (PEMFC), cf. [9]. One of the computationally intensive tasks in this modeling is to solve the steady-state incompressible Navier-Stokes equations describing the gas flow in the long channels, characterised by very stretched rectangular regions linked by relatively small bends, cf. Figure 1. Typically, one prescribes the Poiseuille flow as boundary conditions on the inlet and outlet boundaries of the channel, and no-slip boundary condition on the wall, cf. for example [15]. In order to speed up the computations, the idea is to develop a coupled model as follows (here, in the 2D setting): in the rectangular regions of the domain, the flow is approximated by simple analytical solutions, namely the Poiseuille flow which is very accurate sufficiently far from the bends (we call this 0D model), while keeping the original governing equations (the 2D model) in the bend regions. We refer to this as the 0D/2D model. There exists different ways to derive a coupled model. In [12, 5, 4], 0D/3D coupling is obtained for the time dependent Navier-Stokes system by integrating the governing equations on a section and substituting an appropriate closure approximation. In [6] and in [10] an asymptotic analysis is used to get respectively the 1D/3D and 1D/2D coupled simplified models for time dependent Navier-

<footnote>**Funding:** This work was supported by the French 'Investissements d'Avenir' program, project Agence Nationale de la Recherche (ISITE-BFC) (contract ANR-15-IDEX-0003).</footnote>

<footnote>†Université de Franche-Comté, CNRS, LmB, F-25000 Besançon, France (hussein.albazzal@univ-fcomte.fr).</footnote>

<footnote>‡Université de Franche-Comté, CNRS, LmB, F-25000 Besançon, France (alexei.lozinski@univ-fcomte.fr).</footnote>

<footnote>§ SUPMICROTECH, CNRS, institut FEMTO-ST, F-25000 Besançon, France (roberta.tittarelli@ens2m.fr).</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Gaz channels in a flow plate of a fuel cell [16].</figcaption></figure>

Stokes equations.

In this paper, we consider a simple toy model, replacing the steady Navier-Stokes system by the Poisson equation, to test the basic ideas of the 0D/2D model coupling. A 0D/2D coupled model for Poisson equation is derived in [13] and [11] by the asymptotic analysis. In the stretched rectangular portion of the domain $\Omega$, which we call in this paper $\Omega'$ and refer to as the 0D domain, cf. Fig. 2, the solution is approximated by a simple explicit expression (the analogue of the Poiseuille velocity). In the remaining part of the domain, denoted by $\tilde{\Omega} = \Omega \setminus \overline{\Omega'}$ and referred to as the 2D domain, we approximate the solution of the original Poisson equation (the 2D model) using the Finite Element (FE) method. The 0D and 2D domains are separated by an interface $\Gamma$. We use a simple coupling condition on the interface $\Gamma$: we impose the continuity of the coupled solution on $\Gamma$. The error between the “original 2D model” and the “0D/2D coupled discretized model” comes thus from two sources: the FE discretization in $\tilde{\Omega}$ and the coupling error provoked by the simple approximation in $\Omega'$ and controlled by the position of the interface $\Gamma$.

In the present paper, we propose an *a posteriori* error estimator for the 0D/2D coupled model which can be interpreted as a sum of two contributions: the first one measuring the error due to the simplification introduced by the 0D domain $\Omega'$ and the second one measuring the discretization error in the 2D domain $\tilde{\Omega}$. Equilibrating the two contributions, under a given tolerance, enables us to determine the position of the interface $\Gamma$, reducing the size of the 2D domain $\tilde{\Omega}$ to be discretized, and to construct an optimized mesh on $\tilde{\Omega}$. Our error estimator is based on the flux reconstruction technique as in [3, 2]. The originality of our work consists in a new flux reconstruction $\sigma_h$, which is defined on the whole domain $\Omega = \Omega' \cup \Gamma \cup \tilde{\Omega}$, i.e. both on the discretized (2D) and the non-discretized (0D) regions.

The article is organized as follows. In the next section, we introduce the governing equations, the geometry of the domain, and advocate for coupled simplified model. Our *a posteriori*
error is presented in section 3. In section 4, we prove the global upper and lower bounds for the error with respect to the estimator, called respectively the global reliability and efficiency of the estimator. The reliability is guaranteed, i.e. the upper bound does not contain any unknown constants. Its proof is quite straightforward. The main technical hurdle is in the proof of the global efficiency of the estimator. Finally, we propose in section 5 an algorithm that uses our estimator both to choose the interface position and to make the local mesh refinement. The section is concluded by numerical tests. In section 6 we provide conclusions and forthcoming works.

**2. Original problem and approximated coupled problem.** Let $\Omega \subset \mathbb{R}^2$ be the polygonal bounded domain, as shown in Figure 2, resembling a portion of the channel in a real fuel cell from Figure 1. Domain $\Omega$ is split into two parts $\Omega'$ and $\tilde{\Omega}$, separated by the interface $\Gamma$, which is placed at coordinate $x_{\Gamma} \in (0, L)$, with $L > 0$ representing the length of the straight part of the channel. Thus, $\Omega' = (0, x_{\Gamma}) \times (0, R)$, with $R > 0$ standing for the width of the channel. The boundary $\partial\Omega$ of $\Omega$ is partitioned in the inlet $\Gamma_{\text{in}}$ and outlet $\Gamma_{\text{out}}$ parts, where we impose a Poiseuille-like profile, and in the remaining part $\Gamma_{\text{wall}}$, where we impose the homogeneous Dirichlet conditions. We consider the problem with the scalar unknown $u$ on $\Omega$ such that

$$ (2.1a) \qquad - \Delta u = f \quad \text{in } \Omega, $$

$$ (2.1b) \qquad u = g \quad \text{on } \partial\Omega, $$

where $f = \frac{12u_{\text{av}}}{R^2}$, and

$$ (2.2) \qquad g = \begin{cases} S(y) & \text{on } \Gamma_{\text{in}}, \\ 0 & \text{on } \Gamma_{\text{wall}}, \\ S(y + W + R) & \text{on } \Gamma_{\text{out}}, \end{cases} $$

where $W > 0$ stands for the straight part of the bend, cf. Figure 2, and

$$ S(y) = \frac{6u_{\text{av}}}{R^2} y(R - y), $$

with $u_{\text{av}}$ the average value of the solution on the inlet/outlet. The weak formulation is then given by: find $u \in H_g^1(\Omega) := \{v \in H^1(\Omega) : u = g$ on $\partial\Omega$} $ such that

$$ (2.3) \qquad \int_{\Omega} \nabla u \cdot \nabla v = \int_{\Omega} fv \quad \text{for all } v \in H_0^1(\Omega). $$

A typical solution to this problem is shown at Figure 3. We observe that it is indeed of simple form, independent of $x$, in the rectangular portion of the channel, to the left of the bend. It is where we place the subdomain $\Omega'$. Let us now derive a simplified model in $\Omega'$ using techniques as in [11, 13]. Let us choose as characteristic constants $R, L, u_{\text{av}}$ for $x, y, u$, and consider the dimensionless quantities: $\bar{u} = \frac{u}{u_{\text{av}}}$, $\bar{x} = \frac{x}{L}$, $\bar{y} = \frac{y}{R}$. Assuming $R \ll L$, let us also define $\epsilon = \frac{R}{L}$. Then, rewriting (2.1a) on $\Omega'$ in dimensionless variables, we get

$$ (2.4) \qquad - \epsilon^2 \partial_{\bar{x}\bar{x}}^2 \bar{u} - \partial_{\bar{y}\bar{y}}^2 \bar{u} = \bar{f} \quad \text{for } \bar{y} \in (0, 1), $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Partition of the domain $\Omega$ in domains for the coupled 0D/2D model. $x_\Gamma$ is the $x$-coordinate of the interface position, $L, R, W, L_{out} > 0$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Solution to problem (2.1) approximated by $\mathbb{P}_2$-FE on a mesh of size approximately $h = 0.005$, used as reference solution for numerical tests in subsection 5.1.</figcaption></figure>

with $\bar{f} = \frac{R^2 f}{u_{av}} = 12$. Neglecting the term multiplied by $\epsilon^2 \ll 1$ and recovering the variables with dimensions, we get the approximate equation on $\Omega'$ whose solution is denoted by $u'$:

$$ (2.5) \qquad - \partial_{yy}^2 u' = 12 \frac{u_{av}}{R^2} \quad \text{in } \Omega'. $$

Combining this with the boundary conditions on $\Gamma_{wall}$, i.e. $u' = 0$ for $y = 0$ and $y = R$, leads to the approximation

$$ (2.6) \qquad u' = S(y) \quad \text{in } \Omega'. $$

This approximation is evidently compatible with the inlet boundary condition on $\Gamma_{in}$.

In $\tilde{\Omega}$, we shall solve numerically (2.1) by the FE method. Let $\mathcal{T}_h$ be a triangular mesh on $\tilde{\Omega}$ and $V_h$ the FE space of piecewise $\mathbb{P}_k$ polynomials with $k \ge 2$, i.e. $V_h = \{v_h \in H^1(\tilde{\Omega}) : v_{h|_K} \in \mathbb{P}_k \quad \forall K \in \mathcal{T}_h\}$. The simplest coupling condition between the approximation $u'$ on $\Omega'$ and the numerical approximation $\tilde{u}_h$ on $\tilde{\Omega}$ is $\tilde{u}_h = u'$ on $\Gamma$. We thus search $\tilde{u}_h \in V_h$ such that $\tilde{u}_h = g$ on $\partial\tilde{\Omega}$, and

$$ (2.7) \qquad \int_{\tilde{\Omega}} \nabla \tilde{u}_h \cdot \nabla v_h = \int_{\tilde{\Omega}} f v_h, \quad \forall v_h \in V_h, v_h = 0 \quad \text{on } \partial \tilde{\Omega}. $$
Finally, the coupled approximated solution $u_h^c$ is defined as

$$ (2.8) \qquad u_h^c = \begin{cases} u' \text{ in } \Omega', \\ \tilde{u}_h \text{ in } \tilde{\Omega}. \end{cases} $$

**3. A posteriori error estimator.** In this section, we introduce an *a posteriori* error estimator with guaranteed reliability and provable efficiency for the error defined by

$$ (3.1) \qquad e = \| \nabla u - \nabla u_h^c \|_{\Omega} $$

in the spirit of [3]. To this end, we should define a reconstructed flux $\sigma_h^c \in H(\operatorname{div}, \Omega)$ such that $\nabla \cdot \sigma_h^c = f$ in $\Omega$ and such that its construction be local. This is done below, cf. (3.7), as the sum of local contributions from the patches around the mesh nodes and a patch attached to the interface $\Gamma$. Then, the error estimator will be defined as

$$ (3.2) \qquad \eta = \| \sigma_h^c + \nabla u_h^c \| . $$

Let $\mathcal{V}_h$ and $\mathcal{V}_h^\Gamma$ be the sets respectively of all the nodes of the mesh $\mathcal{T}_h$ of $\tilde{\Omega}$ and of all the nodes lying on the interface $\Gamma$. Let $\psi^\Gamma$ be the continuous function on $\Omega$ defined on rectangular portion $[0, x_\Gamma + R] \times [0, R]$ of $\Omega$, supposing that $x_\Gamma + R \le L$, by

$$ \psi^\Gamma(x, y) = \begin{cases} 1, & \text{for } x < x_\Gamma, 0 \le y \le R, \\ \frac{x_\Gamma + R - x}{R}, & \text{for } x \in [x_\Gamma, x_\Gamma + R], 0 \le y \le R, \\ 0, & \text{for } x > x_\Gamma + R, 0 \le y \le R, \end{cases} $$

and extended by 0 everywhere else. We also introduce a piecewise affine on mesh $\mathcal{T}_h$ version of $\psi^\Gamma$, defining

$$ \psi_h^\Gamma(x, y) = \begin{cases} \sum_{a \in \mathcal{V}_h} \psi^\Gamma(a) \psi^a(x, y) & \text{if } (x, y) \in \tilde{\Omega} \cup \Gamma, \\ 1 & \text{if } (x, y) \in \Omega', \end{cases} $$

where $\psi^a$ is the hat function associated to node $a \in \mathcal{V}_h$, i.e. the piecewise $\mathbb{P}_1$ function taking the value 1 at a and vanishing at all the other mesh nodes. It follows that we have the partition of unity on $\Omega$:

$$ (3.3) \qquad 1 = \psi_h^\Gamma + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} (1 - \psi^\Gamma(a)) \psi^a. $$

**Remark 3.1.** If $x_\Gamma + R > L$, then the function $\psi^\Gamma$ defined above becomes discontinuous and this prevents us from proving the optimality estimate (4.2) in the manner of Theorem 4.1. However, $\psi_h^\Gamma$ remains continuous even if $x_\Gamma + R > L$. Partition of unity also remains valid, so that the construction of the flux can be carried out as explained below. Hence, this hypothesis $x_\Gamma + R \le L$ is not necessary for the implementation of our error estimator. The numerical results show that it remains robust even if this theoretical assumption is violated.
Let us define the patches $\omega^a = \text{supp } \psi^a$ for all the nodes $a \in \mathcal{V}_h$ and the patch attached to $\Gamma$: $\omega^\Gamma = \text{supp}(\psi_h^\Gamma) \cap \tilde{\Omega}$. Consider the three following cases to define the $\mathbf{H}(\text{div})$-conforming Raviart-Thomas piecewise polynomials spaces on patches. We remind that $\mathbf{n}$ denotes the outward unit normal vector of patches. For the case where $a$ is an internal node of $\tilde{\Omega}$, that is $a \in \mathcal{V}_h \cap \tilde{\Omega}$, we define

$$ 
\begin{aligned}
\Sigma_h^a &= \{\tau_h \in \mathbf{H}(\text{div}, \omega^a); \tau_{h|T} \in \mathbf{RT}_k(T) \ \forall T \in \omega^a, \tau_h \cdot \mathbf{n} = 0 \text{ on } \partial \omega^a \}, \\
Q_h^a &= \{q_h \in L^2(\omega^a); q_{h|T} \in \mathbb{P}_k(T) \ \forall T \in \omega^a, \int_{\omega^a} q_h = 0 \},
\end{aligned}
$$

where $\mathbf{RT}_k(T) = [\mathbb{P}_k(T)]^2 + \mathbf{x}\mathbb{P}_k(T)$ with $T$ a mesh element. For the case where $a$ is a node on the boundary of $\tilde{\Omega}$ excluded the interface $\Gamma$, that is $a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma \cap \partial \tilde{\Omega}$, we define

$$ 
\begin{aligned}
\Sigma_h^a &= \{\tau_h \in \mathbf{H}(\text{div}, \omega^a); \tau_{h|T} \in \mathbf{RT}_k(T) \ \forall T \in \omega^a, \tau_h \cdot \mathbf{n} = 0 \text{ on } \partial \omega^a \cap \tilde{\Omega}\}, \\
Q_h^a &= \{q_h \in L^2(\omega^a); q_{h|T} \in \mathbb{P}_k(T) \ \forall T \in \omega^a\}.
\end{aligned}
$$

Finally, on the patch $\omega^\Gamma$ attached to the interface $\Gamma$, we define

$$ 
\begin{aligned}
\Sigma_h^\Gamma &= \{\tau_h \in \mathbf{H}(\text{div}, \omega^\Gamma); \tau_{h|T} \in \mathbf{RT}_k(T) \ \forall T \in \omega^\Gamma, \tau_h \cdot \mathbf{n} = 0 \text{ on } \Gamma \cup \tilde{\Gamma}\}, \\
Q_h^\Gamma &= \{q_h \in L^2(\omega^\Gamma); q_{h|T} \in \mathbb{P}_k(T) \ \forall T \in \omega^\Gamma\},
\end{aligned}
$$

with $\tilde{\Gamma} = \partial \omega^\Gamma \cap \tilde{\Omega}$. These are the appropriate spaces to define reconstructed fluxes, cf. the following Definition 3.2. This definition is completely standard, cf. [3], in the case of patches centered at the mesh nodes, but we extend it also to the patch attached to the interface $\Gamma$. This is a novelty of the present work.

**Definition 3.2.** For any mesh node $a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma$ or for $a = \Gamma$, recalling the corresponding spaces $\Sigma_h^a$ and $Q_h^a$ defined above, let $\sigma_h^a \in \Sigma_h^a$ and $p_h^a \in Q_h^a$ be the solution to

$$ (3.4a) \qquad \int_{\omega^a} \sigma_h^a \cdot \tau_h + \int_{\omega^a} p_h^a \operatorname{div} \tau_h = - \int_{\omega^a} (\nabla \tilde{u}_h) \psi^a \cdot \tau_h \qquad \forall \tau_h \in \Sigma_h^a, $$

$$ (3.4b) \qquad \int_{\omega^a} q_h \operatorname{div} \sigma_h^a = \int_{\omega^a} (f \psi^a - \nabla \tilde{u}_h \cdot \nabla \psi^a) q_h \qquad \forall q_h \in Q_h^a. $$

Each $\sigma_h^a$ is extended by 0 on $\tilde{\Omega}$ outside $\omega^a$.

**Remark 3.3.** Relation (3.4b) is actually satisfied for all the test functions $q_h$ in the discontinuous $\mathbb{P}_k$ space, without the constraint of vanishing average, even at the internal nodes $a \in \mathcal{V}_h \cap \tilde{\Omega}$. Indeed, taking any internal node $a$, we have thanks to the divergence theorem and $\sigma_h^a \in \Sigma_h^a$,

$$ (3.5) \qquad \int_{\omega^a} \operatorname{div} \sigma_h^a = - \int_{\partial \omega^a} \sigma_h^a \cdot \mathbf{n} = 0. $$

Thanks to (2.7) with $v_h = \psi^a$, we also have

$$ (3.6) \qquad \int_{\omega^a} (f \psi^a - \nabla \tilde{u}_h \cdot \nabla \psi^a) 1 = 0. $$

Hence, (3.4b) is satisfied for $q_h = 1$ and thus for any $q_h \in Q_h^a + \operatorname{span}(1)$.
We now define the flux $\sigma_h^c$ on $\Omega$ as

$$ (3.7) \qquad \sigma_h^c = \begin{cases} \sigma' & \text{in } \Omega', \\ \tilde{\sigma}_h & \text{in } \tilde{\Omega}, \end{cases} $$

where

$$ (3.8) \qquad \sigma' = -\nabla u' \text{ in } \Omega', $$

and

$$ (3.9) \qquad \tilde{\sigma}_h = \sigma_h^\Gamma + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} (1 - \psi^\Gamma(a)) \sigma_h^a \text{ in } \tilde{\Omega}, $$

with $\sigma_h^\Gamma$ and $\sigma_h^a$ from Definition 3.2.

**Proposition 3.4.** The flux $\sigma_h^c$ defined in (3.7) is such that $\sigma_h^c \in \mathbf{H}(\text{div}, \Omega)$ and $\text{div}\sigma_h^c = f$ in $\Omega$, where $f$ is the source term of (2.1).

*Proof.* On $\Omega'$, we obviously have

$$ (3.10) \qquad \text{div}\sigma' = -\Delta u' = f \text{ in } \Omega'. $$

In particular, $\sigma' \in \mathbf{H}(\text{div}, \Omega')$.

On $\tilde{\Omega}$, we have $\sigma_h^a \in \mathbf{H}(\text{div}, \tilde{\Omega})$ for any $a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma$ or $a = \Gamma$ thanks to the boundary conditions in the definition of $\Sigma_h^a$ and to the extension by 0 outside $\omega^a$. Moreover, thanks to Remark 3.3, equation (3.4b) can be rewritten as

$$ \int_{\omega^a} q_h \text{div} \sigma_h^a = \int_{\omega^a} (f \psi^a - \nabla \tilde{u}_h \cdot \nabla \psi^a) q_h \quad \forall q_h \in \mathcal{Q}_h, $$

where $\mathcal{Q}_h = \mathbb{P}_k(\mathcal{T}_h)$ is the broken polynomial space of degree $k \ge 2$ on the mesh $\mathcal{T}_h$. Thus, by (3.9),

$$ \begin{align*} \int_{\tilde{\Omega}} (\text{div}\tilde{\sigma}_h) q_h &= \int_{\tilde{\Omega}} (\text{div}\sigma_h^\Gamma) q_h + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} (1 - \psi^\Gamma(a)) \int_{\tilde{\Omega}} (\text{div}\sigma_h^a) q_h \\ &= \int_{\tilde{\Omega}} \left( f(\psi_h^\Gamma + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} (1 - \psi^\Gamma(a))\psi^a) - \nabla \tilde{u}_h \cdot \nabla (\psi_h^\Gamma + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} (1 - \psi^\Gamma(a))\psi^a) \right) q_h \\ &= \int_{\tilde{\Omega}} f q_h, \end{align*} $$

where to pass from the second to the third line, we have used the partition of unity (3.3). Since $\text{div}\mathbf{R}\mathbf{T}_k(T) = \mathbb{P}_k(T)$ on any mesh element $T$, this implies

$$ (3.11) \qquad \text{div}\tilde{\sigma}_h = f \text{ in } \tilde{\Omega}. $$

In view of (3.10) and (3.11) and the definition (3.7) of $\sigma_h^c$, it remains to prove that the normal component of $\sigma_h^c$ is continuous across $\Gamma$. This is evident by construction since $\sigma' \cdot \mathbf{n} = \tilde{\sigma}_h \cdot \mathbf{n} = 0$ on $\Gamma$. ■
4. Upper and lower bounds.

Theorem 4.1 (Global reliability and efficiency). Let $u$ be the weak solution of (2.3), $u_h^c$ be the solution (2.8) of the coupled 0D/2D problem, $e$ and $\eta$ be the error and the error estimator defined respectively by (3.1) and (3.2). Then

$$
(4.1) \qquad e \leq \eta
$$

and, assuming $x_{\Gamma} + R \leq L$,

$$
(4.2) \qquad \eta \leq Ce
$$

with a constant C depending only on the mesh regularity.

*Proof.* The proof of the upper bound (4.1) is completely standard: we set $e = u - u_h^c$,
observe that $e \in H_0^1(\Omega)$, use the weak formulation (2.3), Proposition 3.4 and integration by
parts (with the fact that $e = 0$ on $\partial\Omega$), to get

$$
\begin{align*}
\|\nabla u - \nabla u_h^c\|_{\Omega} &= \int_{\Omega} (\nabla u - \nabla u_h^c) \cdot (\nabla u - \nabla u_h^c) = \int_{\Omega} (\nabla u - \nabla u_h^c) \cdot \nabla e = \int_{\Omega} f e - \int_{\Omega} \nabla u_h^c \cdot \nabla e \\
&= \int_{\Omega} \operatorname{div} \sigma_h^c e - \int_{\Omega} \nabla u_h^c \cdot \nabla e = \int_{\Omega} (-\sigma_h^c - \nabla u_h^c) \cdot \nabla e \le \|\sigma_h^c + \nabla u_h^c\|_{\Omega} \|\nabla e\|_{\Omega},
\end{align*}
$$

where for the inequality we apply Cauchy-Schwartz inequality. Hence (4.1).

The proof of the lower bound (4.2) needs the following four steps. In this proof, we shall
use the letter *C* for various constants that depend only on the mesh regularity (in particular,
independent of the geometrical parameter *R* and of the polynomial degree *k*).

Step 1: error caused by the interface, prior to discretization. Let us begin with a “continuous” version of the coupled 0D/2D model:

$$
(4.4a) \qquad - \Delta \tilde{u} = f \text{ in } \tilde{\Omega},
$$

$$
(4.4b) \qquad \tilde{u} = \begin{cases} u' & \text{on } \Gamma, \\ 0 & \text{on } \Gamma_{\text{wall}}, \\ u_{\text{out}} & \text{on } \Gamma_{\text{out}}, \end{cases}
$$

and set

$$
u^c = \begin{cases} u' & \text{on } \Omega', \\ \tilde{u}, & \text{on } \tilde{\Omega}. \end{cases}
$$

We want to study $\|\nabla u - \nabla u^c\|_{\Omega}$ which represents the error introduced by the interface alone, without discretizing the problem in $\tilde{\Omega}$. We also introduce the continuous version of $\sigma_h^\Gamma$: $\sigma^\Gamma \in H_\Gamma(\text{div}, \omega^\Gamma)$, $p^\Gamma \in L^2(\omega^\Gamma)$ with $H_\Gamma(\text{div}, \omega^\Gamma) = \{\tau \in H(\text{div}, \omega^\Gamma) : \tau \cdot n = 0 \text{ on } \Gamma\}$, such that

$$
(4.5a) \qquad \int_{\omega^\Gamma} \sigma^\Gamma \cdot \tau^\Gamma + \int_{\omega^\Gamma} p^\Gamma \operatorname{div} \tau^\Gamma = - \int_{\omega^\Gamma} \nabla \tilde{u} \cdot \tau^\Gamma \qquad \forall \tau^\Gamma \in H_\Gamma(\operatorname{div}, \omega^\Gamma),
$$

$$
(4.5b) \qquad \int_{\omega^\Gamma} q^\Gamma \operatorname{div} \sigma^\Gamma = \int_{\omega^\Gamma} f q^\Gamma \qquad \forall q^\Gamma \in L^2(\omega^\Gamma).
$$
Let us prove

$$ (4.6) \qquad \|\sigma^{\Gamma} + \nabla \tilde{u}\|_{\omega^{\Gamma}} \le C \|\nabla u - \nabla u^{c}\|_{\Omega}. $$

Let $\omega_R^\Gamma = \text{supp}(\psi^\Gamma) \cap \tilde{\Omega}$ and $\Gamma_R = \partial\omega_R^\Gamma \cap \{x = x_\Gamma + R\}$. In Lemma A.1 we show that there exists $\theta \in H^1(\omega_R^\Gamma)$ such that

$$ (4.7) \qquad \begin{aligned} \Delta\theta &= 0 && \text{in } \omega_R^\Gamma, \\ \nabla\theta \cdot \mathbf{n} &= -\nabla\tilde{u} \cdot \mathbf{n} && \text{on } \Gamma, \\ \nabla\theta \cdot \mathbf{n} &= 0 && \text{on } \Gamma_R, \\ \theta &= 0 && \text{on } \tilde{\Gamma}_{\text{wall}} \cap \partial\omega_R^\Gamma, \end{aligned} $$

and $\|\nabla\theta\|_{\omega_R^\Gamma} \le C_1 \|\nabla\tilde{u} \cdot \mathbf{n}\|_{-1/2, \Gamma}$, with $C_1$ independent of $R$. Let

$$ \tau^c = \begin{cases} \nabla\theta & \text{in } \omega_R^\Gamma, \\ 0 & \text{on } \omega^\Gamma \setminus \omega_R^\Gamma. \end{cases} $$

By construction, $\text{div}\tau^c = 0$ in $\omega^\Gamma$. This implies

$$ \text{div}(\sigma^{\Gamma} + \nabla \tilde{u} + \tau^{c}) = 0 \quad \text{in } \omega^{\Gamma} $$

thanks to (4.5b), and

$$ (4.8) \qquad (\sigma^{\Gamma} + \nabla \tilde{u} + \tau^{c}) \cdot \mathbf{n} = 0 \quad \text{on } \Gamma $$

since $\sigma^\Gamma \cdot \mathbf{n} = 0$ and $\tau^c \cdot \mathbf{n} = -\nabla \tilde{u} \cdot \mathbf{n}$ by (4.7) on $\Gamma$. This allows us to take $\tau^\Gamma = \sigma^\Gamma + \nabla \tilde{u} + \tau^c$ as a test function in (4.5a). Hence, by Cauchy-Schwartz,

$$ (4.9) \qquad \|\sigma^{\Gamma} + \nabla \tilde{u}\|_{\omega^{\Gamma}} \le \|\tau^{c}\|_{\omega^{\Gamma}} \le C_1 \|\nabla \tilde{u} \cdot \mathbf{n}\|_{-1/2, \Gamma}. $$

Coming back to the error, denoting by $[\cdot]$ the jump on $\Gamma$, we have that for all $v \in H_0^1(\Omega)$

$$ \int_{\Omega} \nabla(u - u^c) \cdot \nabla v = \int_{\Gamma} [\nabla(u - u^c) \cdot \mathbf{n}] v = \int_{\Gamma} \nabla \tilde{u} \cdot \mathbf{n} v, $$

where we have split the domain $\Omega$ in $\Omega'$ and $\tilde{\Omega}$, integrated by parts, using $\Delta(u - u^c) = 0$ on $\Omega'$ and $\tilde{\Omega}$, $v \equiv 0$ on $\partial\Omega$, and $\nabla u' \cdot \mathbf{n} = 0$ on $\Gamma$. Now, from Lemma A.2 we know that for each $w \in H^{1/2}(\Gamma)$ there exists $v \in H_0^1(\Omega)$ such that $v|_{\Gamma} = w$ and such that $\|\nabla v\|_{\Omega} \le C_2 \|w\|_{1/2, \Gamma}$; with $C_2$ independent of geometrical parameters. We can take $v = \theta$ in $\omega_R^\Gamma$, with $\theta$ the solution of Lemma A.2 and extend it to all the domain taking e.g. $v$ the mirror of $\theta$ in $\Omega'$. Thus

$$ \int_{\Gamma} \nabla \tilde{u} \cdot \mathbf{n} v = \int_{\Omega} \nabla(u - u^c) \cdot \nabla v \le \|\nabla(u - u^c)\|_{\Omega} \|\nabla v\|_{\Omega} \le C_2 \|w\|_{1/2, \Gamma} \|\nabla(u - u^c)\|_{\Omega} $$

that implies

$$ \|\nabla \tilde{u} \cdot \mathbf{n}\|_{-1/2, \Gamma} \le C_2 \|\nabla(u - u^c)\|_{\Omega}. $$

Coming back to (4.9) and using this latter inequality, we get the desired estimate (4.6) with $C = C_1 C_2$.
Step 2: error caused by the interface, adding the discretization in $\tilde{\Omega}$. Here we prove that

$$ (4.10) \qquad \|\sigma_h^\Gamma + (\nabla \tilde{u}_h) \psi_h^\Gamma\|_{\omega^\Gamma} \le C(\|\nabla(\tilde{u} - \tilde{u}_h)\|_{\omega^\Gamma} + \|\nabla(u - u^c)\|_{\Omega}). $$

This can be viewed as a discrete analogue of (4.6). The proof is based on Theorem 1.2 of [2], but the latter cannot be applied directly due to a mismatch in boundary conditions: $\sigma_h^\Gamma$ is required to vanish on both $\Gamma$ and $\tilde{\Gamma}$, while $\psi_h^\Gamma$ vanishes only on $\tilde{\Gamma}$. To circumvent this difficulty, we enlarge $\omega^\Gamma$ to $\omega^{\Gamma,m}$ by adding to $\omega^\Gamma$ its mirror image with respect to $\Gamma$. Similarly, we extend $\psi_h^\Gamma$, $\tilde{u}_h$, $\sigma_h^\Gamma$ from $\omega^\Gamma$ to $\psi_h^{\Gamma,m}$, $\tilde{u}_h^m$, $\sigma_h^{\Gamma,m}$ on $\omega^{\Gamma,m}$ as functions symmetric with respect to $\Gamma$. Let $\tilde{\Gamma}^m$ denote the mirror image of $\tilde{\Gamma}$. Theorem 1.2 of [2] can be formulated on $\omega^{\Gamma,m}$ as

$$ \begin{array}{ll} \min_{\mathbf{v}_h \in H(\mathrm{div}, \omega^{\Gamma,m}) \cap \mathrm{RT}_k} & \left\| \mathbf{v}_h + \psi_h^{\Gamma,m} \nabla \tilde{u}_h^m \right\|_{\omega^{\Gamma,m}} \le C(\omega^{\Gamma,m}, \psi_h^{\Gamma,m}) \min_{\substack{\mathbf{v} \in H(\mathrm{div}, \omega^{\Gamma,m}) \\ \mathrm{div} \, \mathbf{v} = f}} \left\| \mathbf{v} + \nabla \tilde{u}_h^m \right\|_{\omega^{\Gamma,m}} \\ \mathrm{div} \, \mathbf{v}_h = \psi_h^{\Gamma,m} f - \nabla \psi_h^{\Gamma,m} \cdot \nabla \tilde{u}_h^m & \mathbf{v}_h = 0 \text{ on } \tilde{\Gamma} \cup \tilde{\Gamma}^m \end{array} $$

where $C(\omega^{\Gamma,m}, \psi_h^{\Gamma,m}) \le C(\|\psi_h^{\Gamma,m}\|_{\infty} + C_P \|\nabla \psi_h^{\Gamma,m}\|_{\infty})$ with $C_P$ the Poincaré constant of the space $H^1(\omega^{\Gamma,m})$ under the constraint of functions vanishing on $\tilde{\Gamma} \cup \tilde{\Gamma}^m$, i.e. $\|v\|_{\omega^{\Gamma,m}} \le C_P \|\nabla v\|_{\omega^{\Gamma,m}}$ for all $v \in H^1(\omega^{\Gamma,m})$ with $v = 0$ on $\tilde{\Gamma} \cup \tilde{\Gamma}^m$ (note that these constraints are imposed on the part of $\partial\omega^{\Gamma,m}$ where $\psi_h^{\Gamma,m}$ vanishes). By our geometrical assumptions, $C_P$ is of order $R$ and $\|\nabla \psi_h^\Gamma\|_{\infty} \le \frac{C}{R}$ so that $C(\omega^{\Gamma,m}, \psi_h^{\Gamma,m}) \le C$. Comparing with the definition of $\sigma_h^\Gamma$ we see that the minimum on the left-hand side of (4.11) is attained on $\sigma_h^{\Gamma,m}$ (note in particular that $\sigma_h^{\Gamma,m} \cdot n = 0$ on $\Gamma$ by symmetry). In order to identify the minimum on the other side, we introduce $\hat{\sigma}^\Gamma \in H_\Gamma(\mathrm{div}, \omega^\Gamma)$, $\hat{p}^\Gamma \in L^2(\omega^\Gamma)$ such that

$$ (4.12a) \qquad \int_{\omega^\Gamma} \hat{\sigma}^\Gamma \cdot \tau^\Gamma + \int_{\omega^\Gamma} \hat{p}^\Gamma \mathrm{div} \tau^\Gamma = - \int_{\omega^\Gamma} \nabla \tilde{u}_h \cdot \tau^\Gamma \qquad \forall \tau^\Gamma \in H_\Gamma(\mathrm{div}, \omega^\Gamma), $$

$$ (4.12b) \qquad \int_{\omega^\Gamma} q^\Gamma \mathrm{div} \hat{\sigma}^\Gamma = \int_{\omega^\Gamma} f q^\Gamma \qquad \forall q^\Gamma \in L^2(\omega^\Gamma). $$

We see then that the minimum on the right-hand side of (4.11) is attained on $\hat{\sigma}^{\Gamma,m}$, which is the mirror extension of $\hat{\sigma}^\Gamma$ to $\omega^{\Gamma,m}$. Going back in (4.11) to the subdomain $\omega^\Gamma$ of $\omega^{\Gamma,m}$ and using the symmetry gives

$$ \|\sigma_h^\Gamma + (\nabla \tilde{u}_h) \psi_h^\Gamma\|_{\omega^\Gamma} \le C \|\hat{\sigma}^\Gamma + \nabla \tilde{u}_h\|_{\omega^\Gamma}. $$

This entails by the triangle inequality

$$ (4.13) \qquad \|\sigma_h^\Gamma + (\nabla \tilde{u}_h) \psi_h^\Gamma\|_{\omega^\Gamma} \le C(\|\hat{\sigma}^\Gamma - \sigma^\Gamma\|_{\omega^\Gamma} + \|\sigma^\Gamma + \nabla \tilde{u}\|_{\omega^\Gamma} + \|-\nabla \tilde{u} + \nabla \tilde{u}_h\|_{\omega^\Gamma}). $$

To bound the first term in (4.13), we take the difference between (4.5) and (4.12) setting $\tau^\Gamma = \sigma^\Gamma - \hat{\sigma}^\Gamma$. Noting that $\mathrm{div} \tau^c = 0$, this yields

$$ (4.14) \qquad \|\sigma^\Gamma - \hat{\sigma}^\Gamma\|_{\omega^\Gamma} \le \|\nabla \tilde{u} - \nabla \tilde{u}_h\|_{\omega^\Gamma}. $$

The second term in (4.13) is bounded by (4.6). Finally, (4.13) gives (4.10).
Step 3: discretization error inside $\tilde{\Omega}$. We have at all the nodes $a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma$

$$ (4.15) \qquad \|\sigma_h^a + (\nabla \tilde{u}_h)\psi^a\|_{\omega^a} \le C \|\nabla \tilde{u} - \nabla \tilde{u}_h\|_{\omega^a}. $$

This well known estimate follows, for example, from Theorem 1.2 of [2] applied on each patch $\omega^a$.

Step 4: putting everything together. From definition of $\sigma_h^c$, see (3.7), (3.8) (3.9), and the partition of unity (3.3) we obtain on every mesh element $T \in \mathcal{T}_h$

$$ \|\tilde{\sigma}_h + \nabla \tilde{u}_h\|_T \le \|\sigma_h^\Gamma + (\nabla \tilde{u}_h)\psi_h^\Gamma\|_T + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} (1 - \psi^\Gamma(a)) \|\sigma_h^a + (\nabla \tilde{u}_h)\psi^a\|_T $$

where $\sigma_h^\Gamma$ and $\sigma_h^a$ are extended by 0 outside of their domains of definitions $\omega^\Gamma$ and $\omega^a$ respectively. The number of non-zero terms in the sum above is thus uniformly bounded by a constant that depends only on the regularity of the mesh. Taking the squares on both sides of the inequality above leads to

$$ \|\tilde{\sigma}_h + \nabla \tilde{u}_h\|_T^2 \le C \left( \|\sigma_h^\Gamma + (\nabla \tilde{u}_h)\psi_h^\Gamma\|_T^2 + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} \|\sigma_h^a + (\nabla \tilde{u}_h)\psi^a\|_T^2 \right). $$

Taking the sum over $T \in \mathcal{T}_h$, noting $\sigma' + \nabla u' = 0$ on $\Omega'$, and then using the bounds (4.10), (4.15) leads to

$$ \begin{aligned} \|\sigma_h^c + \nabla u_h^c\|_{\Omega}^2 &\le C \left( \|\sigma_h^\Gamma + (\nabla \tilde{u}_h)\psi_h^\Gamma\|_{\omega^\Gamma}^2 + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} \|\sigma_h^a + (\nabla \tilde{u}_h)\psi^a\|_{\omega^a}^2 \right) \\ &\le C \left( \|\nabla(u-u^c)\|_{\Omega}^2 + \|\nabla \tilde{u} - \nabla \tilde{u}_h\|_{\omega^\Gamma}^2 + \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} \|\nabla \tilde{u} - \nabla \tilde{u}_h\|_{\omega^a}^2 \right) \\ &\le C \left( \|\nabla(u-u^c)\|_{\Omega}^2 + \|\nabla \tilde{u} - \nabla \tilde{u}_h\|_{\tilde{\Omega}}^2 \right) \end{aligned} $$

since the number of possible overlaps between the different patches $\omega^\Gamma$ and $\omega^a$ is uniformly bounded.

By integration by parts,

$$ \int_{\tilde{\Omega}} (\nabla u - \nabla \tilde{u}) \cdot (\nabla \tilde{u} - \nabla \tilde{u}_h) = - \int_{\tilde{\Omega}} (\Delta u - \Delta \tilde{u})(\tilde{u} - \tilde{u}_h) + \int_{\partial \tilde{\Omega}} (\nabla u - \nabla \tilde{u}) \cdot \mathbf{n}(\tilde{u} - \tilde{u}_h) = 0. $$

Hence

$$ \|\nabla(u - u_h^c)\|_{\Omega}^2 = \|\nabla(u - u')\|_{\Omega'}^2 + \|\nabla(u - \tilde{u}) + \nabla(\tilde{u} - \tilde{u}_h)\|_{\tilde{\Omega}}^2 = \|\nabla(u - u^c)\|_{\Omega}^2 + \|\nabla \tilde{u} - \nabla \tilde{u}_h\|_{\tilde{\Omega}}^2, $$

so that

$$ \|\sigma_h^c + \nabla u_h^c\|_{\Omega}^2 \le C \|\nabla(u - u_h^c)\|_{\Omega}, $$

i.e. (4.2).

<figure><img src="image_5.png" /></figure>
**5. Numerical results.** We report here the results obtained using the 0D/2D model (2.7)–(2.8) for the problem (2.1)–(2.2) in the domain presented in Figure 2 with $L = 5.1$, $R = 0.5$, $W = 0.9$, and $L_{out} = 0.8$. All the computations are performed in FreeFEM [8] and we use $\mathbb{P}_2$-Lagrange FE for the 2D model in $\tilde{\Omega}$. Since we do not dispose of an analytical solution to (2.1)–(2.2), we use a reference solution obtained with $\mathbb{P}_2$-Lagrange FEM on a fine quasi-uniform mesh on $\Omega$ with mesh size $h = 0.005$ for the tests in subsection 5.1, and on a very fine adapted mesh on $\Omega$ with 1540177 Degrees of Freedom (DoF) for the tests in subsection 5.2. The reference solution on the quasi-uniform mesh is shown in Figure 3.

**5.1. Numerical validation on quasi-uniform meshes.** In Figure 4, we report on a series of numerical experiments varying the interface position $x_{\Gamma}$ from 0.1 (near the inlet) to 5.08 (practically in the corner), and employing quasi-uniform meshes on $\tilde{\Omega}$, composed of almost equilateral triangles of approximately the same size. In Figures 4a, 4c, and 4e, we plot the total error $e$ (3.1) and the total estimator $\eta$ (3.2) vs. the interface positions $x_{\Gamma}$, using the meshes of maximal sizes $h \approx 0.08, 0.04, 0.02$ respectively. As expected, the “modeling” error caused by 0D model simplification is negligible over a wide range of interface positions. It becomes predominant only when the interface is placed very near the corner $x_{\Gamma} \to 5.1$. In Figures 4b, 4d, and 4f, we plot the total error $e$ and the total estimator $\eta$ with respect to different mesh sizes $h$ and for an interface position $x_{\Gamma}$ fixed respectively to 0.1, 4.6 and 5.08. We choose these values because they represent the three typical cases: i) the interface placed too far from the corner so that the meshed domain $\tilde{\Omega}$ is unnecessarily big; ii) an optimal position of the interface providing a good balance between the modeling error and the discretization error, thus minimizing the size of the meshed domain without compromising the overall accuracy; (iii) interface placed too close to the corner, i.e. in the range of $x_{\Gamma}$ values where the error/estimator grows up abruptly in the plots on Figures 4a, 4c, and 4e. We observe the convergence under the mesh refinement in the first 2 situations ($x_{\Gamma} = 0.1$ and 4.6), contrary to the last situation with $x_{\Gamma} = 5.08$, where the convergence is lost, due to the bad choice of interface. When the interface is well chosen, like in Figures 4b and 4d, the rate of convergence is sub-optimal with respect to what one would expect on a smooth benchmark solution. This is not surprising because the actual solution is singular near the reentrant corners. In the next subsection 5.2, we shall employ a doubly adaptive strategy (for interface position and for the mesh), which will enable us to restore an optimal rate of convergence.

All the plots of Figure 4 confirm that the estimator $\eta$ provides indeed an upper bound for the error $e$ and both quantities are always of the same order of magnitude, even if $x_{\Gamma} + R > L$, cf. Remark 3.1. A study of the index of efficiency that indicates the optimality of the estimator will be carried out in the next subsection in the case of adapted meshes.

**5.2. Adaptive algorithm.** We now inspire ourselves from the construction of the estimator (3.2) to detect a suitable interface position $x_{\Gamma}$. Guided by the proof of the upper and lower bounds for the error (Steps 2 and 3 of the proof of Theorem 4.1 in particular), we conjecture that the quantity

$$ (5.1) \qquad \eta^{\Gamma} := \|\sigma_h^{\Gamma} + (\nabla \tilde{u}_h)\psi_h^{\Gamma}\|_{\omega^{\Gamma}}, $$

can indicate the modelling error, caused by the replacement of the 2D model by the 0D one, and can thus drive the choice of the interface position. A theoretical motivation for this
<figure><img src="image_6.png" /><figcaption>(a) e, η vs. x<sub>Γ</sub>, h = 0.0789.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) e, η vs. x<sub>Γ</sub>, h = 0.0410.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) e, η vs. x<sub>Γ</sub>, h = 0.02159.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) e and η vs. h, x<sub>Γ</sub> = 0.1.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) e and η vs. h, x<sub>Γ</sub> = 4.6.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(f) e and η vs. h, x<sub>Γ</sub> = 5.08.</figcaption></figure>

Figure 4: Comparison between error *e* and error estimator *η*: with respect to the interface position *x*<sub>Γ</sub> for a fixed uniform mesh with mesh size *h*, see Figures 4a, 4c, and 4e, and with respect to the mesh size *h* once the interface position *x*<sub>Γ</sub> is fixed, see Figures 4b, 4d, and 4f.

conjecture stems from the estimate (4.10) established in the proof of our main theorem. ηΓ is
bounded there by two terms, the second one being precisely the part of the error caused by
the introduction of the coupling interface. This contribution is global. The first contribution
<figure><img src="image_12.png" /><figcaption>Figure 5: Evolution of $\eta^\Gamma$ with respect to the interface position $x_\Gamma$ and a fixed mesh size $h$ and with respect to the mesh refinement with a fixed $x_\Gamma$.</figcaption></figure>

in (4.10) represents a discretization error, locally on the patch $\omega^\Gamma$. This local contribution should be negligible in most practical cases, in comparison with the global one. A numerical evidence for the pertinence of $\eta^\Gamma$ is provided in Figure 5. We plot there $\eta^\Gamma$ vs. the interface position (using quasi-uniform meshes of the same size), and $\eta^\Gamma$ vs. the mesh size, once the interface position has been fixed. We observe that $\eta^\Gamma$ is indeed almost independent of the mesh refinement, and the range of interface positions where it grows abruptly coincides with the similar region for the actual error, cf. Figures 4a and 5a. We conclude that $\eta^\Gamma$ can be used to detect the part of the error caused by the interface position.

Once the interface is fixed, it remains to refine the mesh on $\tilde{\Omega}$ locally. We use for this the error indicators defined on each mesh cell as

$$ (5.2) \qquad \eta_T^D := \|(\sigma_h^D + \nabla \tilde{u}_h)|_T\|_T \quad \text{for all } T \in \mathcal{T}_h, \text{ where } \sigma_h^D := \sum_{a \in \mathcal{V}_h \setminus \mathcal{V}_h^\Gamma} \sigma_h^a $$

(the superscript “D” here refers to the “discretization” part of the error). These are the usual equilibrated flux error indicators, normally used for a problem with Dirichlet boundary conditions everywhere. This seems a reasonable choice since the “modeling” error is supposed to be controlled at this stage by a proper choice of the coupling interface.

Our approach is summarised in Algorithm 5.1. We choose first a suitable interface using the error indicator $\eta^\Gamma$ (5.1) in the loop in Line 3, employing coarse quasi-uniform meshes, cf. Line 2. The goal is to find an interface position such that the error indicator $\eta^\Gamma$ (5.1) is smaller than the prescribed tolerance $\text{tol}^\Gamma$, which is in turn sufficiently smaller than the desired tolerance $\text{tol}$ for the overall error. At the second stage of the Algorithm, Line 8, we perform the mesh adaptation loop without ever moving the interface again. Note that the mesh adaptation is guided by the local error indicators $\{\eta_T^D\}_{T \in \mathcal{T}_h}$ from (5.2), while the overall error is estimated by the global error indicator $\eta$ from our main Theorem, which controls effectively both the modeling and discretization errors. The details of our approach to the
mesh adaptation (i.e. constructing the new mesh using the error indicators on the current mesh) are given in Appendix B.

<caption>**Algorithm 5.1** Adaptive interface position and mesh refinement</caption>

**Require:** A tolerance for the global error tol, a tolerance for the error due to the interface position (modelling error) tol<sup>Γ</sup>, a step to move the interface δx > 0.
1: Set the interface position at coordinate $x_Γ = L - R$, cf. Figure 2
2: Construct a coarse mesh on $\tilde{\Omega}$, compute the coupled solution (2.8) and $\eta^Γ$, cf. (5.1)
3: **while** $\eta^Γ > \text{tol}^Γ$ and $x_Γ > 0$ **do**
4:     Redefine $x_Γ := x_Γ - \delta x$
5:     Construct a coarse mesh on the new $\tilde{\Omega}$, compute the coupled solution and $\eta^Γ$
6: **end while**
7: Compute $\eta$, as in (3.2), and $\eta_T^D$, as in (5.2)
8: **while** $\eta > \text{tol}$ **do**
9:     Construct a new mesh on $\tilde{\Omega}$ based on the error indicators $\eta_T^D$ (cf. Appendix B)
10:    Compute the coupled solution $u_h^c$, $\eta$, and $\eta_T^D$ on the new mesh
11: **end while**
12: **return** Coupled solution and $x_Γ$

We have tested Algorithm 5.1 for two values of the desired tolerance: tol = 1e-2 and tol = 1e-4. The other parameters were set to δx = 0.1 and tolΓ = 0.1 tol (numerical experiments with other choices of tolΓ, ex. tolΓ = 0.5 tol, give similar results but are not reported here). As expected, the Algorithm puts the interface further from the corner when we decrease the tolerance: xΓ = 4.4 for tol = 1e-2 and xΓ = 3.8 for tol = 1e-4. Convergence of the error and the estimator with respect to the number of DoF on the iterations of the mesh adaptation loop (once the interface has been chosen) is presented in Table 1 and Figure 6. We report in Table 1, on subsequent adapted meshes, the numbers of DoF, the error e, the estimator η, and the effectivity index I = η/e. The algorithm needed 4 iterations with tol = 1e-2 and 7 iterations for tol = 1e-4, with I close to 1 in both cases. We also observe in Figure 6 that the error and the estimator converge optimally with respect to the number of DoF, i.e. at the rate that would be expected on quasi-uniform meshes if the solution were smooth.

In Figure 7, we plot the initial coarse mesh and one of the adaptive meshes constructed
by (Appendix B) with tol = 1e-4. As expected, the algorithm refines mostly near the corner
where the solution is singular. In Figure 8, we present the distributions of the error and of
the estimator on the same adapted mesh given in Figure 7. We observe that they are locally
equivalent, confirming that our choice of the local estimators is reasonable. At Figure 7b,
we also observe that the error is sufficiently well equi-distributed on mesh cells (the error of
about 1e-5 per cell) with the exception of the cells near the reentrant corners, where more
refinement is clearly needed. In particular, this confirms that the relatively fine meshes in the
outer corners are indeed a reasonable choice at this iteration of the mesh adaptation.

**6. Conclusions.** In this paper we have proposed a coupled approach to solve a 2D Poisson problem in a domain with long channels, splitting the domain into a portion where the solution
<table><caption>Table 1: Results of Algorithm 5.1 with tol = 1e-2 and tol = 1e-4.</caption><thead><tr><th>DoF</th><th>e</th><th>η</th><th>I</th></tr></thead><tbody><tr><td>84</td><td>1.4968e-1</td><td>1.7087e-1</td><td>1.1</td></tr><tr><td>302</td><td>5.1694e-2</td><td>6.5485e-2</td><td>1.3</td></tr><tr><td>417</td><td>2.9590e-2</td><td>3.5823e-2</td><td>1.2</td></tr><tr><td>814</td><td>1.3262e-2</td><td>1.5647e-2</td><td>1.2</td></tr><tr><td>2395</td><td>4.7947e-3</td><td>5.6709e-3</td><td>1.2</td></tr></tbody></table><p>(a) tol = 1e-2, x<sub>Γ</sub> = 4.4.</p><table><thead><tr><th>DoF</th><th>e</th><th>η</th><th>I</th></tr></thead><tbody><tr><td>98</td><td>1.5490e-1</td><td>1.7616e-1</td><td>1.1</td></tr><tr><td>347</td><td>5.3231e-2</td><td>6.6758e-2</td><td>1.3</td></tr><tr><td>512</td><td>2.6044e-2</td><td>3.1548e-2</td><td>1.2</td></tr><tr><td>1269</td><td>1.0420e-2</td><td>1.2632e-2</td><td>1.2</td></tr><tr><td>2721</td><td>3.8437e-3</td><td>4.5794e-3</td><td>1.2</td></tr><tr><td>8388</td><td>1.2095e-3</td><td>1.4192e-3</td><td>1.2</td></tr><tr><td>27531</td><td>3.3055e-4</td><td>3.9007e-4</td><td>1.2</td></tr><tr><td>106489</td><td>8.42367e-5</td><td>9.9955e-5</td><td>1.2</td></tr></tbody></table><p>(b) tol = 1e-4, x<sub>Γ</sub> = 3.8.</p>

<figure><img src="image_13.png" /><figcaption>Figure 6: Evolution of e and η vs. the number of DoF for a fixed $x_Γ$.</figcaption></figure>

is approximated by a simple analytical expression and another one where it is approximated numerically. We have developed an *a posteriori* error estimator that, for a given tolerance on the error, allows us both to choose the position of the coupling interface and to adapt the mesh. Our theoretical results, namely that the estimator is guaranteed and efficient, are confirmed and enhanced by numerical experiments on automatically adapted meshes. Some partial results for the Stokes equation in a similar setting are available in [1]. We hope to adapt this approach to more complex governing equations in a future work.

**Appendix A. Some technical Lemmas.** We recall here two well-known lemmas, needed for the proof of Theorem 4.1, and give their proofs for completeness of exposition. Let $H_{00}^s(\Gamma)$, with $s \in \{-\frac{1}{2}, \frac{1}{2}\}$ be the spaces of functions (distributions) on $\Gamma$ of the form

$$ (A.1) \qquad \eta = \sum_{k \ge 1} \eta_k \sin \left( \frac{k\pi}{R} y \right), $$
<figure><img src="image_14.png" /><figcaption>Figure 7: Two meshes on $\tilde{\Omega}$ constructed by Algorithm 5.1 with tol = 1e-4.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 8: Distributions of the error and the estimator on iteration 5 of Algorithm 5.1 with tol = 1e-4.</figcaption></figure>
with the norm $||\eta||_{s,\Gamma} = (\sum_{k \ge 1} \frac{R}{2} \eta_k^2 (\frac{k\pi}{R})^{2s})^{1/2}$.

**Lemma A.1.** For any $\eta \in H_{00}^{-1/2}(\Gamma)$, let $\theta$ be the solution to

$$ \Delta\theta = 0 \qquad \text{in } \omega_R^\Gamma, $$

$$ \nabla\theta \cdot \mathbf{n} = \eta \qquad \text{on } \Gamma, $$

$$ \nabla\theta \cdot \mathbf{n} = 0 \qquad \text{on } \Gamma_R, $$

$$ \theta = 0 \qquad \text{on } \tilde{\Gamma}_{\text{wall}} \cap \partial\omega_R^\Gamma. $$

Then $||\nabla\theta||_{\omega_R^\Gamma} \le C_1 ||\eta||_{-1/2,\Gamma}$ with $C_1 > 0$ which does not depend on $R$.

*Proof.* For any $\eta \in H_{00}^{-1/2}(\Gamma)$ written as (A.1), by direct calculation, $\theta$ is given by

$$ \theta = \sum_{k \ge 1} \eta_k \sin \left( \frac{k\pi}{R} y \right) \frac{\cosh \left( \frac{k\pi}{R} (x - x_{\Gamma} - R) \right)}{\frac{k\pi}{R} \sinh(k\pi)}. $$

Thus, denoting $x_R^\Gamma = x_\Gamma + R$,

$$ \begin{aligned} ||\nabla\theta||_{\omega_R^\Gamma}^2 &= \sum_{k \ge 1} \left( \frac{\eta_k}{\sinh(k\pi)} \right)^2 \times \\ & \qquad \left( \left\| \sin\left(\frac{k\pi}{R}y\right) \sinh\left(\frac{k\pi}{R}(x-x_\Gamma - R)\right) \right\|_{\omega_R^\Gamma}^2 + \left\| \cos\left(\frac{k\pi}{R}y\right) \cosh\left(\frac{k\pi}{R}(x-x_\Gamma^R)\right) \right\|_{\omega_R^\Gamma}^2 \right) \\ &= \sum_{k \ge 1} \eta_k^2 \frac{R^2}{2k\pi \tanh(k\pi)} \le \frac{1}{\tanh(\pi)} \sum_{k \ge 1} \eta_k^2 \frac{R^2}{2k\pi} := C_1^2 ||\eta||_{-1/2,\Gamma}^2, \end{aligned} $$

since $\tanh(k\pi) \ge \tanh(\pi)$ for $k \ge 1$.

**Lemma A.2.** For any $\eta \in H_{00}^{1/2}(\Gamma)$, let $\theta$ be the solution to

$$ \Delta\theta = 0 \quad \text{in } \omega_R^\Gamma, $$

$$ \theta = \eta \quad \text{on } \Gamma, $$

$$ \theta = 0 \quad \text{on } \Gamma_R \cup (\tilde{\Gamma}_{\text{wall}} \cap \partial\omega_R^\Gamma). $$

Then $||\nabla\theta||_{\omega_R^\Gamma} \le C_2 ||\eta||_{-1/2,\Gamma}$ with $C_2 > 0$ which does not depend on $R$.

*Proof.* For any $\eta \in H_{00}^{1/2}(\Gamma)$ written as (A.1), by direct calculation the solution $\theta$ is given by

$$ \theta = \sum_{k \ge 1} \eta_k \sin \left( \frac{k\pi}{R} y \right) \frac{\sinh \left( \frac{k\pi}{R} (x_\Gamma + R - x) \right)}{\sinh(k\pi)}. $$

Calculations similar to those in the proof of Lemma A.1 lead to the desired result with $C_2 = C_1$ of the previous lemma.
**Appendix B. Mesh adaptation in FreeFEM.** The traditional strategy for the mesh adaptation (see, for example, [14, Chapter 2] for a review) consists in marking certain mesh cells to refine according to the error indicator, and then splitting them in smaller cells. The library FreeFEM, used in our numerical experiments, does not provide tools for such a procedure: one cannot ask it to split only some specific cells, without touching the others. FreeFEM proposes instead the function **adaptmesh** that creates an entirely new mesh with the mesh size prescribed (approximately) at every point of the entire computation domain (it is also possible to generate anisotropic meshes, but we deal here only with the isotropic version, option **IsMetric=1**, where all the mesh cells are close to equilateral triangles).

We describe here an alternative mesh adaptation strategy from [7], which we name “hopt” (for “h optimal”). We start from the following (admittedly not always realistic) assumption: the error between the exact solution $\tilde{u}$ and the approximate solution $\tilde{u}_h$ on any mesh $\mathcal{T}_h$ of $\tilde{\Omega}$ is approximately given by

$$ (B.1) \qquad \|\tilde{u} - \tilde{u}_h\|_{1,\tilde{\Omega}}^2 \approx \int_{\tilde{\Omega}} h^{2\delta}(x) c^2(x) \, dx $$

where $h(x)$ is the mesh size distribution, i.e. $h(x) = h_T$ on any $T \in \mathcal{T}_h$, $\delta$ is a fixed parameter, and $c(x)$ is some a priori unknown function. Note also that the number of DoF is approximately given in the 2D case by

$$ N_{\text{DoF}} \sim \int_{\tilde{\Omega}} \frac{dx}{h^2(x)}, $$

since a regular triangle of diameter $h$ occupies the area of order $h^2$. Let us imagine momentarily that we know $c(x)$ and we want to construct an optimal mesh with the minimal possible $N_{\text{DoF}}$ to achieve a given error $\epsilon$, i.e. $\|u - \tilde{u}_h\|_{1,\tilde{\Omega}} = \epsilon$. This is a constrained minimization problem for the mesh size distribution $h(x)$, i.e. minimize $\int_{\tilde{\Omega}} \frac{dx}{h^2(x)}$ under $\int_{\tilde{\Omega}} h^{2\delta}(x) c^2(x) \, dx = \epsilon^2$, which gives the following optimal mesh size distribution

$$ h_{\text{opt}}(x) = \frac{\epsilon^{1/\delta}}{(\int_{\tilde{\Omega}} c^{2/(\delta+1)}(x) \, dx)^{1/(2\delta)}} \frac{1}{(c(x))^{1/(\delta+1)}}. $$

Of course, $c(x)$ is not known in practice, but, on a given mesh $\tilde{\mathcal{T}}_h$, we have the estimators $\eta_T^D$, see (5.2). It seems thus reasonable to expect that

$$ \|\tilde{u} - \tilde{u}_h\|_{1,\tilde{\Omega}}^2 \approx (\eta^D)^2 := \sum_{T \in \tilde{\mathcal{T}}_h} (\eta_T^D)^2 = \sum_{T \in \tilde{\mathcal{T}}_h} \|\sigma_h^D + \nabla \tilde{u}_h\|_T^2. $$

Reinterpreting the error in the form (B.1), and localizing to each triangle of the current mesh $\tilde{\mathcal{T}}_h$, suggests

$$ \int_T h^{2\delta}(x) c^2(x) \, dx \sim (\eta_T^D)^2, \quad \forall T \in \tilde{\mathcal{T}}_h. $$

We can thus approximate $c(x)$ on any triangle $T \in \tilde{\mathcal{T}}_h$ by $c(x) \approx \frac{\eta_T^D}{h_T^\delta \sqrt{|T|}}$, which gives

$$ (B.2) \qquad h_{\text{opt}}(x) = \frac{\epsilon^{1/\delta}}{\left(\sum_{T \in \tilde{\mathcal{T}}_h} (\eta_T^D)^{2/(\delta+1)} h_T^{-2\delta/(\delta+1)} |T|^{\delta/(\delta+1)}\right)^{1/(2\delta)}} \frac{h_T^{\delta/(\delta+1)} |T|^{1/(2\delta+2)}}{(\eta_T^D)^{1/(\delta+1)}} \quad \text{for } x \in T. $$
In practice, given a mesh $\tilde{T}_h$, we want to reduce the error $R_e$ times (with $R_e > 1$ a parameter to be specified). We thus put $\epsilon = \eta^D/R_e$ in (B.2) and give the resulting $h_{opt}$ mesh size distribution to the mesh-generating function of FreeFEM. This is the procedure we have used for mesh adaptation in Algorithm 5.1. There are two parameters ($\delta$ and $R_e$) to be specified. Based on some numerical experimentation (not reported here, but given in [1]) we have chosen $\delta = 1$ and $R_e = 4$.

**Acknowledgments.** We are grateful to Fei Gao for giving the initial impetus to this work and for sharing with us some of his expertise on fuel cells, and to Martin Vohralík for his interest in this work and for several stimulating and enlightening discussions that has helped us to improve it.

## REFERENCES

[1] H. ALBAZZAL, *Simulations efficaces de piles à combustible à l'aide de l'analyse d'erreur a posteriori*, PhD thesis, Univ. Bourgogne Franche-Comté, 2023.

[2] A. ERN, I. SMEARS, AND M. VOHRALÍK, *Discrete p-robust H(div)-lifting s and a posteriori estimates for elliptic problems with $H^{-1}$ source terms*, Calcolo, 54 (2017), pp. 1009–1025.

[3] A. ERN AND M. VOHRALÍK, *Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous galerkin, and mixed discretizations*, SIAM Journal on Numerical Analysis, 53 (2015), pp. 1058–1081.

[4] L. FORMAGGIA, J.-F. GERBEAU, F. NOBILE, AND A. QUARTERONI, *On the coupling of 3d and 1d navier-stokes equations for flow problems in compliant vessels*, Computer methods in applied mechanics and engineering, 191 (2001), pp. 561–582.

[5] L. FORMAGGIA, F. NOBILE, A. QUARTERONI, AND A. VENEZIANI, *Multiscale modelling of the circulatory system: a preliminary analysis*, Computing and visualization in science, 2 (1999), pp. 75–83.

[6] J.-F. GERBEAU AND B. PERTHAME, *Derivation of Viscous Saint-Venant System for Laminar Shallow Water; Numerical Validation*, Research Report RR-4084, INRIA, 2000. Projet M3N.

[7] O. GORYNINA, *Eléments finis adaptatifs pour l'équation des ondes instationnaire*, PhD thesis, Univ. Bourgogne Franche-Comté, 2018.

[8] F. HECHT, *New development in freefem++*, J. Numer. Math., 20 (2012), pp. 251–265, https://freefem.org/.

[9] S. KARVONEN, T. HOTTINEN, J. SAARINEN, AND O. HIMANEN, *Modeling of flow field in polymer electrolyte membrane fuel cell*, Journal of power sources, 161 (2006), pp. 876–884.

[10] E. MIGLIO, S. PEROTTO, AND F. SALERI, *Model coupling techniques for free-surface flow problems: Part i*, Nonlinear Analysis: Theory, Methods & Applications, 63 (2005), pp. e1885–e1896.

[11] G. PANASENKO, *Method of asymptotic partial decomposition of domain*, Mathematical Models and Methods in Applied Sciences, 8 (1998), pp. 139–156.

[12] A. QUARTERONI AND L. FORMAGGIA, *Mathematical modelling and numerical simulation of the cardiovascular system*, Handbook of numerical analysis, 12 (2004), pp. 3–127.

[13] M. TAYACHI, A. ROUSSEAU, E. BLAYO, N. GOUTAL, AND V. MARTIN, *Design and analysis of a schwarz coupling method for a dimensionally heterogeneous problem*, International Journal for Numerical Methods in Fluids, 75 (2014), pp. 446–465.

[14] R. VERFÜRTH, *A posteriori error estimation techniques for finite element methods*, OUP Oxford, 2013.

[15] C.-Y. WEN, A. D. LE, K.-T. JENG, AND B.-T. TSAI, *A numerical model of the cathode of a proton exchange membrane fuel cell with experimental validation*, International Journal Of Numerical Analysis And Modeling, Series B, 1 (2010), pp. 123–146.

[16] D. ZHOU, F. GAO, E. BREAZ, A. RAVEY, AND A. MIRAOUI, *Tridiagonal matrix algorithm for real-time simulation of a two-dimensional pem fuel cell model*, IEEE Transactions on Industrial Electronics, 65 (2017), pp. 7106–7118.