# Sharp convergence rates for the homogenization of the Stokes equations in a perforated domain

Loïc Balazi, Grégoire Allaire, Pascal Omnes

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

## To cite this version:

Loïc Balazi, Grégoire Allaire, Pascal Omnes. Sharp convergence rates for the homogenization of the Stokes equations in a perforated domain. 2024. hal-04541828v2

HAL Id: hal-04541828

https://hal.science/hal-04541828v2

Preprint submitted on 7 May 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sharp convergence rates for the homogenization of the Stokes equations in a perforated domain

Loïc Balazi¹,², Grégoire Allaire¹, Pascal Omnes²,³

May 7, 2024

¹ CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France.

² Université Paris-Saclay, CEA, Service de Génie Logiciel pour la Simulation, 91191, Gif-sur-Yvette, France.

³ Université Sorbonne Paris Nord, LAGA, CNRS UMR 7539, Institut Galilée, 99 Av. J.-B. Clément, 93430, Villetaneuse, France.

## Abstract

This paper is concerned with the homogenization of the Stokes equations in a periodic perforated domain. The homogenized model is known to be Darcy's law in the full domain. We establish a sharp convergence rate $O(\sqrt{\varepsilon})$ for the energy norm of the difference of the velocities, where $\varepsilon$ represents the size of the solid obstacles. This is achieved by using a two-scale asymptotic expansion of the Stokes equations and a new construction of a cut-off function which avoids the introduction of boundary layers. The main novelty is that our analysis applies for the physically relevant case of a porous medium where each of the fluid and solid parts is a connected subdomain.

## 1 Introduction

In this paper, we deal with the Stokes problem for an incompressible viscous fluid in a porous medium which is represented as a periodic perforated domain (see for example Figure 1 or Figure 2). This work is motivated by the derivation of an error estimate for a Multiscale Finite Element Method (MsFEM) applied to the Stokes equations in perforated domains [11, 7] which is based on the homogenization theory. Let $\Omega \subset \mathbb{R}^d$ be a regular bounded open set. As shown in Figure 1, we divide the domain $\Omega$ into a fixed solid part $B^\varepsilon$ and its complementary fluid part $\Omega^\varepsilon$. Here $\varepsilon$ denotes a small parameter equal to the ratio between the characteristic length of the periodic heterogeneities and the characteristic length of the domain. A first typical example of such a porous domain $\Omega^\varepsilon$ is given by Figure 1 where the solid obstacles $B^\varepsilon$ are a collection of isolated and periodically repeated obstacles. A second typical case, in dimension $d \ge 3$, is a porous domain $\Omega^\varepsilon$ perforated by a regular lattice of interconnected solid cylinders as presented in Figure 2. A precise definition of $\Omega^\varepsilon$ is given in Section 2.
# Sharp convergence rates for the homogenization of the Stokes equations in a perforated domain

Loïc Balazi¹,², Grégoire Allaire¹, Pascal Omnes²,³

May 7, 2024

¹ CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France.

² Université Paris-Saclay, CEA, Service de Génie Logiciel pour la Simulation, 91191, Gif-sur-Yvette, France.

³ Université Sorbonne Paris Nord, LAGA, CNRS UMR 7539, Institut Galilée, 99 Av. J.-B. Clément, 93430, Villetaneuse, France.

## Abstract

This paper is concerned with the homogenization of the Stokes equations in a periodic perforated domain. The homogenized model is known to be Darcy's law in the full domain. We establish a sharp convergence rate $O(\sqrt{\varepsilon})$ for the energy norm of the difference of the velocities, where $\varepsilon$ represents the size of the solid obstacles. This is achieved by using a two-scale asymptotic expansion of the Stokes equations and a new construction of a cut-off function which avoids the introduction of boundary layers. The main novelty is that our analysis applies for the physically relevant case of a porous medium where each of the fluid and solid parts is a connected subdomain.

## 1 Introduction

In this paper, we deal with the Stokes problem for an incompressible viscous fluid in a porous medium which is represented as a periodic perforated domain (see for example Figure 1 or Figure 2). This work is motivated by the derivation of an error estimate for a Multiscale Finite Element Method (MsFEM) applied to the Stokes equations in perforated domains [11, 7] which is based on the homogenization theory. Let $\Omega \subset \mathbb{R}^d$ be a regular bounded open set. As shown in Figure 1, we divide the domain $\Omega$ into a fixed solid part $B^\varepsilon$ and its complementary fluid part $\Omega^\varepsilon$. Here $\varepsilon$ denotes a small parameter equal to the ratio between the characteristic length of the periodic heterogeneities and the characteristic length of the domain. A first typical example of such a porous domain $\Omega^\varepsilon$ is given by Figure 1 where the solid obstacles $B^\varepsilon$ are a collection of isolated and periodically repeated obstacles. A second typical case, in dimension $d \ge 3$, is a porous domain $\Omega^\varepsilon$ perforated by a regular lattice of interconnected solid cylinders as presented in Figure 2. A precise definition of $\Omega^\varepsilon$ is given in Section 2.
<figure><img src="image_1.png" /><figcaption>Figure 1: Porous medium $\Omega^\epsilon$, obtained from a domain $\Omega$ perforated by a set of solid isolated obstacles $B^\epsilon$</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: Regular lattice of interconnected cylinders</figcaption></figure>

The steady-state Stokes problem with homogeneous Dirichlet boundary condition in a perforated domain is to find the velocity $u_\epsilon : \Omega^\epsilon \to \mathbb{R}^d$ and the pressure $p_\epsilon : \Omega^\epsilon \to \mathbb{R}$ solution to:

$$ egin{cases} -\nu \Delta u_{\epsilon} + \nabla p_{\epsilon} = \mathbf{f} & \text{in } \Omega^{\epsilon}, \\ \text{div } u_{\epsilon} = 0 & \text{in } \Omega^{\epsilon}, \\ u_{\epsilon} = \mathbf{0} & \text{on } \partial B^{\epsilon} \cap \partial \Omega^{\epsilon}, \\ u_{\epsilon} = \mathbf{0} & \text{on } \partial \Omega \cap \partial \Omega^{\epsilon}, \end{cases} \quad (1) $$

where $\nu > 0$ is the viscosity and $\mathbf{f}$ the applied force. Vector valued functions are written in bold. In the literature, some authors consider a different scaling of (1), replacing the velocity $u_\epsilon$ by $u_\epsilon = \epsilon^2 \tilde{u}_\epsilon$. This does not change anything to our methodology and our results, up to this $\epsilon^2$ factor.

We denote by $(\cdot, \cdot)$ the usual scalar product in $L^2(\Omega^\epsilon)$ for scalar and vector-valued functions. We introduce the classical velocity space $V = H_0^1(\Omega^\epsilon)^d = \{u \in H^1(\Omega^\epsilon)^d \text{ s.t. } u|_{\partial\Omega^\epsilon} = \mathbf{0}\}$ and pressure space $M = L_0^2(\Omega^\epsilon) = \{q \in L^2(\Omega^\epsilon), \text{ s.t. } \int_{\Omega^\epsilon} q = 0\}$. In $V$ we shall often use the semi-norm $|\cdot|_{H^1}$, equivalent to the usual norm and defined by:

$$ |v|_{H^1(\Omega^\epsilon)} = \|\nabla v\|_{L^2(\Omega^\epsilon)}, \quad v \in H^1(\Omega^\epsilon)^d. $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Porous medium $\Omega^\epsilon$, obtained from a domain $\Omega$ perforated by a set of solid isolated obstacles $B^\epsilon$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Regular lattice of interconnected cylinders</figcaption></figure>

The steady-state Stokes problem with homogeneous Dirichlet boundary condition in a perforated domain is to find the velocity $u_\epsilon: \Omega^\epsilon \to \mathbb{R}^d$ and the pressure $p_\epsilon: \Omega^\epsilon \to \mathbb{R}$ solution to:

$$ egin{cases} -\nu \Delta u_{\epsilon} + \nabla p_{\epsilon} = \mathbf{f} & \text{in } \Omega^{\epsilon}, \\ \operatorname{div} u_{\epsilon} = 0 & \text{in } \Omega^{\epsilon}, \\ u_{\epsilon} = \mathbf{0} & \text{on } \partial B^{\epsilon} \cap \partial \Omega^{\epsilon}, \\ u_{\epsilon} = \mathbf{0} & \text{on } \partial \Omega \cap \partial \Omega^{\epsilon}, \end{cases} \quad (1) $$

where $\nu > 0$ is the viscosity and $\mathbf{f}$ the applied force. Vector valued functions are written in bold. In the literature, some authors consider a different scaling of (1), replacing the velocity $u_\epsilon$ by $u_\epsilon = \epsilon^2 \tilde{u}_\epsilon$. This does not change anything to our methodology and our results, up to this $\epsilon^2$ factor.

We denote by $(\cdot, \cdot)$ the usual scalar product in $L^2(\Omega^\epsilon)$ for scalar and vector-valued functions. We introduce the classical velocity space $V = H_0^1(\Omega^\epsilon)^d = \{u \in H^1(\Omega^\epsilon)^d \text{ s.t. } u|_{\partial\Omega^\epsilon} = \mathbf{0}\}$ and pressure space $M = L_0^2(\Omega^\epsilon) = \{q \in L^2(\Omega^\epsilon), \text{ s.t. } \int_{\Omega^\epsilon} q = 0\}$. In $V$ we shall often use the semi-norm $|\cdot|_{H^1}$, equivalent to the usual norm and defined by:

$$ |v|_{H^1(\Omega^\epsilon)} = \|\nabla v\|_{L^2(\Omega^\epsilon)}, \quad v \in H^1(\Omega^\epsilon)^d. $$
We introduce the bilinear forms $a: V \times V \to \mathbb{R}$ and $b: M \times V \to \mathbb{R}$ by:

$$a(\mathbf{u}, \mathbf{v}) = \nu(\nabla \mathbf{u}, \nabla \mathbf{v}), \quad b(p, \mathbf{v}) = (p, \operatorname{div} \mathbf{v}).$$

Then, assuming that $\mathbf{f} \in L^2(\Omega)^d$, a weak formulation of the Stokes problem (1) reads as follows: find $\mathbf{u}_\varepsilon \in V$ and $p_\varepsilon \in M$ such that

$$\left\{ \begin{array}{l} a(\mathbf{u}_\varepsilon, \mathbf{v}) + b(p_\varepsilon, \mathbf{v}) = (\mathbf{f}, \mathbf{v}) \quad \forall \mathbf{v} \in V, \\ b(q, \mathbf{u}_\varepsilon) = 0 \quad \forall q \in M. \end{array} \right. \qquad (2)$$

It is well known that, if $\Omega^\varepsilon$ is connected (to ensure that the zero-average condition of the pressure is enough to remove the undetermined constant of the pressure), there exists a unique weak solution to (2) [9].

The homogenization of the Stokes equations, i.e. finding the limit system satisfied by the limit of $(\mathbf{u}_\varepsilon, p_\varepsilon)$ as $\varepsilon$ goes to zero, was first investigated by [13, 15, 1]. A review of these results can be found in [10], [5]. The homogenized or effective equations for the Stokes system, in a periodic perforated domain, is Darcy's law for the velocity $\mathbf{u}$ and the homogenized pressure $p^*$:

$$\left\{ \begin{array}{ll} \mathbf{u} = \frac{1}{\nu} A^*(\mathbf{f} - \nabla p^*) & \text{in } \Omega, \\ \operatorname{div} \mathbf{u} = 0 & \text{in } \Omega, \\ \mathbf{u} \cdot \mathbf{n} = 0 & \text{on } \partial\Omega, \end{array} \right.$$

where $A^*$ is a constant permeability tensor (see Proposition 3.2).

The first convergence result for the homogenization of the Stokes equations was established in [15], where the author proved the weak $L^2$ convergence of the velocity and the strong $L^2$ convergence of the pressure. The strong $L^2$ convergence of the velocity with a corrector was proven later in [1]. The first quantitative result of convergence in the $H^1$-norm was obtained in [12] with a relative error estimate of order $\varepsilon^{1/6}$ in a two-dimensional domain ($d=2$). This result has been improved to $\sqrt{\varepsilon}$ in [11] for the case of isolated solid obstacles (still in two space dimension $d=2$). This $\sqrt{\varepsilon}$ error estimate was independently and further improved in [14] for any space dimension $d \ge 2$, still under the assumption of isolated solid obstacles (as in Figure 1). The main technical idea in [14] is to construct boundary correctors which control the boundary layers appearing in the proof of the error estimate. Note that the assumption of isolated solid obstacles is not physically realistic in dimension $d \ge 3$.

The goal of the present paper is to extend this $\sqrt{\varepsilon}$ error estimate to the case of connected solid obstacles (as in Figure 2). Before stating our main result (in loose terms, see Theorem 4.1 for a precise statement) we need to introduce some notations. Actually the Darcy velocity $\mathbf{u}$ does not see the solid obstacles and thus cannot be a good approximation of the original velocity $\mathbf{u}_\varepsilon$. It must be corrected by introducing a better approximation which is

$$\mathbf{u}_2(x, \frac{x}{\varepsilon}) = \frac{1}{\nu} \sum_{i=1}^{d} \omega_i \left(\frac{x}{\varepsilon}\right) \left(f_i(x) - \frac{\partial p^*}{\partial x_i}(x)\right)$$

where the local velocities $\omega_i$ are solutions of cell problems (6). Actually the link between $\mathbf{u}_2$ and the Darcy velocity $\mathbf{u}$ is that $\mathbf{u}(x)$ is the average of $\mathbf{u}_2(x, y)$ with respect to the periodic variable $y$.

The main result we shall prove is:

**Theorem 1.1.** Let $\mathbf{u}_\varepsilon, p_\varepsilon$ be the solution to the Stokes problem (1) and $\mathbf{u}_2, p^*$, their homogenized approximations. Assuming that $\mathbf{f}$ is smooth enough, there exists a constant $C$, independent of $\varepsilon$ and $f$, such that

$$\frac{\|p_\varepsilon - p^*\|_{L^2(\Omega^\varepsilon)}}{\|p^*\|_{L^2(\Omega^\varepsilon)}} \le C\varepsilon^{\frac{1}{2}} \|\mathbf{f}\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})},$$
We introduce the bilinear forms $a: V \times V \to \mathbb{R}$ and $b: M \times V \to \mathbb{R}$ by:

$$a(\mathbf{u}, \mathbf{v}) = \nu(\nabla \mathbf{u}, \nabla \mathbf{v}), \quad b(p, \mathbf{v}) = (p, \operatorname{div} \mathbf{v}).$$

Then, assuming that $\mathbf{f} \in L^2(\Omega)^d$, a weak formulation of the Stokes problem (1) reads as follows: find $\mathbf{u}_\varepsilon \in V$ and $p_\varepsilon \in M$ such that

$$\left\{
\begin{array}{l}
a(\mathbf{u}_\varepsilon, \mathbf{v}) + b(p_\varepsilon, \mathbf{v}) = (\mathbf{f}, \mathbf{v}) \quad \forall \mathbf{v} \in V, \\
b(q, \mathbf{u}_\varepsilon) = 0 \quad \forall q \in M.
\end{array}
\right.
\qquad (2)$$

It is well known that, if $\Omega^\varepsilon$ is connected (to ensure that the zero-average condition of the pressure is enough to remove the undetermined constant of the pressure), there exists a unique weak solution to (2) [9].

The homogenization of the Stokes equations, i.e. finding the limit system satisfied by the limit of $(\mathbf{u}_\varepsilon, p_\varepsilon)$ as $\varepsilon$ goes to zero, was first investigated by [13, 15, 1]. A review of these results can be found in [10], [5]. The homogenized or effective equations for the Stokes system, in a periodic perforated domain, is Darcy's law for the velocity $\mathbf{u}$ and the homogenized pressure $p^*$:

$$\left\{
\begin{array}{ll}
\mathbf{u} = \frac{1}{\nu} A^*(\mathbf{f} - \nabla p^*) & \text{in } \Omega, \\
\operatorname{div} \mathbf{u} = 0 & \text{in } \Omega, \\
\mathbf{u} \cdot \mathbf{n} = 0 & \text{on } \partial\Omega,
\end{array}
\right.
$$

where $A^*$ is a constant permeability tensor (see Proposition 3.2).

The first convergence result for the homogenization of the Stokes equations was established in [15], where the author proved the weak $L^2$ convergence of the velocity and the strong $L^2$ convergence of the pressure. The strong $L^2$ convergence of the velocity with a corrector was proven later in [1]. The first quantitative result of convergence in the $H^1$-norm was obtained in [12] with a relative error estimate of order $\varepsilon^{1/6}$ in a two-dimensional domain ($d=2$). This result has been improved to $\sqrt{\varepsilon}$ in [11] for the case of isolated solid obstacles (still in two space dimension $d=2$). This $\sqrt{\varepsilon}$ error estimate was independently and further improved in [14] for any space dimension $d \ge 2$, still under the assumption of isolated solid obstacles (as in Figure 1). The main technical idea in [14] is to construct boundary correctors which control the boundary layers appearing in the proof of the error estimate. Note that the assumption of isolated solid obstacles is not physically realistic in dimension $d \ge 3$.

The goal of the present paper is to extend this $\sqrt{\varepsilon}$ error estimate to the case of connected solid obstacles (as in Figure 2). Before stating our main result (in loose terms, see Theorem 4.1 for a precise statement) we need to introduce some notations. Actually the Darcy velocity $\mathbf{u}$ does not see the solid obstacles and thus cannot be a good approximation of the original velocity $\mathbf{u}_\varepsilon$. It must be corrected by introducing a better approximation which is

$$\mathbf{u}_2(x, \frac{x}{\varepsilon}) = \frac{1}{\nu} \sum_{i=1}^{d} \omega_i \left(\frac{x}{\varepsilon}\right) \left(f_i(x) - \frac{\partial p^*}{\partial x_i}(x)\right)$$

where the local velocities $\omega_i$ are solutions of cell problems (6). Actually the link between $\mathbf{u}_2$ and the Darcy velocity $\mathbf{u}$ is that $\mathbf{u}(x)$ is the average of $\mathbf{u}_2(x, y)$ with respect to the periodic variable $y$.

The main result we shall prove is:

**Theorem 1.1.** Let $\mathbf{u}_\varepsilon, p_\varepsilon$ be the solution to the Stokes problem (1) and $\mathbf{u}_2, p^*$, their homogenized approximations. Assuming that $\mathbf{f}$ is smooth enough, there exists a constant $C$, independent of $\varepsilon$ and $f$, such that

$$\frac{\|p_\varepsilon - p^*\|_{L^2(\Omega^\varepsilon)}}{\|p^*\|_{L^2(\Omega^\varepsilon)}} \le C\varepsilon^{\frac{1}{2}} \|\mathbf{f}\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})},$$
$$ \frac{|u_{\varepsilon} - \varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})}}{|\varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})}} \leq C \varepsilon^{\frac{1}{2}} \|f\|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})}, $$

$$ \frac{\|u_{\varepsilon} - \varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})}}{\|\varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})}} \leq C \varepsilon^{\frac{1}{2}} \|f\|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})}. $$

Theorem 1.1 is stated in terms of relative errors since $\varepsilon^2 u_2$ (and thus $u_{\varepsilon}$) is small. Its proof follows the same strategy as that in [11] (but extends it to dimensions larger than 2) and differs significantly from the proof in [14]. Actually, our proof is simpler in the sense that no boundary layers (and thus correctors of their effect) appear. It turns out that we are able to use a cut-off argument as in the standard elliptic case [3] but without compromising the divergence-free condition for the velocities. Of course, the main novelty is that Theorem 1.1 applies to any case of solide obstacles, isolated or connected. In dimension $d \geq 3$, the physically relevant case is that of connected obstacles and not isolated obstacles.

The content of this paper is the following. Section 2 gives a precise definition of the geometrical assumptions of a periodic porous medium. Section 3 is devoted to the formal two-scale asymptotic expansion method which delivers the homogenized problem, as well as some cell correctors, including the velocities $\omega_i$ appearing in the definition of $u_2$. In Section 4 is stated our main result (Theorem 4.1 which is just a rigorous statement of Theorem 1.1) and its proof is decomposed in several technical lemmas.

# 2 Geometric modeling of the porous medium

As usual in periodic homogenization theory [3, 10, 13], we consider a porous medium obtained by the periodic repetition of an elementary cell of size $\varepsilon$, in a bounded domain of $\mathbb{R}^d$. We first define the corresponding dimensionless elementary cell $Y$.

## 2.1 Definition of the elementary cell $Y$

Let $Y=]0, 1]^d$ be the open unit cube of $\mathbb{R}^d$, $d \geq 2$. Let $Y_S$ be a closed subset of $\bar{Y}$. We define $Y_F$, an open set of $\mathbb{R}^d$, by $Y_F = Y \setminus Y_S$, where $Y_S$ represents the part of $Y$ occupied by the solid and $Y_F$ represents the part of $Y$ occupied by the fluid. The fluid part $Y_F$ and the solid part $Y_S$ satisfy $Y_F \cup Y_S = Y$ and $Y_F \cap Y_S = \emptyset$ (see Figure 3). The closed set $Y_S$ is repeated by $Y$-periodicity and fills the entire space $\mathbb{R}^d$, in order to obtain a closed set of $\mathbb{R}^d$, noted $E_S$. Let the open set $E_F$ be the complementary of $E_S$ in $\mathbb{R}^d$, i.e. $E_F = \mathbb{R}^d \setminus E_S$. We assume the following hypotheses on $Y_F$ and $E_F$:

1. $Y_F$ and $Y_S$ have strictly positive measure on $\bar{Y}$ (the elementary cell $Y$ contains fluid and solid together).

2. $Y_F$ is an open connected set with a locally Lipschitz boundary.

3. $E_F$ and the interior of $E_S$ are open sets with smooth boundaries of class $C^{m+2}$, with $m > \frac{d}{2}$, and are locally located on one side of their boundary. Moreover $E_F$ is connected.

Remark that the solid part $E_S$ can be connected or not, corresponding to the two different geometric cases of Figures 1 and 2 which are later called isolated obstacles or connected obstacles. The assumption that $E_F$ is smooth and connected implies that, on each face of $Y$, there is a fluid sub-domain of non-zero (surface) measure.
$$ \frac{|u_{\varepsilon} - \varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})}}{|\varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})}} \leq C \varepsilon^{\frac{1}{2}} \|f\|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})}, $$

$$ \frac{\|u_{\varepsilon} - \varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})}}{\|\varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})}} \leq C \varepsilon^{\frac{1}{2}} \|f\|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})}. $$

Theorem 1.1 is stated in terms of relative errors since $\varepsilon^2 u_2$ (and thus $u_{\varepsilon}$) is small. Its proof follows the same strategy as that in [11] (but extends it to dimensions larger than 2) and differs significantly from the proof in [14]. Actually, our proof is simpler in the sense that no boundary layers (and thus correctors of their effect) appear. It turns out that we are able to use a cut-off argument as in the standard elliptic case [3] but without compromising the divergence-free condition for the velocities. Of course, the main novelty is that Theorem 1.1 applies to any case of solide obstacles, isolated or connected. In dimension $d \geq 3$, the physically relevant case is that of connected obstacles and not isolated obstacles.

The content of this paper is the following. Section 2 gives a precise definition of the geometrical assumptions of a periodic porous medium. Section 3 is devoted to the formal two-scale asymptotic expansion method which delivers the homogenized problem, as well as some cell correctors, including the velocities $\omega_i$ appearing in the definition of $u_2$. In Section 4 is stated our main result (Theorem 4.1 which is just a rigorous statement of Theorem 1.1) and its proof is decomposed in several technical lemmas.

# 2 Geometric modeling of the porous medium

As usual in periodic homogenization theory [3, 10, 13], we consider a porous medium obtained by the periodic repetition of an elementary cell of size $\varepsilon$, in a bounded domain of $\mathbb{R}^d$. We first define the corresponding dimensionless elementary cell $Y$.

## 2.1 Definition of the elementary cell $Y$

Let $Y=]0, 1]^d$ be the open unit cube of $\mathbb{R}^d$, $d \geq 2$. Let $Y_S$ be a closed subset of $\bar{Y}$. We define $Y_F$, an open set of $\mathbb{R}^d$, by $Y_F = Y \setminus Y_S$, where $Y_S$ represents the part of $Y$ occupied by the solid and $Y_F$ represents the part of $Y$ occupied by the fluid. The fluid part $Y_F$ and the solid part $Y_S$ satisfy $Y_F \cup Y_S = Y$ and $Y_F \cap Y_S = \emptyset$ (see Figure 3). The closed set $Y_S$ is repeated by $Y$-periodicity and fills the entire space $\mathbb{R}^d$, in order to obtain a closed set of $\mathbb{R}^d$, noted $E_S$. Let the open set $E_F$ be the complementary of $E_S$ in $\mathbb{R}^d$, i.e. $E_F = \mathbb{R}^d \setminus E_S$. We assume the following hypotheses on $Y_F$ and $E_F$:

1. $Y_F$ and $Y_S$ have strictly positive measure on $\bar{Y}$ (the elementary cell $Y$ contains fluid and solid together).

2. $Y_F$ is an open connected set with a locally Lipschitz boundary.

3. $E_F$ and the interior of $E_S$ are open sets with smooth boundaries of class $C^{m+2}$, with $m > \frac{d}{2}$, and are locally located on one side of their boundary. Moreover $E_F$ is connected.

Remark that the solid part $E_S$ can be connected or not, corresponding to the two different geometric cases of Figures 1 and 2 which are later called isolated obstacles or connected obstacles. The assumption that $E_F$ is smooth and connected implies that, on each face of $Y$, there is a fluid sub-domain of non-zero (surface) measure.
<figure><img src="image_3.png" /><figcaption>Figure 3: Isolated (left) or connected (right) solid obstacles in the unit cell $Y = Y_F \cup Y_S$</figcaption></figure>

## 2.2 Definition of the open set $\Omega^\varepsilon$

Let $\Omega$ be a bounded and connected open set of $\mathbb{R}^d$ ($d \ge 2$) which is assumed to have a smooth boundary of class $C^{3,\alpha}$, for some $0 < \alpha < 1$. The set $\Omega$ is covered with a regular square mesh of size $\varepsilon$, each cell being a cube $Y_i^\varepsilon$. We define for any $i \in \mathbb{Z}^d$:

$$ 
\begin{gathered}
Y_i = Y + i, \quad Y_{i,F} = Y_F + i, \quad Y_{i,S} = Y_S + i, \\
Y_i^\varepsilon = \varepsilon(Y + i), \quad Y_{i,F}^\varepsilon = \varepsilon(Y_F + i), \quad Y_{i,S}^\varepsilon = \varepsilon(Y_S + i).
\end{gathered}
$$

The periodic set of solid obstacles (or perforations) inside $\Omega$ is defined as:

$$ B^\varepsilon = \bigcup_{i \in \mathcal{I}} Y_{i,S}^\varepsilon, $$

where $\mathcal{I} \subset \mathbb{Z}^d$ is a set of indices which is precisely defined in Assumption 2.1 below. The role of $\mathcal{I}$ is to remove obstacles from $\Omega$ only if they are sufficiently away from the boundary $\partial\Omega$. Finally, the fluid part $\Omega^\varepsilon$ of the perforated medium is defined by $\Omega^\varepsilon = \Omega \setminus B^\varepsilon$, i.e.

$$ \Omega^\varepsilon = \Omega \setminus \bigcup_{i \in \mathcal{I}} Y_{i,S}^\varepsilon. $$

The definition of $\mathcal{I}$ depends on the two cases of isolated or connected solid obstacles (see Figure 3).

**Assumption 2.1.** We consider two different sets of assumptions depending on the two cases of isolated or connected solid obstacles.

1. If the solid part $Y_S$ is isolated, namely strictly included in $Y$ (not touching its boundary), then $\mathcal{I} = \{i \in \mathbb{Z}^d \text{ s.t. } Y_i^\varepsilon \subset \Omega\}$. In other words, only entire obstacles $Y_{i,S}^\varepsilon$ are removed from $\Omega$ and thus no obstacles are cut by the boundary $\partial\Omega$.

2. If the solid part $Y_S$ touches the boundary of $Y$, meaning that obstacles are connected, then we define a first set of indices $\mathcal{I}_1 = \{i \in \mathbb{Z}^d \text{ s.t. } Y_i^\varepsilon \subset \Omega\}$ and an open subset of $\Omega$ defined by its closure

$$ \overline{C}_1^\varepsilon = \bigcup_{i \in \mathcal{I}_1} \overline{Y}_i^\varepsilon. $$

Then, a second smaller set of indices is

$$ \mathcal{I}_2 = \{i \in \mathcal{I}_1 \text{ such that } \overline{Y}_i^\varepsilon \cap \partial C_1^\varepsilon = \emptyset\}, $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Isolated (left) or connected (right) solid obstacles in the unit cell $Y = Y_F \cup Y_S$</figcaption></figure>

## 2.2 Definition of the open set $\Omega^\epsilon$

Let $\Omega$ be a bounded and connected open set of $\mathbb{R}^d$ ($d \ge 2$) which is assumed to have a smooth boundary of class $C^{3,\alpha}$, for some $0 < \alpha < 1$. The set $\Omega$ is covered with a regular square mesh of size $\epsilon$, each cell being a cube $Y_i^\epsilon$. We define for any $i \in \mathbb{Z}^d$:

$$ 
\begin{gathered}
Y_i = Y + i, \quad Y_{i,F} = Y_F + i, \quad Y_{i,S} = Y_S + i, \\
Y_i^\epsilon = \epsilon(Y + i), \quad Y_{i,F}^\epsilon = \epsilon(Y_F + i), \quad Y_{i,S}^\epsilon = \epsilon(Y_S + i).
\end{gathered}
$$

The periodic set of solid obstacles (or perforations) inside $\Omega$ is defined as:

$$ B^\epsilon = \bigcup_{i \in \mathcal{I}} Y_{i,S}^\epsilon, $$

where $\mathcal{I} \subset \mathbb{Z}^d$ is a set of indices which is precisely defined in Assumption 2.1 below. The role of $\mathcal{I}$ is to remove obstacles from $\Omega$ only if they are sufficiently away from the boundary $\partial\Omega$. Finally, the fluid part $\Omega^\epsilon$ of the perforated medium is defined by $\Omega^\epsilon = \Omega \setminus B^\epsilon$, i.e.

$$ \Omega^\epsilon = \Omega \setminus \bigcup_{i \in \mathcal{I}} Y_{i,S}^\epsilon. $$

The definition of $\mathcal{I}$ depends on the two cases of isolated or connected solid obstacles (see Figure 3).

**Assumption 2.1.** We consider two different sets of assumptions depending on the two cases of isolated or connected solid obstacles.

1. If the solid part $Y_S$ is isolated, namely strictly included in $Y$ (not touching its boundary), then $\mathcal{I} = \{i \in \mathbb{Z}^d \text{ s.t. } Y_i^\epsilon \subset \Omega\}$. In other words, only entire obstacles $Y_{i,S}^\epsilon$ are removed from $\Omega$ and thus no obstacles are cut by the boundary $\partial\Omega$.

2. If the solid part $Y_S$ touches the boundary of $Y$, meaning that obstacles are connected, then we define a first set of indices $\mathcal{I}_1 = \{i \in \mathbb{Z}^d \text{ s.t. } Y_i^\epsilon \subset \Omega\}$ and an open subset of $\Omega$ defined by its closure

$$ \overline{C}_1^\epsilon = \bigcup_{i \in \mathcal{I}_1} \overline{Y}_i^\epsilon. $$

Then, a second smaller set of indices is

$$ \mathcal{I}_2 = \{i \in \mathcal{I}_1 \text{ such that } \overline{Y}_i^\epsilon \cap \partial C_1^\epsilon = \emptyset\}, $$
meaning that no cell $Y_i^\epsilon$ in $\mathcal{I}_2$ touches the boundary $\partial C_1^\epsilon$. Similarly, we define an open subset of $\Omega$ defined by its closure

$$ \overline{C}_2^\epsilon = \bigcup_{i \in \mathcal{I}_2} \overline{Y}_i^\epsilon. $$

Finally, we choose $\mathcal{I} = \mathcal{I}_2$, meaning that the obstacles are not cut by the boundary $\partial\Omega$ and stay away from $\partial\Omega$ and $\partial C_1^\epsilon$ at a distance of the order of $\epsilon$ (see Figure 4, which is 2-d sketch of a 3-d situation).

**Remark 2.2.** The second item of Assumption 2.1 allows us to treat the realistic case of a connected solid part, which is possible only for space dimensions $d \ge 3$ because in 2-d $E_F$ and $E_S$ cannot be connected simultaneously. The second item of Assumption 2.1 is similar to an assumption in [5] and is slightly different from the assumption made in [1]. It is a necessary assumption for several technical reasons, including the fact that, when $E_S$ is connected, there may be some pathological cut obstacles near the boundary $\partial\Omega$. In particular, the fluid domain may be not connected.

<figure><img src="image_4.png" /><figcaption>Figure 4: Porous medium $\Omega^\epsilon$ when the solid part $Y_S$ touches the boundary of $Y$</figcaption></figure>

## 3 Two-scale asymptotic expansion

This section recalls how to formally obtain the homogenized problem for the Stokes equations (1), as well as the definitions of several corrector terms in the unit cell which are required for the sequel. This is achieved by applying the classical method of two-scale asymptotic expansion [3, 13]. The following results are already well-known but are recalled since the notations are required for the next sections. The formal method of two-scale asymptotic expansion starts with the following ansatz for the velocity $u_\epsilon$ and pressure $p_\epsilon$:

$$ \mathbf{u}_{\epsilon}(x) = \sum_{k=2}^{+\infty} \epsilon^k \mathbf{u}_k \left(x, \frac{x}{\epsilon}\right), \quad p_{\epsilon}(x) = \sum_{k=0}^{+\infty} \epsilon^k p_k \left(x, \frac{x}{\epsilon}\right). \qquad (3) $$

All functions $u_k(x, y)$, $p_k(x, y)$ are assumed $\mathbb{Z}^d$-periodic in $y$ i.e. 1-periodic with respect to each component $y_1, \dots, y_d$.

**Remark 3.1.** We begin the asymptotic expansion of the velocity $u_\epsilon$ at $k = 2$ since it is well-known in the literature that the first terms $u_0$ and $u_1$ vanish.
meaning that no cell $Y_i^\epsilon$ in $\mathcal{I}_2$ touches the boundary $\partial C_1^\epsilon$. Similarly, we define an open subset of $\Omega$ defined by its closure

$$ \overline{C}_2^\epsilon = \bigcup_{i \in \mathcal{I}_2} \overline{Y}_i^\epsilon. $$

Finally, we choose $\mathcal{I} = \mathcal{I}_2$, meaning that the obstacles are not cut by the boundary $\partial\Omega$ and stay away from $\partial\Omega$ and $\partial C_1^\epsilon$ at a distance of the order of $\epsilon$ (see Figure 4, which is 2-d sketch of a 3-d situation).

**Remark 2.2.** The second item of Assumption 2.1 allows us to treat the realistic case of a connected solid part, which is possible only for space dimensions $d \ge 3$ because in 2-d $E_F$ and $E_S$ cannot be connected simultaneously. The second item of Assumption 2.1 is similar to an assumption in [5] and is slightly different from the assumption made in [1]. It is a necessary assumption for several technical reasons, including the fact that, when $E_S$ is connected, there may be some pathological cut obstacles near the boundary $\partial\Omega$. In particular, the fluid domain may be not connected.

<figure><img src="image_6.png" /><figcaption>Figure 4: Porous medium $\Omega^\epsilon$ when the solid part $Y_S$ touches the boundary of $Y$</figcaption></figure>

## 3 Two-scale asymptotic expansion

This section recalls how to formally obtain the homogenized problem for the Stokes equations (1), as well as the definitions of several corrector terms in the unit cell which are required for the sequel. This is achieved by applying the classical method of two-scale asymptotic expansion [3, 13]. The following results are already well-known but are recalled since the notations are required for the next sections. The formal method of two-scale asymptotic expansion starts with the following ansatz for the velocity $u_\epsilon$ and pressure $p_\epsilon$:

$$ \mathbf{u}_{\epsilon}(x) = \sum_{k=2}^{+\infty} \epsilon^k \mathbf{u}_k \left(x, \frac{x}{\epsilon}\right), \quad p_{\epsilon}(x) = \sum_{k=0}^{+\infty} \epsilon^k p_k \left(x, \frac{x}{\epsilon}\right). \qquad (3) $$

All functions $u_k(x, y)$, $p_k(x, y)$ are assumed $\mathbb{Z}^d$-periodic in $y$ i.e. 1-periodic with respect to each component $y_1, \dots, y_d$.

**Remark 3.1.** We begin the asymptotic expansion of the velocity $u_\epsilon$ at $k = 2$ since it is well-known in the literature that the first terms $u_0$ and $u_1$ vanish.
We recall the following derivation rule:

$$ \nabla \left( \mathbf{u}_k \left( x, \frac{x}{\epsilon} \right) \right) = \left( \epsilon^{-1} \nabla_y \mathbf{u}_k + \nabla \mathbf{u}_k \right) \left( x, \frac{x}{\epsilon} \right), \quad (4) $$

where $\nabla$ and $\nabla_y$ denote the partial derivatives with respect to the macroscopic variable $x$ and the microscopic variable $y$. Then, introducing the series (3) in the Stokes equations (1), using the derivation rule (4), we obtain for the momentum equation:

$$ \begin{aligned} -\nu \Delta \mathbf{u}_{\varepsilon}(x) + \nabla p_{\varepsilon}(x) ={}& \epsilon^{-1} \nabla_y p_0 + \epsilon^0 \{-\nu \operatorname{div}_y \nabla_y \mathbf{u}_2 + \nabla p_0 + \nabla_y p_1\} \\ & + \epsilon \{-\nu (\operatorname{div} \nabla_y \mathbf{u}_2 + \operatorname{div}_y \nabla_y \mathbf{u}_3) + \nabla p_1 + \nabla_y p_2\} \\ & + \sum_{k=2}^{+\infty} \epsilon^k \{-\nu (\operatorname{div} \nabla \mathbf{u}_k + \operatorname{div} \nabla_y \mathbf{u}_{k+1} + \operatorname{div}_y \nabla \mathbf{u}_{k+1} + \operatorname{div}_y \nabla_y \mathbf{u}_{k+2}) \\ & \qquad + \nabla p_k + \nabla_y p_{k+1}\} = \mathbf{f} \end{aligned} $$

and for the incompressibility condition:

$$ \operatorname{div} \mathbf{u}_{\varepsilon}(x) = \epsilon \operatorname{div}_y \mathbf{u}_2 + \sum_{k=2}^{+\infty} \epsilon^k \{\operatorname{div} \mathbf{u}_k + \operatorname{div}_y \mathbf{u}_{k+1}\} = 0. $$

Now, we identify the quantities associated to the different orders of $\epsilon$. The leading order $\epsilon^{-1}$ of the momentum equation is:

$$ \left\{ \begin{array}{ll} \nabla_y p_0 = 0 & \text{in } Y_F, \\ y \to p_0(x, y) & \text{Y-periodic.} \end{array} \right. \qquad (5) $$

From (5) one deduces that $p_0$ does not depend on $y$ and there exists a function $p^*(x)$ such that $p_0(x, y) = p^*(x)$. At the next order $\epsilon^0$ for the momentum equation and $\epsilon$ for the incompressibility condition, we get:

$$ \left\{ \begin{array}{ll} -\nu \Delta_y \mathbf{u}_2 + \nabla_y p_1 = \mathbf{f} - \nabla p^* & \text{in } Y_F, \\ \operatorname{div}_y \mathbf{u}_2 = 0 & \text{in } Y_F, \\ \mathbf{u}_2 = \mathbf{0} & \text{on } \partial Y_S, \\ y \to \mathbf{u}_2(x, y), p_1(x, y) & \text{Y-periodic.} \end{array} \right. $$

It follows that $\mathbf{u}_2$ and $p_1$ are linear combinations of the solutions to the following cell Stokes problems: for $i = 1, \dots, d$ find $\omega_i : Y \to \mathbb{R}^d$ and $\pi_i : Y \to \mathbb{R}$, $\mathbb{Z}^d$-periodic and solutions to:

$$ \left\{ \begin{array}{ll} -\Delta_y \omega_i + \nabla_y \pi_i = e_i & \text{in } Y_F, \\ \operatorname{div}_y \omega_i = 0 & \text{in } Y_F, \\ \omega_i = \mathbf{0} & \text{on } \partial Y_S, \\ \int_Y \pi_i = 0, \\ y \to \omega_i, \pi_i & \text{Y-periodic.} \end{array} \right. \qquad (6) $$

Velocity and pressure are extended by 0 in $Y_S$. Using the Einstein summation convention on repeated indices, it follows that:

$$ \left\{ \begin{array}{l} \mathbf{u}_2(x, y) = \frac{1}{\nu} \omega_i(y) \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right), \\ p_1(x, y) = \pi_i(y) \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right). \end{array} \right. \qquad (7) $$
We recall the following derivation rule:

$$ \nabla \left( \mathbf{u}_k \left( x, \frac{x}{\epsilon} \right) \right) = \left( \epsilon^{-1} \nabla_y \mathbf{u}_k + \nabla \mathbf{u}_k \right) \left( x, \frac{x}{\epsilon} \right), \quad (4) $$

where $\nabla$ and $\nabla_y$ denote the partial derivatives with respect to the macroscopic variable $x$ and the microscopic variable $y$. Then, introducing the series (3) in the Stokes equations (1), using the derivation rule (4), we obtain for the momentum equation:

$$ \begin{aligned} -\nu \Delta \mathbf{u}_{\varepsilon}(x) + \nabla p_{\varepsilon}(x) ={}& \epsilon^{-1} \nabla_y p_0 + \epsilon^0 \left\{ -\nu \operatorname{div}_y \nabla_y \mathbf{u}_2 + \nabla p_0 + \nabla_y p_1 \right\} \\ & + \epsilon \left\{ -\nu \left( \operatorname{div} \nabla_y \mathbf{u}_2 + \operatorname{div}_y \nabla_y \mathbf{u}_3 \right) + \nabla p_1 + \nabla_y p_2 \right\} \\ & + \sum_{k=2}^{+\infty} \epsilon^k \left\{ -\nu \left( \operatorname{div} \nabla \mathbf{u}_k + \operatorname{div} \nabla_y \mathbf{u}_{k+1} + \operatorname{div}_y \nabla \mathbf{u}_{k+1} + \operatorname{div}_y \nabla_y \mathbf{u}_{k+2} \right) \right\} \\ & + \nabla p_k + \nabla_y p_{k+1} \} = \mathbf{f} \end{aligned} $$

and for the incompressibility condition:

$$ \operatorname{div} \mathbf{u}_{\varepsilon}(x) = \epsilon \operatorname{div}_y \mathbf{u}_2 + \sum_{k=2}^{+\infty} \epsilon^k \{ \operatorname{div} \mathbf{u}_k + \operatorname{div}_y \mathbf{u}_{k+1} \} = 0. $$

Now, we identify the quantities associated to the different orders of $\epsilon$. The leading order $\epsilon^{-1}$ of the momentum equation is:

$$ \left\{ \begin{array}{ll} \nabla_y p_0 = 0 & \text{in } Y_F, \\ y \to p_0(x, y) & \text{Y-periodic.} \end{array} \right. \qquad (5) $$

From (5) one deduces that $p_0$ does not depend on $y$ and there exists a function $p^*(x)$ such that $p_0(x, y) = p^*(x)$. At the next order $\epsilon^0$ for the momentum equation and $\epsilon$ for the incompressibility condition, we get:

$$ \left\{ \begin{array}{ll} -\nu \Delta_y \mathbf{u}_2 + \nabla_y p_1 = \mathbf{f} - \nabla p^* & \text{in } Y_F, \\ \operatorname{div}_y \mathbf{u}_2 = 0 & \text{in } Y_F, \\ \mathbf{u}_2 = \mathbf{0} & \text{on } \partial Y_S, \\ y \to \mathbf{u}_2(x, y), \quad p_1(x, y) & \text{Y-periodic.} \end{array} \right. $$

It follows that $\mathbf{u}_2$ and $p_1$ are linear combinations of the solutions to the following cell Stokes problems: for $i = 1, \dots, d$ find $\omega_i : Y \to \mathbb{R}^d$ and $\pi_i : Y \to \mathbb{R}$, $\mathbb{Z}^d$-periodic and solutions to:

$$ \left\{ \begin{array}{ll} -\Delta_y \omega_i + \nabla_y \pi_i = e_i & \text{in } Y_F, \\ \operatorname{div}_y \omega_i = 0 & \text{in } Y_F, \\ \omega_i = \mathbf{0} & \text{on } \partial Y_S, \\ \int_Y \pi_i = 0, \\ y \to \omega_i, \quad \pi_i & \text{Y-periodic.} \end{array} \right. \qquad (6) $$

Velocity and pressure are extended by 0 in $Y_S$. Using the Einstein summation convention on repeated indices, it follows that:

$$ \left\{ \begin{array}{l} \mathbf{u}_2(x, y) = \frac{1}{\nu} \omega_i(y) \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right), \\ p_1(x, y) = \pi_i(y) \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right). \end{array} \right. \qquad (7) $$
Furthermore, averaging on $Y_F$ the $\epsilon^2$ term of the incompressibility condition (taking into account the boundary conditions for $u_3$ in (11)) yields

$$ \mathrm{div} \langle u_2 \rangle = 0, \qquad (8) $$

where $\langle \cdot \rangle$ stands for the average over $Y_F$:

$$ \langle v \rangle = \frac{1}{|Y_F|} \int_{Y_F} v(y) dy. $$

Inserting formula (7) for $u_2$ in (8) leads to the Darcy equation for the pressure $p^*$. Furthermore, a similar asymptotic analysis can be performed on the boundary $\partial\Omega$ but the homogenized boundary condition turns out to be weaker, namely only the normal component of $u$ vanishes on $\partial\Omega$ [1, 10, 11].

**Proposition 3.2.** *The Darcy velocity is defined as*

$$ u = \langle u_2 \rangle = \frac{1}{\nu} \langle \omega_i \rangle \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) \qquad (9) $$

and the homogenized pressure $p^*$ is the solution to the following Darcy problem:

$$ \left\{ \begin{array}{ll} \mathrm{div} A^*(f - \nabla p^*) = 0 & \text{in } \Omega, \\ n \cdot (A^*(f - \nabla p^*)) = 0 & \text{on } \partial\Omega, \end{array} \right. \qquad (10) $$

where the constant tensor $A^*$ is defined by its columns

$$ A^* e_i = \langle \omega_i \rangle = \frac{1}{|Y_F|} \int_{Y_F} \omega_i(y) dy. $$

There exists a solution $p^* \in H^1(\Omega)$ of (10), which is unique up to an additive constant. Furthermore, under Assumption 2.1 and assuming $f \in C^{1,\alpha}(\bar{\Omega})^d$ with $0 < \alpha < 1$, then $\nabla p^*$ is of class $C^{1,\alpha}(\bar{\Omega})$ and $p^*$ is of class $C^{2,\alpha}(\bar{\Omega})$. Similarly, if $f \in H^2(\Omega)^d$, then $\nabla p^* \in H^2(\Omega)^d$ too.

*Proof.* Problem (10) is well-posed in $H^1(\Omega)/\mathbb{R}$ because it is just a second-order elliptic equation for the pressure $p^*$, complemented by a Neumann boundary condition (one can check that $A^*$ is positive definite). Recall that $C^{1,\alpha}(\bar{\Omega})$ is the space of functions in $C^1(\bar{\Omega})$ with derivatives which are $\alpha$-Hölderian. Since $f \in C^{1,\alpha}(\bar{\Omega})$ then it is known that $\nabla p^*$ is of class $C^{1,\alpha}(\bar{\Omega})$ (see for example [8, Theorems 8.33 and 8.34]). Finally, $f \in H^2(\Omega)^d$ implies $\nabla p^* \in H^2(\Omega)^d$ by standard Sobolev regularity theory for elliptic equations. $\square$

At the next order $\epsilon$ for the momentum equation and $\epsilon^2$ for the incompressibility condition, we get:

$$ \left\{ \begin{array}{ll} -\nu \Delta_y u_3 + \nabla_y p_2 = \nu \mathrm{div} \, \nabla_y u_2 + \nu \mathrm{div}_y \nabla u_2 - \nabla p_1 & \text{in } Y_F, \\ \mathrm{div}_y u_3 = -\mathrm{div} \, u_2 & \text{in } Y_F, \\ u_3 = \mathbf{0} & \text{on } \partial Y_S, \\ y \to u_3(x, y), \quad p_2(x, y) & \text{Y-periodic.} \end{array} \right. \qquad (11) $$

Substituting the expressions (7) of $u_2$ and $p_1$ in (11) leads to

$$ -\nu \Delta_y u_3 + \nabla_y p_2 = (2 \nabla_y \omega_i(y) - \pi_i(y) I_d) \nabla \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right). $$

By linearity, we deduce that

$$ \left\{ \begin{array}{l} u_3(x, y) = \frac{1}{\nu} \gamma_{ij}(y) \frac{\partial}{\partial x_j} \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right), \\ p_2(x, y) = \theta_{ij}(y) \frac{\partial}{\partial x_j} \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right), \end{array} \right. $$
Furthermore, averaging on $Y_F$ the $\varepsilon^2$ term of the incompressibility condition (taking into account the boundary conditions for $u_3$ in (11)) yields

$$ \mathrm{div} \langle u_2 \rangle = 0, \qquad (8) $$

where $\langle \cdot \rangle$ stands for the average over $Y_F$:

$$ \langle v \rangle = \frac{1}{|Y_F|} \int_{Y_F} v(y) dy. $$

Inserting formula (7) for $u_2$ in (8) leads to the Darcy equation for the pressure $p^*$. Furthermore, a similar asymptotic analysis can be performed on the boundary $\partial\Omega$ but the homogenized boundary condition turns out to be weaker, namely only the normal component of $u$ vanishes on $\partial\Omega$ [1, 10, 11].

**Proposition 3.2.** The Darcy velocity is defined as

$$ u = \langle u_2 \rangle = \frac{1}{\nu} \langle \omega_i \rangle \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) \qquad (9) $$

and the homogenized pressure $p^*$ is the solution to the following Darcy problem:

$$ \left\{ \begin{array}{ll} \mathrm{div} A^*(f - \nabla p^*) = 0 & \text{in } \Omega, \\ n \cdot (A^*(f - \nabla p^*)) = 0 & \text{on } \partial\Omega, \end{array} \right. \qquad (10) $$

where the constant tensor $A^*$ is defined by its columns

$$ A^* e_i = \langle \omega_i \rangle = \frac{1}{|Y_F|} \int_{Y_F} \omega_i(y) dy. $$

There exists a solution $p^* \in H^1(\Omega)$ of (10), which is unique up to an additive constant. Furthermore, under Assumption 2.1 and assuming $f \in C^{1,\alpha}(\bar{\Omega})^d$ with $0 < \alpha < 1$, then $\nabla p^*$ is of class $C^{1,\alpha}(\bar{\Omega})$ and $p^*$ is of class $C^{2,\alpha}(\bar{\Omega})$. Similarly, if $f \in H^2(\Omega)^d$, then $\nabla p^* \in H^2(\Omega)^d$ too.

**Proof.** Problem (10) is well-posed in $H^1(\Omega)/\mathbb{R}$ because it is just a second-order elliptic equation for the pressure $p^*$, complemented by a Neumann boundary condition (one can check that $A^*$ is positive definite). Recall that $C^{1,\alpha}(\bar{\Omega})$ is the space of functions in $C^1(\bar{\Omega})$ with derivatives which are $\alpha$-Hölderian. Since $f \in C^{1,\alpha}(\bar{\Omega})$ then it is known that $\nabla p^*$ is of class $C^{1,\alpha}(\bar{\Omega})$ (see for example [8, Theorems 8.33 and 8.34]). Finally, $f \in H^2(\Omega)^d$ implies $\nabla p^* \in H^2(\Omega)^d$ by standard Sobolev regularity theory for elliptic equations. $\square$

At the next order $\varepsilon$ for the momentum equation and $\varepsilon^2$ for the incompressibility condition, we get:

$$ \left\{ \begin{array}{ll} -\nu \Delta_y u_3 + \nabla_y p_2 = \nu \mathrm{div} \, \nabla_y u_2 + \nu \mathrm{div}_y \nabla u_2 - \nabla p_1 & \text{in } Y_F, \\ \mathrm{div}_y u_3 = -\mathrm{div} \, u_2 & \text{in } Y_F, \\ u_3 = 0 & \text{on } \partial Y_S, \\ y \to u_3(x, y), \ p_2(x, y) & \text{Y-periodic.} \end{array} \right. \qquad (11) $$

Substituting the expressions (7) of $u_2$ and $p_1$ in (11) leads to

$$ -\nu \Delta_y u_3 + \nabla_y p_2 = (2 \nabla_y \omega_i(y) - \pi_i(y) I_d) \nabla \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right). $$

By linearity, we deduce that

$$ \left\{ \begin{array}{l} u_3(x, y) = \frac{1}{\nu} \gamma_{ij}(y) \frac{\partial}{\partial x_j} \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right), \\ p_2(x, y) = \theta_{ij}(y) \frac{\partial}{\partial x_j} \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right), \end{array} \right. $$
where, for $i, j = 1, \dots, d$, $\gamma_{ij} : Y \to \mathbb{R}^d$ and $\theta_{ij} : Y \to \mathbb{R}$ are $\mathbb{Z}^d$-periodic solutions to another cell Stokes problem:

$$ egin{cases} 
-\Delta_y \gamma_{ij} + \nabla_y \theta_{ij} = 2 \frac{\partial \omega_i}{\partial y_j} - \pi_i e_j & \text{in } Y_F, \\
\text{div}_y \gamma_{ij} = -\omega_i \cdot e_j + \langle \omega_i \cdot e_j \rangle & \text{in } Y_F, \\
\gamma_{ij} = \mathbf{0} & \text{on } \partial Y_S, \\
\int_Y \theta_{ij} = 0, \\
y \to \gamma_{ij}, \theta_{ij} & \text{Y-periodic.}
\end{cases} 
\quad (12) $$

We recall some notations on Sobolev spaces of periodic functions (see [10, Appendix B.2.2] for details). The $L^2$ space of Y-periodic functions is

$$ L_{\#}^2(Y) = \{ \varphi \in L_{\text{loc}}^2(\mathbb{R}^d), \varphi \text{ is Y-periodic} \}, $$

endowed with the usual scalar product and norm of $L^2(Y)$. It turns out that $L_{\#}^2(Y)$ coincides with $L^2(Y)$. Higher-order Sobolev spaces are defined, for $m \ge 1$, by

$$ H_{\#}^m(Y) = \{ \varphi \in H_{\text{loc}}^m(\mathbb{R}^d), \varphi \text{ is Y-periodic} \}, $$

endowed with the usual scalar product and norm of $H^m(Y)$. Let us recall a classical regularity result [4].

**Lemma 3.3.** Let $m \ge 0$ and $Y_F$ be of class $C^{m+2}$, consider $f(y) \in H_{\#}^m(Y)^d$ and $g(y) \in H_{\#}^{m+1}(Y)$ with zero mean on Y. The problem

$$ egin{cases} 
-\Delta_y u + \nabla_y p = f(y) & \text{in } Y_F, \\
\text{div}_y u = g(y) & \text{in } Y_F, \\
u = \mathbf{0} & \text{on } \partial Y_S, \\
y \to u(y), p(y) & \text{Y-periodic,}
\end{cases} $$

admits a unique solution $u \in H_{0,\#}^{m+2}(Y_F)^d$ and $p \in H_{\#}^{m+1}(Y_F)/\mathbb{R}$, where the space $H_{0,\#}^k(Y_F)$ is defined by

$$ H_{0,\#}^k(Y_F) = \{ \varphi \in H_{\#}^k(Y_F) \text{ such that } \varphi = \mathbf{0} \text{ on } \partial Y_S \} $$

Using Lemma 3.3, it is clear that (6) and (12) admit each a unique solution since the right-hand side of the divergence condition has zero-average in $Y_F$. In particular $\omega_i$ and $\gamma_{i,j}$ are in $H_{0,\#}^{m+2}(Y_F)^d$, and, $\pi_i$ and $\theta_{i,j}$ are in $H_{\#}^{m+1}(Y_F)/\mathbb{R}$.

The conclusion of this section is that we obtained an approximation of order 3 of the velocity

$$ u_{\varepsilon}(x) \approx \varepsilon^2 u_2 \left(x, \frac{x}{\varepsilon}\right) + \varepsilon^3 u_3 \left(x, \frac{x}{\varepsilon}\right), \quad (13) $$

and of order 1 of the pressure

$$ p_{\varepsilon}(x) \approx p^{*}(x) + \varepsilon p_{1} \left(x, \frac{x}{\varepsilon}\right). \quad (14) $$

The derivation of (13) and (14) is merely formal but these approximations are key ingredients in the rigorous proof of our main result of quantitative homogenization.
where, for $i, j = 1, \dots, d$, $\gamma_{ij} : Y \to \mathbb{R}^d$ and $\theta_{ij} : Y \to \mathbb{R}$ are $\mathbb{Z}^d$-periodic solutions to another cell Stokes problem:

$$ egin{cases} 
-\Delta_y \gamma_{ij} + \nabla_y \theta_{ij} = 2 \frac{\partial \omega_i}{\partial y_j} - \pi_i e_j & \text{in } Y_F, \\
\text{div}_y \gamma_{ij} = -\omega_i \cdot e_j + \langle \omega_i \cdot e_j \rangle & \text{in } Y_F, \\
\gamma_{ij} = \mathbf{0} & \text{on } \partial Y_S, \\
\int_Y \theta_{ij} = 0, \\
y \to \gamma_{ij}, \theta_{ij} & \text{Y-periodic.} 
\end{cases} 
\quad (12) $$

We recall some notations on Sobolev spaces of periodic functions (see [10, Appendix B.2.2] for details). The $L^2$ space of Y-periodic functions is

$$ L_{\#}^2(Y) = \{\varphi \in L_{\text{loc}}^2(\mathbb{R}^d), \varphi \text{ is Y-periodic}\}, $$

endowed with the usual scalar product and norm of $L^2(Y)$. It turns out that $L_{\#}^2(Y)$ coincides with $L^2(Y)$. Higher-order Sobolev spaces are defined, for $m \ge 1$, by

$$ H_{\#}^m(Y) = \{\varphi \in H_{\text{loc}}^m(\mathbb{R}^d), \varphi \text{ is Y-periodic}\}, $$

endowed with the usual scalar product and norm of $H^m(Y)$. Let us recall a classical regularity result [4].

**Lemma 3.3.** Let $m \ge 0$ and $Y_F$ be of class $C^{m+2}$, consider $f(y) \in H_{\#}^m(Y)^d$ and $g(y) \in H_{\#}^{m+1}(Y)$ with zero mean on Y. The problem

$$ egin{cases} 
-\Delta_y u + \nabla_y p = f(y) & \text{in } Y_F, \\
\text{div}_y u = g(y) & \text{in } Y_F, \\
u = \mathbf{0} & \text{on } \partial Y_S, \\
y \to u(y), p(y) & \text{Y-periodic,} 
\end{cases} $$

admits a unique solution $u \in H_{0,\#}^{m+2}(Y_F)^d$ and $p \in H_{\#}^{m+1}(Y_F)/\mathbb{R}$, where the space $H_{0,\#}^k(Y_F)$ is defined by

$$ H_{0,\#}^k(Y_F) = \{\varphi \in H_{\#}^k(Y_F) \text{ such that } \varphi = \mathbf{0} \text{ on } \partial Y_S\} $$

Using Lemma 3.3, it is clear that (6) and (12) admit each a unique solution since the right-hand side of the divergence condition has zero-average in $Y_F$. In particular $\omega_i$ and $\gamma_{i,j}$ are in $H_{0,\#}^{m+2}(Y_F)^d$, and, $\pi_i$ and $\theta_{i,j}$ are in $H_{\#}^{m+1}(Y_F)/\mathbb{R}$.

The conclusion of this section is that we obtained an approximation of order 3 of the velocity

$$ u_{\varepsilon}(x) \approx \varepsilon^2 u_2 \left(x, \frac{x}{\varepsilon}\right) + \varepsilon^3 u_3 \left(x, \frac{x}{\varepsilon}\right), \quad (13) $$

and of order 1 of the pressure

$$ p_{\varepsilon}(x) \approx p^{*}(x) + \varepsilon p_{1} \left(x, \frac{x}{\varepsilon}\right). \quad (14) $$

The derivation of (13) and (14) is merely formal but these approximations are key ingredients in the rigorous proof of our main result of quantitative homogenization.
# 4 A quantitative homogenization result

## 4.1 Main result

Our main result is a quantitative error estimate for the homogenization of Stokes equations in a periodic porous medium, which was loosely announced as Theorem 1.1 in the introduction. A more precise statement of this result is given below.

**Theorem 4.1.** Let $u_{\varepsilon}$, $p_{\varepsilon}$ be the solution to the Stokes equations (1), $p^*$ be the solution to the Darcy equation (10) and $u_2$ be defined by (7). Under Assumption 2.1 and assuming that $f$ belongs to $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$ with $0 < \alpha < 1$, there holds:

$$ \|p_{\varepsilon}(x) - p^{*}(x)\|_{L^2(\Omega^{\varepsilon})} \leq C\varepsilon^{\frac{1}{2}}\|f(x) - \nabla p^{*}(x)\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (15) $$

$$ |u_{\varepsilon}(x) - \varepsilon^2 u_2\left(x, \frac{x}{\varepsilon}\right)|_{H^1(\Omega^{\varepsilon})} \leq C\varepsilon^{\frac{3}{2}}\|f(x) - \nabla p^{*}(x)\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (16) $$

$$ \|u_{\varepsilon}(x) - \varepsilon^2 u_2\left(x, \frac{x}{\varepsilon}\right)\|_{L^2(\Omega^{\varepsilon})} \leq C\varepsilon^{\frac{5}{2}}\|f(x) - \nabla p^{*}(x)\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (17) $$

where $C$ is independent of $\varepsilon$ and the data.

Theorem 4.1 was already proved in the easier case of isolated obstacles in [11] when $d=2$ and in [14] for any dimension. The main novelty of the present paper is that Theorem 4.1 applies to the physically relevant case of connected obstacles (in dimension $d \geq 3$).

**Remark 4.2.** In Theorem 4.1, the regularity of $f$ ensures that $\nabla p^*$ also belongs to $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$ (see Proposition 3.2). We assume that $f \in H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$, but all what follows only requires the assumption that the difference $(f - \nabla p^*)$ belongs to $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$. The space $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$ is equipped with the norm $\|\cdot\|_{H^2 \cap C^{1,\alpha}}$ defined by $\|\cdot\|_{H^2 \cap C^{1,\alpha}} = \|\cdot\|_{H^2} + \|\cdot\|_{C^{1,\alpha}}$.

**Remark 4.3.** The two estimates for the velocity essentially mean that the relative error for the velocity in the $H^1$-norm is of the order $\sqrt{\varepsilon}$. Indeed, the velocity is of order $\varepsilon^2$ (be it $u_{\varepsilon}$ or $\varepsilon^2 u_2$), but its gradient is of order $\varepsilon$ (see the derivative rule (4)) since they oscillate on the length scale $\varepsilon$. Finally the relative error for the pressure is also of order $\sqrt{\varepsilon}$. Eventually, upon rescaling the velocity, we can write:

$$ \varepsilon^{-2}\|u_{\varepsilon} - \varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})} + \varepsilon^{-1}|u_{\varepsilon} - \varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})} + \|p_{\varepsilon} - p^*\|_{L^2(\Omega^{\varepsilon})} \leq C\sqrt{\varepsilon}\|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} $$

The relative error of order $\sqrt{\varepsilon}$ is worse than the expected $\varepsilon$. It is a classical phenomenon in homogenization [3, 2, 11] due to boundary layer effects caused by the Dirichlet boundary condition on $\partial\Omega$. Technically, this is taken into account by introducing a suitably chosen cut-off function $\eta^{\varepsilon}$ in the proof. If the boundary layers were absent, which would be the case, for example, with periodic boundary conditions over a rectangular box $\Omega = \prod_{i=1}^d (0, \varepsilon n_i)$ with $n_i \in \mathbb{N}$, the error estimate would improve and give a relative error of order $\varepsilon$:

$$ \varepsilon^{-2}\|u_{\varepsilon} - \varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})} + \varepsilon^{-1}|u_{\varepsilon} - \varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})} + \|p_{\varepsilon} - p^*\|_{L^2(\Omega^{\varepsilon})} \leq C\varepsilon\|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} $$

The proof of Theorem 4.1 is given in Subsection 4.5. Beforehand, it requires three technical lemmas concerning (i) solving a non-homogeneous divergence equation in a porous domain (Lemma 4.9), (ii) correcting non-homogeneous boundary conditions (Lemma 4.13) and (iii) representing periodic divergence-free fields (Lemma 4.15).
# 4 A quantitative homogenization result

## 4.1 Main result

Our main result is a quantitative error estimate for the homogenization of Stokes equations in a periodic porous medium, which was loosely announced as Theorem 1.1 in the introduction. A more precise statement of this result is given below.

**Theorem 4.1.** Let $u_{\varepsilon}$, $p_{\varepsilon}$ be the solution to the Stokes equations (1), $p^*$ be the solution to the Darcy equation (10) and $u_2$ be defined by (7). Under Assumption 2.1 and assuming that $f$ belongs to $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$ with $0 < \alpha < 1$, there holds:

$$ \|p_{\varepsilon}(x) - p^{*}(x)\|_{L^2(\Omega^{\varepsilon})} \leq C\varepsilon^{\frac{1}{2}}\|f(x) - \nabla p^{*}(x)\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (15) $$

$$ |u_{\varepsilon}(x) - \varepsilon^2 u_2\left(x, \frac{x}{\varepsilon}\right)|_{H^1(\Omega^{\varepsilon})} \leq C\varepsilon^{\frac{3}{2}}\|f(x) - \nabla p^{*}(x)\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (16) $$

$$ \|u_{\varepsilon}(x) - \varepsilon^2 u_2\left(x, \frac{x}{\varepsilon}\right)\|_{L^2(\Omega^{\varepsilon})} \leq C\varepsilon^{\frac{5}{2}}\|f(x) - \nabla p^{*}(x)\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (17) $$

where $C$ is independent of $\varepsilon$ and the data.

Theorem 4.1 was already proved in the easier case of isolated obstacles in [11] when $d=2$ and in [14] for any dimension. The main novelty of the present paper is that Theorem 4.1 applies to the physically relevant case of connected obstacles (in dimension $d \geq 3$).

**Remark 4.2.** In Theorem 4.1, the regularity of $f$ ensures that $\nabla p^*$ also belongs to $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$ (see Proposition 3.2). We assume that $f \in H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$, but all what follows only requires the assumption that the difference $(f - \nabla p^*)$ belongs to $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$. The space $H^2(\Omega)^d \cap C^{1,\alpha}(\bar{\Omega})^d$ is equipped with the norm $\|\cdot\|_{H^2 \cap C^{1,\alpha}}$ defined by $\|\cdot\|_{H^2 \cap C^{1,\alpha}} = \|\cdot\|_{H^2} + \|\cdot\|_{C^{1,\alpha}}$.

**Remark 4.3.** The two estimates for the velocity essentially mean that the relative error for the velocity in the $H^1$-norm is of the order $\sqrt{\varepsilon}$. Indeed, the velocity is of order $\varepsilon^2$ (be it $u_{\varepsilon}$ or $\varepsilon^2 u_2$), but its gradient is of order $\varepsilon$ (see the derivative rule (4)) since they oscillate on the length scale $\varepsilon$. Finally the relative error for the pressure is also of order $\sqrt{\varepsilon}$. Eventually, upon rescaling the velocity, we can write:

$$ \varepsilon^{-2}\|u_{\varepsilon} - \varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})} + \varepsilon^{-1}|u_{\varepsilon} - \varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})} + \|p_{\varepsilon} - p^*\|_{L^2(\Omega^{\varepsilon})} \leq C\sqrt{\varepsilon}\|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} $$

The relative error of order $\sqrt{\varepsilon}$ is worse than the expected $\varepsilon$. It is a classical phenomenon in homogenization [3, 2, 11] due to boundary layer effects caused by the Dirichlet boundary condition on $\partial\Omega$. Technically, this is taken into account by introducing a suitably chosen cut-off function $\eta^{\varepsilon}$ in the proof. If the boundary layers were absent, which would be the case, for example, with periodic boundary conditions over a rectangular box $\Omega = \prod_{i=1}^d (0, \varepsilon n_i)$ with $n_i \in \mathbb{N}$, the error estimate would improve and give a relative error of order $\varepsilon$:

$$ \varepsilon^{-2}\|u_{\varepsilon} - \varepsilon^2 u_2\|_{L^2(\Omega^{\varepsilon})} + \varepsilon^{-1}|u_{\varepsilon} - \varepsilon^2 u_2|_{H^1(\Omega^{\varepsilon})} + \|p_{\varepsilon} - p^*\|_{L^2(\Omega^{\varepsilon})} \leq C\varepsilon\|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} $$

The proof of Theorem 4.1 is given in Subsection 4.5. Beforehand, it requires three technical lemmas concerning (i) solving a non-homogeneous divergence equation in a porous domain (Lemma 4.9), (ii) correcting non-homogeneous boundary conditions (Lemma 4.13) and (iii) representing periodic divergence-free fields (Lemma 4.15).
## 4.2 Correction of the divergence condition

**Notation 4.4.** In what follows, for any $\phi \in H_0^1(\Omega^\varepsilon)$, we also denote by $\phi$ its extension by zero in $\Omega \setminus \Omega^\varepsilon$, which belongs to $H_0^1(\Omega)$.

**Definition 4.5.** We note $\pi_i^\varepsilon$ the linear continuous invertible application, composed of a translation and an homothety of ratio $\frac{1}{\varepsilon}$, which maps any cell $Y_i^\varepsilon$ onto the reference cell Y:

$$ \pi_i^\varepsilon: \begin{array}{ccc} Y_i^\varepsilon & \to & Y \\ x & \mapsto & y = \frac{x}{\varepsilon} + \text{translation.} \end{array} \qquad (18) $$

First, we recall the Poincaré inequality in perforated domains (for a proof, see [10]).

**Lemma 4.6.** Assume that the hypotheses on the unit cell Y (Section 2.1) hold true, then there exists a constant $C > 0$, independent of $\varepsilon$, such that, for any $\phi \in H_0^1(\Omega^\varepsilon)$,

$$ \|\phi\|_{L^2(\Omega^\varepsilon)} \leq C\varepsilon|\phi|_{H^1(\Omega^\varepsilon)}. $$

We now recall a result [15, 1] on a restriction operator for velocities defined on the full domain $\Omega$ which are restricted to the porous domain $\Omega^\varepsilon$.

**Lemma 4.7.** Assume that the hypotheses on the unit cell Y (Section 2.1) hold true. Then there exists a linear continuous operator $R_\varepsilon$ such that

1. $R_\varepsilon \in \mathcal{L}(H_0^1(\Omega)^d; H_0^1(\Omega^\varepsilon)^d)$;

2. $u|_{\Omega^\varepsilon} \in H_0^1(\Omega^\varepsilon)^d$ implies $R_\varepsilon u = u|_{\Omega^\varepsilon}$ in $\Omega^\varepsilon$;

3. $\text{div } u = 0$ in $\Omega \setminus \Omega^\varepsilon$ implies $\text{div}(R_\varepsilon u) = \text{div } u$ in $\Omega^\varepsilon$;

4. there exists a constant $C$, which does not depend on $\varepsilon$, such that, for any $u \in H_0^1(\Omega)^d$, we have:

$$ \|R_\varepsilon u\|_{L^2(\Omega^\varepsilon)} + \varepsilon \|\nabla(R_\varepsilon u)\|_{L^2(\Omega^\varepsilon)} \leq C (\|u\|_{L^2(\Omega)} + \varepsilon \|\nabla u\|_{L^2(\Omega)}). $$

**Remark 4.8.** Lemma 4.7 was first proved in [15] for the case of a solid part $Y_S$ strictly included in the cell Y (see Figure 3a) and was then extended in [1] for the more general case of a solid part $Y_S$ which touches the boundary of Y (see Figure 3b).

*Proof.* Since we shall need specific properties of the operator $R_\varepsilon$ in the sequel, we briefly recall how to construct this operator. We start with the simpler case of isolated obstacles, namely when the solid part $Y_S$ is strictly included in the cell Y (the general case is presented afterwards). Consider a domain $Y_M$, strictly contained in Y, such that $\partial Y_M$ is smooth and encloses $Y_S$ as presented in Figure 5.

<figure><img src="image_5.png" /><figcaption>Figure 5: Schematic description of the domain $Y_M$ around $Y_S$</figcaption></figure>
## 4.2 Correction of the divergence condition

**Notation 4.4.** In what follows, for any $\phi \in H_0^1(\Omega^\epsilon)$, we also denote by $\phi$ its extension by zero in $\Omega \setminus \Omega^\epsilon$, which belongs to $H_0^1(\Omega)$.

**Definition 4.5.** We note $\pi_i^\epsilon$ the linear continuous invertible application, composed of a translation and an homothety of ratio $\frac{1}{\epsilon}$, which maps any cell $Y_i^\epsilon$ onto the reference cell Y:

$$ \pi_i^\epsilon: \begin{array}{ccc} Y_i^\epsilon & \to & Y \\ x & \mapsto & y = \frac{x}{\epsilon} + \text{translation.} \end{array} \qquad (18) $$

First, we recall the Poincaré inequality in perforated domains (for a proof, see [10]).

**Lemma 4.6.** Assume that the hypotheses on the unit cell Y (Section 2.1) hold true, then there exists a constant $C > 0$, independent of $\epsilon$, such that, for any $\phi \in H_0^1(\Omega^\epsilon)$,

$$ \|\phi\|_{L^2(\Omega^\epsilon)} \leq C\epsilon|\phi|_{H^1(\Omega^\epsilon)}. $$

We now recall a result [15, 1] on a restriction operator for velocities defined on the full domain $\Omega$ which are restricted to the porous domain $\Omega^\epsilon$.

**Lemma 4.7.** Assume that the hypotheses on the unit cell Y (Section 2.1) hold true. Then there exists a linear continuous operator $R_\epsilon$ such that

1. $R_\epsilon \in \mathcal{L}(H_0^1(\Omega)^d; H_0^1(\Omega^\epsilon)^d)$;

2. $u|_{\Omega^\epsilon} \in H_0^1(\Omega^\epsilon)^d$ implies $R_\epsilon u = u|_{\Omega^\epsilon}$ in $\Omega^\epsilon$;

3. $\text{div } u = 0$ in $\Omega \setminus \Omega^\epsilon$ implies $\text{div}(R_\epsilon u) = \text{div } u$ in $\Omega^\epsilon$;

4. there exists a constant $C$, which does not depend on $\epsilon$, such that, for any $u \in H_0^1(\Omega)^d$, we have:

$$ \|R_\epsilon u\|_{L^2(\Omega^\epsilon)} + \epsilon \|\nabla(R_\epsilon u)\|_{L^2(\Omega^\epsilon)} \leq C (\|u\|_{L^2(\Omega)} + \epsilon \|\nabla u\|_{L^2(\Omega)}). $$

**Remark 4.8.** Lemma 4.7 was first proved in [15] for the case of a solid part $Y_S$ strictly included in the cell Y (see Figure 3a) and was then extended in [1] for the more general case of a solid part $Y_S$ which touches the boundary of Y (see Figure 3b).

*Proof.* Since we shall need specific properties of the operator $R_\epsilon$ in the sequel, we briefly recall how to construct this operator. We start with the simpler case of isolated obstacles, namely when the solid part $Y_S$ is strictly included in the cell Y (the general case is presented afterwards). Consider a domain $Y_M$, strictly contained in Y, such that $\partial Y_M$ is smooth and encloses $Y_S$ as presented in Figure 5.

<figure><img src="image_7.png" /><figcaption>Figure 5: Schematic description of the domain $Y_M$ around $Y_S$</figcaption></figure>
Then for a given $\mathbf{u} \in H^1(Y)^d$, there exists $\mathbf{v} \in H^1(Y_M \setminus Y_S)^d$, $q \in L^2(Y_M \setminus Y_S)/\mathbb{R}$ satisfying:

$$
\left\{
\begin{array}{ll}
-\Delta \mathbf{v} + \nabla p = -\Delta \mathbf{u} & \text{in } Y_M \setminus Y_S, \\
\operatorname{div} \mathbf{v} = \operatorname{div} \mathbf{u} + \frac{1}{|Y_M \setminus Y_S|} \int_{Y_S} \operatorname{div} \mathbf{u} & \text{in } Y_M \setminus Y_S, \\
\mathbf{v} = \mathbf{u} & \text{on } \partial Y_M, \\
\mathbf{v} = \mathbf{0} & \text{on } \partial Y_S.
\end{array}
\right.
$$

This problem is well posed, since the compatibility condition can easily be checked on the data,
namely

$$
\int_{Y_M \setminus Y_S} \operatorname{div} \mathbf{v} = \int_{\partial(Y_M \setminus Y_S)} \mathbf{v} \cdot \mathbf{n}.
$$

Now, a restriction operator $R$ is constructed as follow: for any $\mathbf{u} \in H^1(Y)^d$,

$$
Ru = \begin{cases} \mathbf{u} & \text{in } Y \setminus Y_M, \\ \mathbf{v} & \text{in } Y_M \setminus Y_S, \\ \mathbf{0} & \text{in } Y_S. \end{cases}
$$

Clearly, $R$ is linear continuous from $H^1(Y)^d$ to its subspace of functions vanishing on $Y_S$, satisfying

$$
\|Ru\|_{H^1(Y)} \le C \|u\|_{H^1(Y)}, \quad (19)
$$

and such that, if $\mathbf{u} = 0$ in $Y_S$, then $Ru = \mathbf{u}$, and, if $\operatorname{div} \mathbf{u} = 0$ in $Y_S$, then $\operatorname{div} Ru = \operatorname{div} \mathbf{u}$ in $Y$. Then the operator $R_\varepsilon$ is defined by applying the mapping (18) to the operator $R$ and, with standard scaling properties, it comes from (19):

$$
\| R_{\varepsilon} \mathbf{u} \|_{L^2(\Omega^{\varepsilon})}^2 + \varepsilon^2 \| \nabla (R_{\varepsilon} \mathbf{u}) \|_{L^2(\Omega^{\varepsilon})}^2 \le C \left( \| \mathbf{u} \|_{L^2(\Omega)}^2 + \varepsilon^2 \| \nabla \mathbf{u} \|_{L^2(\Omega)}^2 \right),
$$

which implies that $R_\epsilon$ satisfies all four properties of Lemma 4.7.

We now turn to the case of connected obstacles, namely when $Y_S$ touches the cell bound-
ary $\partial Y$ as presented in Figure 3b. The construction of $R_\epsilon$ is slightly more involved. First, we
define a linear continuous operator $Q_\epsilon \in \mathcal{L}(H_0^1(\Omega)^d; H_0^1(\Omega^\epsilon)^d)$ by

$$
\left\{
\begin{array}{rcl@{\quad}l}
-\Delta(Q_{\varepsilon}\mathbf{u}) & = & -\Delta\mathbf{u} & \text{in } \Omega^{\varepsilon} \cap C_{1}^{\varepsilon}, \\
Q_{\varepsilon}\mathbf{u} & = & \mathbf{0} & \text{on } \partial\Omega^{\varepsilon} \cap C_{1}^{\varepsilon}, \\
Q_{\varepsilon}\mathbf{u} & = & \mathbf{u} & \text{on } \partial C_{1}^{\varepsilon}, \\
Q_{\varepsilon}\mathbf{u} & = & \mathbf{u} & \text{in } \Omega \setminus C_{1}^{\varepsilon},
\end{array}
\right.
\tag{20}
$$

where $C_1^\epsilon$ is defined by Assumption 2.1 as the union of entire cells inside $\Omega$. Clearly, $Q_\epsilon \mathbf{u} = \mathbf{u}$ if and only if $\mathbf{u}|_{\Omega^\epsilon} \in H_0^1(\Omega^\epsilon)^d$. Thanks to the definition of $C_1^\epsilon$ and $C_2^\epsilon$, there exists a smooth cut-off function $\zeta^\epsilon$ such that $\zeta^\epsilon(x) = 1$ on $\partial C_1^\epsilon$, $\zeta^\epsilon(x) = 0$ in $C_2^\epsilon$ and $\|\nabla\zeta^\epsilon\|_{L^\infty(C_1^\epsilon)} \le \frac{C}{\epsilon}$. One can rewrite (20) for the new unknown $\mathbf{w}^\epsilon = Q_\epsilon \mathbf{u} - \zeta^\epsilon \mathbf{u}$ as

$$
\left\{
\begin{array}{rcl}
-\Delta \mathbf{w}^{\epsilon} & = & -\Delta((1-\zeta^{\epsilon})\mathbf{u}) \quad \text{in } \Omega^{\epsilon} \cap C_{1}^{\epsilon}, \\
\mathbf{w}^{\epsilon} & = & \mathbf{0} \quad \text{on } \partial\Omega^{\epsilon} \cap C_{1}^{\epsilon}, \\
\mathbf{w}^{\epsilon} & = & \mathbf{0} \quad \text{on } \partial C_{1}^{\epsilon}.
\end{array}
\right.
\qquad (21)
$$

Multiplying (21) by $\mathbf{w}^\epsilon$ leads to an a priori estimate, which, combined with the Poincaré in-
equality (Lemma 4.6), implies that there exists a constant $C$, which does not depend on $\epsilon$, such
that, for any $\mathbf{u} \in H_0^1(\Omega)^d$, we have:

$$
\| Q_{\epsilon} \mathbf{u} \|_{L^2(\Omega^{\epsilon})} + \epsilon \| \nabla (Q_{\epsilon} \mathbf{u}) \|_{L^2(\Omega^{\epsilon})} \le C ( \| \mathbf{u} \|_{L^2(\Omega)} + \epsilon \| \nabla \mathbf{u} \|_{L^2(\Omega)} ).
$$
Then for a given $\mathbf{u} \in H^1(Y)^d$, there exists $\mathbf{v} \in H^1(Y_M \setminus Y_S)^d$, $q \in L^2(Y_M \setminus Y_S)/\mathbb{R}$ satisfying:

$$
\left\{
\begin{array}{ll}
-\Delta \mathbf{v} + \nabla p = -\Delta \mathbf{u} & \text{in } Y_M \setminus Y_S, \\
\operatorname{div} \mathbf{v} = \operatorname{div} \mathbf{u} + \frac{1}{|Y_M \setminus Y_S|} \int_{Y_S} \operatorname{div} \mathbf{u} & \text{in } Y_M \setminus Y_S, \\
\mathbf{v} = \mathbf{u} & \text{on } \partial Y_M, \\
\mathbf{v} = \mathbf{0} & \text{on } \partial Y_S.
\end{array}
\right.
$$

This problem is well posed, since the compatibility condition can easily be checked on the data,
namely

$$
\int_{Y_M \setminus Y_S} \operatorname{div} \mathbf{v} = \int_{\partial(Y_M \setminus Y_S)} \mathbf{v} \cdot \mathbf{n}.
$$

Now, a restriction operator $R$ is constructed as follow: for any $\mathbf{u} \in H^1(Y)^d$,

$$
Ru = \begin{cases} \mathbf{u} & \text{in } Y \setminus Y_M, \\ \mathbf{v} & \text{in } Y_M \setminus Y_S, \\ \mathbf{0} & \text{in } Y_S. \end{cases}
$$

Clearly, $R$ is linear continuous from $H^1(Y)^d$ to its subspace of functions vanishing on $Y_S$, satisfying

$$
\|Ru\|_{H^1(Y)} \leq C \|u\|_{H^1(Y)}, \quad (19)
$$

and such that, if $\mathbf{u} = 0$ in $Y_S$, then $Ru = \mathbf{u}$, and, if $\operatorname{div} \mathbf{u} = 0$ in $Y_S$, then $\operatorname{div} Ru = \operatorname{div} \mathbf{u}$ in $Y$. Then the operator $R_\varepsilon$ is defined by applying the mapping (18) to the operator $R$ and, with standard scaling properties, it comes from (19):

$$
\| R_{\varepsilon} \mathbf{u} \|_{L^2(\Omega^{\varepsilon})}^2 + \varepsilon^2 \| \nabla (R_{\varepsilon} \mathbf{u}) \|_{L^2(\Omega^{\varepsilon})}^2 \le C \left( \| \mathbf{u} \|_{L^2(\Omega)}^2 + \varepsilon^2 \| \nabla \mathbf{u} \|_{L^2(\Omega)}^2 \right),
$$

which implies that $R_\epsilon$ satisfies all four properties of Lemma 4.7.

We now turn to the case of connected obstacles, namely when $Y_S$ touches the cell bound-
ary $\partial Y$ as presented in Figure 3b. The construction of $R_\epsilon$ is slightly more involved. First, we
define a linear continuous operator $Q_\epsilon \in \mathcal{L}(H_0^1(\Omega)^d; H_0^1(\Omega^\epsilon)^d)$ by

$$
\left\{
\begin{array}{rcl@{\quad}l}
-\Delta(Q_{\varepsilon}\mathbf{u}) & = & -\Delta\mathbf{u} & \text{in } \Omega^{\varepsilon} \cap C_{1}^{\varepsilon}, \\
Q_{\varepsilon}\mathbf{u} & = & \mathbf{0} & \text{on } \partial\Omega^{\varepsilon} \cap C_{1}^{\varepsilon}, \\
Q_{\varepsilon}\mathbf{u} & = & \mathbf{u} & \text{on } \partial C_{1}^{\varepsilon}, \\
Q_{\varepsilon}\mathbf{u} & = & \mathbf{u} & \text{in } \Omega \setminus C_{1}^{\varepsilon},
\end{array}
\right.
\qquad (20)
$$

where $C_1^\epsilon$ is defined by Assumption 2.1 as the union of entire cells inside $\Omega$. Clearly, $Q_\epsilon \mathbf{u} = \mathbf{u}$ if and only if $\mathbf{u}|_{\Omega^\epsilon} \in H_0^1(\Omega^\epsilon)^d$. Thanks to the definition of $C_1^\epsilon$ and $C_2^\epsilon$, there exists a smooth cut-off function $\zeta^\epsilon$ such that $\zeta^\epsilon(x) = 1$ on $\partial C_1^\epsilon$, $\zeta^\epsilon(x) = 0$ in $C_2^\epsilon$ and $\|\nabla\zeta^\epsilon\|_{L^\infty(C_1^\epsilon)} \le \frac{C}{\epsilon}$. One can rewrite (20) for the new unknown $\mathbf{w}^\epsilon = Q_\epsilon \mathbf{u} - \zeta^\epsilon \mathbf{u}$ as

$$
\left\{
\begin{array}{rcl@{\quad}l}
-\Delta \mathbf{w}^{\epsilon} & = & -\Delta((1-\zeta^{\epsilon})\mathbf{u}) & \text{in } \Omega^{\epsilon} \cap C_{1}^{\epsilon}, \\
\mathbf{w}^{\epsilon} & = & \mathbf{0} & \text{on } \partial\Omega^{\epsilon} \cap C_{1}^{\epsilon}, \\
\mathbf{w}^{\epsilon} & = & \mathbf{0} & \text{on } \partial C_{1}^{\epsilon}.
\end{array}
\right.
\qquad (21)
$$

Multiplying (21) by $\mathbf{w}^\epsilon$ leads to an a priori estimate, which, combined with the Poincaré in-
equality (Lemma 4.6), implies that there exists a constant $C$, which does not depend on $\epsilon$, such
that, for any $\mathbf{u} \in H_0^1(\Omega)^d$, we have:

$$
\| Q_{\epsilon} \mathbf{u} \|_{L^2(\Omega^{\epsilon})} + \epsilon \| \nabla (Q_{\epsilon} \mathbf{u}) \|_{L^2(\Omega^{\epsilon})} \le C ( \| \mathbf{u} \|_{L^2(\Omega)} + \epsilon \| \nabla \mathbf{u} \|_{L^2(\Omega)} ).
$$
Then, for each cell $Y_i^\epsilon$, we define a linear operator $Q_i$, from $H^1(Y)^d$ to its subspace of functions vanishing on $Y_S$, as follows: for any $\mathbf{u} \in H_0^1(\Omega)^d$,

$$ (Q_i(\mathbf{u} \circ (\pi_i^\epsilon)^{-1})) \circ \pi_i^\epsilon = Q_\epsilon \mathbf{u} \text{ in } Y_i^\epsilon. $$

Recall from Assumption 2.1 that we defined two sets of indices $\mathcal{I}_1$ and $\mathcal{I}_2$, corresponding to the cells $Y_i^\epsilon$ which cover $C_1^\epsilon$ and $C_2^\epsilon$, respectively. Define a third set $\mathcal{I}_0 = \{i \in \mathbb{Z}^d \text{ s.t. } Y_i^\epsilon \cap \Omega \neq \emptyset\}$, so that $\mathcal{I}_2 \subset \mathcal{I}_1 \subset \mathcal{I}_0$. We now define a family of restriction operators $R_i$, depending on which set $\mathcal{I}_2, \mathcal{I}_1, \mathcal{I}_0$ the index $i$ belongs to.

For $i \in \mathcal{I}_2$ (meaning that the corresponding cell contains a solid obstacle), $R_i$ is defined, for any $\mathbf{u} \in H^1(Y)^d$, by

$$ \left\{ \begin{array}{ll} -\Delta R_i \mathbf{u} + \nabla p = -\Delta \mathbf{u} & \text{in } Y_F, \\ \operatorname{div} R_i \mathbf{u} = \operatorname{div} \mathbf{u} + \frac{1}{|Y_F|} \int_{Y_S} \operatorname{div} \mathbf{u} & \text{in } Y_F, \\ R_i \mathbf{u} = Q_i \mathbf{u} + \phi_k \left( \int_{\Sigma_k} (\mathbf{u} - Q_i \mathbf{u}) \cdot \mathbf{e}_k \right) \mathbf{e}_k & \text{on } \Sigma_k \cap \bar{Y}_F, \\ R_i \mathbf{u} = \mathbf{0} & \text{on } \partial Y_S, \end{array} \right. \quad (22) $$

where $(\Sigma_k)_{1 \le |k| \le d}$ are the $2d$ faces of the unit cube $(\Sigma_k$ and $\Sigma_{-k}$ are opposite faces), with unit normal vector $\mathbf{e}_k$, and $(\phi_k)_{1 \le |k| \le d}$ is a family of smooth, non-negative functions, compactly supported in the periodic repetition of $\bar{Y}_F$, satisfying

$$ \int_{\Sigma_k} \phi_k = 1, \quad \phi_{k|\Sigma_k} \equiv \phi_{-k|\Sigma_{-k}} \quad \text{and} \quad \phi_{k|\Sigma_{k'}} \equiv 0 \text{ if } k \neq k'. $$

The existence of such a family of functions $\phi_k$ is guaranteed by the assumptions on the unit cell in Subsection 2.1, which imply that the surface measure of $\Sigma_k \cap \bar{Y}_F$ is non-zero (see [1] for more details). Problem (22) is well posed, since the data of the trace of $R_i \mathbf{u}$ belongs to $H^{1/2}(Y_F)^d$ and the compatibility condition on the divergence data is satisfied, namely

$$ \int_{Y_F} \operatorname{div} R_i \mathbf{u} = \int_{\partial Y_F} R_i \mathbf{u} \cdot \mathbf{n}. $$

For $i \in \mathcal{I}_1 \setminus \mathcal{I}_2$ (meaning that the corresponding cell does not contain a solid obstacle but one of its face may touch an obstacle), $R_i$ is defined, for any $\mathbf{u} \in H^1(Y)^d$, by

$$ \left\{ \begin{array}{ll} -\Delta R_i \mathbf{u} + \nabla p = -\Delta \mathbf{u} & \text{in } Y, \\ \operatorname{div} R_i \mathbf{u} = \operatorname{div} \mathbf{u} & \text{in } Y, \\ R_i \mathbf{u} = Q_i \mathbf{u} + \phi_k \left( \int_{\Sigma_k} (\mathbf{u} - Q_i \mathbf{u}) \cdot \mathbf{e}_k \right) \mathbf{e}_k & \text{on } \Sigma_k, \\ R_i \mathbf{u} = \mathbf{u} & \text{on } \tilde{\Sigma}_k, \end{array} \right. \quad (23) $$

where $\Sigma_k$ are the faces of the cell which touch an obstacle and $\tilde{\Sigma}_k$ are the faces of the cell which do not touch any obstacle. Problem (23) is well posed, since the data of the trace of $R_i \mathbf{u}$ belongs to $H^{1/2}(Y)^d$ and the compatibility condition on the divergence data is satisfied, namely

$$ \int_Y \operatorname{div} R_i \mathbf{u} = \int_{\partial Y} R_i \mathbf{u} \cdot \mathbf{n}. $$

For $i \in \mathcal{I}_0 \setminus \mathcal{I}_1$ (implying that the corresponding cell does not contain a solid obstacle and none of its face touches an obstacle), $R_i$ is simply defined as $R_i \mathbf{u} = \mathbf{u}$. Note that this case corresponds to cells $Y_i^\epsilon$ which are cut by the boundary $\partial \Omega$. Inside such a cell, with a possibly wild geometry, it is not clear that a Stokes problem of the type of (23) satisfies a priori estimates, independent of the geometry (this is precisely why the two sets of indices $\mathcal{I}_1$ and $\mathcal{I}_2$ have been introduced).
Then, for each cell $Y_i^\epsilon$, we define a linear operator $Q_i$, from $H^1(Y)^d$ to its subspace of functions vanishing on $Y_S$, as follows: for any $u \in H_0^1(\Omega)^d$,

$$ (Q_i (u \circ (\pi_i^\epsilon)^{-1})) \circ \pi_i^\epsilon = Q_\epsilon u \text{ in } Y_i^\epsilon. $$

Recall from Assumption 2.1 that we defined two sets of indices $\mathcal{I}_1$ and $\mathcal{I}_2$, corresponding to the cells $Y_i^\epsilon$ which cover $C_1^\epsilon$ and $C_2^\epsilon$, respectively. Define a third set $\mathcal{I}_0 = \{i \in \mathbb{Z}^d \text{ s.t. } Y_i^\epsilon \cap \Omega \neq \emptyset\}$, so that $\mathcal{I}_2 \subset \mathcal{I}_1 \subset \mathcal{I}_0$. We now define a family of restriction operators $R_i$, depending on which set $\mathcal{I}_2, \mathcal{I}_1, \mathcal{I}_0$ the index $i$ belongs to.

For $i \in \mathcal{I}_2$ (meaning that the corresponding cell contains a solid obstacle), $R_i$ is defined, for any $u \in H^1(Y)^d$, by

$$ \left\{ \begin{array}{ll} -\Delta R_i u + \nabla p = -\Delta u & \text{in } Y_F, \\ \operatorname{div} R_i u = \operatorname{div} u + \frac{1}{|Y_F|} \int_{Y_S} \operatorname{div} u & \text{in } Y_F, \\ R_i u = Q_i u + \phi_k \left( \int_{\Sigma_k} (u - Q_i u) \cdot e_k \right) e_k & \text{on } \Sigma_k \cap \bar{Y}_F, \\ R_i u = \mathbf{0} & \text{on } \partial Y_S, \end{array} \right. \quad (22) $$

where $(\Sigma_k)_{1 \le |k| \le d}$ are the $2d$ faces of the unit cube $(\Sigma_k$ and $\Sigma_{-k}$ are opposite faces), with unit normal vector $e_k$, and $(\phi_k)_{1 \le |k| \le d}$ is a family of smooth, non-negative functions, compactly supported in the periodic repetition of $\bar{Y}_F$, satisfying

$$ \int_{\Sigma_k} \phi_k = 1, \quad \phi_{k|\Sigma_k} \equiv \phi_{-k|\Sigma_{-k}} \quad \text{and} \quad \phi_{k|\Sigma_{k'}} \equiv 0 \text{ if } k \neq k'. $$

The existence of such a family of functions $\phi_k$ is guaranteed by the assumptions on the unit cell in Subsection 2.1, which imply that the surface measure of $\Sigma_k \cap \bar{Y}_F$ is non-zero (see [1] for more details). Problem (22) is well posed, since the data of the trace of $R_i u$ belongs to $H^{1/2}(Y_F)^d$ and the compatibility condition on the divergence data is satisfied, namely

$$ \int_{Y_F} \operatorname{div} R_i u = \int_{\partial Y_F} R_i u \cdot n. $$

For $i \in \mathcal{I}_1 \setminus \mathcal{I}_2$ (meaning that the corresponding cell does not contain a solid obstacle but one of its face may touch an obstacle), $R_i$ is defined, for any $u \in H^1(Y)^d$, by

$$ \left\{ \begin{array}{ll} -\Delta R_i u + \nabla p = -\Delta u & \text{in } Y, \\ \operatorname{div} R_i u = \operatorname{div} u & \text{in } Y, \\ R_i u = Q_i u + \phi_k \left( \int_{\Sigma_k} (u - Q_i u) \cdot e_k \right) e_k & \text{on } \Sigma_k, \\ R_i u = \mathbf{u} & \text{on } \tilde{\Sigma}_k, \end{array} \right. \quad (23) $$

where $\Sigma_k$ are the faces of the cell which touch an obstacle and $\tilde{\Sigma}_k$ are the faces of the cell which do not touch any obstacle. Problem (23) is well posed, since the data of the trace of $R_i u$ belongs to $H^{1/2}(Y)^d$ and the compatibility condition on the divergence data is satisfied, namely

$$ \int_Y \operatorname{div} R_i u = \int_{\partial Y} R_i u \cdot n. $$

For $i \in \mathcal{I}_0 \setminus \mathcal{I}_1$ (implying that the corresponding cell does not contain a solid obstacle and none of its face touches an obstacle), $R_i$ is simply defined as $R_i u = u$. Note that this case corresponds to cells $Y_i^\epsilon$ which are cut by the boundary $\partial\Omega$. Inside such a cell, with a possibly wild geometry, it is not clear that a Stokes problem of the type of (23) satisfies a priori estimates, independent of the geometry (this is precisely why the two sets of indices $\mathcal{I}_1$ and $\mathcal{I}_2$ have been introduced).
A similar estimate to (19) holds true for all operators $R_i$ with a constant $C$ independent of $i$ and $\epsilon$.

Finally, the restriction operator $R_\epsilon \in \mathcal{L}(H_0^1(\Omega)^d; H_0^1(\Omega^\epsilon)^d)$ is defined, for any $\mathbf{u} \in H_0^1(\Omega)^d$, by

$$R_{\epsilon}\mathbf{u} = (R_i(\mathbf{u} \circ (\pi_i^{\epsilon})^{-1})) \circ \pi_i^{\epsilon} \text{ in } Y_i^{\epsilon}.$$

It is easy again to check that $R_\epsilon$ would satisfy all four properties of Lemma 4.7. Note however that, here, the definition of $R_\epsilon$ is not local in each cell since they are all coupled through the operator $Q_\epsilon$. $\square$

We now state a key lemma for our analysis. It was already presented in [11] in two dimensions for the case when $Y_S$ is strictly included in $Y$. It was also stated as Theorem 14.2 in [5] for the general case of any $Y_S$ and any space dimension but its proof was rather technical, decomposed in several exercices. Here, we present a new and detailed proof, using different arguments.

**Lemma 4.9.** For any $q \in L_0^2(\Omega^\epsilon)$ there exists $v_{d\epsilon} \in H_0^1(\Omega^\epsilon)^d$ such that

$$\operatorname{div} v_{d\epsilon} = q \text{ on } \Omega^\epsilon \text{ and } |v_{d\epsilon}|_{H^1(\Omega^\epsilon)} \le \frac{C}{\epsilon} \|q\|_{L^2(\Omega^\epsilon)} \quad (24)$$

where $C > 0$ is a constant independent of $\epsilon$.

*Proof.* For any $q \in L_0^2(\Omega^\epsilon)$ we still denote by $q$ its extension by zero inside the obstacles $B^\epsilon$. One can easily check that its average on $\Omega$ still vanishes and thus $q \in L_0^2(\Omega)$ with $q = 0$ on $B^\epsilon$. Using [9, corollary 2.4], there exists a (non unique) $\mathbf{w} \in H_0^1(\Omega)^d$ such that

$$\operatorname{div} \mathbf{w} = q \text{ on } \Omega \quad \text{and} \quad \|\mathbf{w}\|_{H^1(\Omega)} \le C \|q\|_{L^2(\Omega)}.$$

Since $q = 0$ on $B^\epsilon$, we have

$$\operatorname{div} \mathbf{w} = 0 \text{ in } \Omega \setminus \Omega^\epsilon.$$

Using the restriction operator $R_\epsilon$ defined in Lemma 4.7, we have

$$\operatorname{div}(R_{\epsilon}\mathbf{w}) = q \text{ in } \Omega^{\epsilon}.$$

Considering $\epsilon < 1$ in the fourth property of Lemma 4.7, it follows

$$\|R_{\epsilon}\mathbf{w}\|_{L^2(\Omega^{\epsilon})} + \epsilon \|\nabla R_{\epsilon}\mathbf{w}\|_{L^2(\Omega^{\epsilon})} \le C \|\mathbf{w}\|_{H^1(\Omega)}$$

and consequently

$$\epsilon \|\nabla(R_{\epsilon}\mathbf{w})\|_{L^2(\Omega^{\epsilon})} \le C \|\mathbf{w}\|_{H^1(\Omega)} \le C \|q\|_{L^2(\Omega^{\epsilon})}.$$

Therefore, $v_{d\epsilon} = R_\epsilon \mathbf{w}$ is a solution of (24). $\square$

## 4.3 Correction of the boundary condition

In this section, we build a divergence-free field which satisfies a given boundary condition. Recall the definition of the curl* operator, proposed in [6].

**Definition 4.10.** The operator curl* is defined as:

$$\operatorname{curl}^* : C^1(\mathbb{R}^d)^{d(d-1)/2} \to C^0(\mathbb{R}^d)^d \\ \mathbf{A} \mapsto \left( \sum_{j=1}^{i-1} \frac{\partial A_{ji}}{\partial x_j} - \sum_{j=i+1}^{d} \frac{\partial A_{ij}}{\partial x_j} \right)_{1 \le i \le d}$$

where $\mathbf{A} = (A_{ij})_{1 \le i < j \le d} \in \mathbb{R}^{d(d-1)/2}$.
A similar estimate to (19) holds true for all operators $R_i$ with a constant $C$ independent of $i$ and $\varepsilon$.

Finally, the restriction operator $R_\varepsilon \in \mathcal{L}(H_0^1(\Omega)^d; H_0^1(\Omega^\varepsilon)^d)$ is defined, for any $u \in H_0^1(\Omega)^d$, by

$$R_\varepsilon u = (R_i (u \circ (\pi_i^\varepsilon)^{-1})) \circ \pi_i^\varepsilon \text{ in } Y_i^\varepsilon.$$

It is easy again to check that $R_\varepsilon$ would satisfy all four properties of Lemma 4.7. Note however that, here, the definition of $R_\varepsilon$ is not local in each cell since they are all coupled through the operator $Q_\varepsilon$. $\square$

We now state a key lemma for our analysis. It was already presented in [11] in two dimensions for the case when $Y_S$ is strictly included in $Y$. It was also stated as Theorem 14.2 in [5] for the general case of any $Y_S$ and any space dimension but its proof was rather technical, decomposed in several exercices. Here, we present a new and detailed proof, using different arguments.

**Lemma 4.9.** For any $q \in L_0^2(\Omega^\varepsilon)$ there exists $v_{d\varepsilon} \in H_0^1(\Omega^\varepsilon)^d$ such that

$$\operatorname{div} v_{d\varepsilon} = q \text{ on } \Omega^\varepsilon \text{ and } |v_{d\varepsilon}|_{H^1(\Omega^\varepsilon)} \le \frac{C}{\varepsilon} \|q\|_{L^2(\Omega^\varepsilon)} \quad (24)$$

where $C > 0$ is a constant independent of $\varepsilon$.

*Proof.* For any $q \in L_0^2(\Omega^\varepsilon)$ we still denote by $q$ its extension by zero inside the obstacles $B^\varepsilon$. One can easily check that its average on $\Omega$ still vanishes and thus $q \in L_0^2(\Omega)$ with $q = 0$ on $B^\varepsilon$. Using [9, corollary 2.4], there exists a (non unique) $w \in H_0^1(\Omega)^d$ such that

$$\operatorname{div} w = q \text{ on } \Omega \quad \text{and} \quad \|w\|_{H^1(\Omega)} \le C \|q\|_{L^2(\Omega)}.$$

Since $q = 0$ on $B^\varepsilon$, we have

$$\operatorname{div} w = 0 \text{ in } \Omega \setminus \Omega^\varepsilon.$$

Using the restriction operator $R_\varepsilon$ defined in Lemma 4.7, we have

$$\operatorname{div}(R_\varepsilon w) = q \text{ in } \Omega^\varepsilon.$$

Considering $\varepsilon < 1$ in the fourth property of Lemma 4.7, it follows

$$\|R_\varepsilon w\|_{L^2(\Omega^\varepsilon)} + \varepsilon \|\nabla R_\varepsilon w\|_{L^2(\Omega^\varepsilon)} \le C \|w\|_{H^1(\Omega)}$$

and consequently

$$\varepsilon \|\nabla (R_\varepsilon w)\|_{L^2(\Omega^\varepsilon)} \le C \|w\|_{H^1(\Omega)} \le C \|q\|_{L^2(\Omega^\varepsilon)}.$$

Therefore, $v_{d\varepsilon} = R_\varepsilon w$ is a solution of (24). $\square$

## 4.3 Correction of the boundary condition

In this section, we build a divergence-free field which satisfies a given boundary condition. Recall the definition of the curl* operator, proposed in [6].

**Definition 4.10.** The operator curl* is defined as:

$$\operatorname{curl}^* : C^1(\mathbb{R}^d)^{d(d-1)/2} \to C^0(\mathbb{R}^d)^d \\ \mathbf{A} \mapsto \left( \sum_{j=1}^{i-1} \frac{\partial A_{ji}}{\partial x_j} - \sum_{j=i+1}^{d} \frac{\partial A_{ij}}{\partial x_j} \right)_{1 \le i \le d}$$

where $\mathbf{A} = (A_{ij})_{1 \le i < j \le d} \in \mathbb{R}^{d(d-1)/2}$.
A simple computation shows that the image of this curl* operator is divergence-free, namely $\operatorname{div} \operatorname{curl}^* \mathcal{A} = 0$ for any $\mathcal{A} \in C^1(\mathbb{R}^d)^{d(d-1)/2}$. This is an expected property since $\operatorname{curl}^*$ is an extension, to any dimension $d$, of the standard curl operator in 2-d and 3-d. We now state a preliminary lemma, also adapted from [6]

**Lemma 4.11.** Let $\Omega \subset \mathbb{R}^d$ be a bounded connected open set with a smooth boundary of class $C^{3,\alpha}$, for some $0 < \alpha < 1$. Let $g$ in $C^{1,\alpha}(\bar{\Omega})^d$ with $g \cdot n = 0$ on $\partial\Omega$. Then there exists $\mathcal{A}$ in $C^{2,\alpha}(\bar{\Omega})^{d(d-1)/2}$ such that

$$ \operatorname{curl}^*(\mathcal{A}) = g \text{ on } \partial\Omega \text{ and } \mathcal{A} = 0 \text{ on } \partial\Omega. $$

In other words, the vector field $\operatorname{curl}^*(\mathcal{A})$ is divergence-free and its trace on $\partial\Omega$ coincides with $g$.

*Proof.* [6, Lemma 1] ensures that there exists $\mathcal{A} \in C^{2,\alpha}(\bar{\Omega})^{d(d-1)/2}$, with $\mathcal{A} = 0$ on $\partial\Omega$, such that

$$ \nabla \mathcal{A}_{ij} = (g_j n_i - g_i n_j) n \text{ on } \partial \Omega, $$

where $n$ is the outward unit normal vector, $g_i$ and $n_i$ are respectively the $i$th coordinates of $g$ and $n$ in the canonical basis of $\mathbb{R}^d$. It remains to check that $\operatorname{curl}^*(\mathcal{A}) = g$ on $\partial\Omega$. Indeed, we have on $\partial\Omega$:

$$ \begin{aligned} (\operatorname{curl}^* \mathcal{A})_i &= \sum_{j=1}^{i-1} (g_j n_j - g_i n_j) n_j - \sum_{j=i+1}^{d} (g_j n_i - g_i n_j) n_j \\ &= \sum_{j=1, j \neq i}^{d} g_i n_j^2 - \sum_{j=1, j \neq i}^{d} g_j n_j n_i = \sum_{j=1}^{d} g_i n_j^2 = g_i, \end{aligned} $$

where we used the assumption that $g \cdot n = 0$ on $\partial\Omega$. $\square$

The next step is to localize in the vicinity of the boundary $\partial\Omega$ the divergence-free field which lifts a given boundary condition. To do this, we introduce the following classical cut-off function as in [3, 11].

**Definition 4.12 (Cut-off function).** For a small enough constant $c > 0$, define a neighborhood of the boundary $\partial\Omega$, denoted by $O^\varepsilon := \{x \in \Omega : \operatorname{dist}(x, \partial\Omega) \le c\varepsilon\}$ (band of width $c\varepsilon$ near $\partial\Omega$). Then we choose a cut-off function $\theta^\varepsilon \in C^\infty(\Omega)$ such that $\theta^\varepsilon(x) = 1$ in $O^\varepsilon/2$, $\theta^\varepsilon(x) = 0$ in $\Omega \setminus O^\varepsilon$, and

$$ \|\nabla \theta^\varepsilon\|_{L^\infty(\Omega)} \le \frac{C}{\varepsilon}, \quad \|\nabla^2 \theta^\varepsilon\|_{L^\infty(\Omega)} \le \frac{C}{\varepsilon^2}. $$

The main result of this section is the following lemma.

**Lemma 4.13.** For any $g \in C^{1,\alpha}(\bar{\Omega})^d$ with $\operatorname{div} g = 0$ in $\Omega$, $g \cdot n = 0$ on $\partial\Omega$, there exists $v_{b\varepsilon} \in H^1(\Omega)^d$ such that $\operatorname{supp} v_{b\varepsilon} \subset O^\varepsilon := \{x \in \Omega : \operatorname{dist}(x, \partial\Omega) < c\varepsilon\}$ and

$$ v_{b\varepsilon} = g \text{ on } \partial\Omega, \operatorname{div} v_{b\varepsilon} = 0 \text{ in } \Omega \text{ and } |v_{b\varepsilon}|_{H^1(\Omega)} \le \frac{C}{\sqrt{\varepsilon}} \|g\|_{C^{1,\alpha}(\Omega)} $$

where $C > 0$ is a constant independent of $\varepsilon$.

**Remark 4.14.** Lemma 4.13 was already proved in [11] in two dimensions. Here, we extend this result to any dimension.
A simple computation shows that the image of this curl* operator is divergence-free, namely div curl* $\mathcal{A}$ = 0 for any $\mathcal{A} \in C^1(\mathbb{R}^d)^{d(d-1)/2}$. This is an expected property since curl* is an extension, to any dimension $d$, of the standard curl operator in 2-d and 3-d. We now state a preliminary lemma, also adapted from [6]

**Lemma 4.11.** Let $\Omega \subset \mathbb{R}^d$ be a bounded connected open set with a smooth boundary of class $C^{3,\alpha}$, for some $0 < \alpha < 1$. Let $g$ in $C^{1,\alpha}(\bar{\Omega})^d$ with $g \cdot n = 0$ on $\partial\Omega$. Then there exists $\mathcal{A}$ in $C^{2,\alpha}(\bar{\Omega})^{d(d-1)/2}$ such that

$$ \mathrm{curl}^*(\mathcal{A}) = g \text{ on } \partial\Omega \text{ and } \mathcal{A} = 0 \text{ on } \partial\Omega. $$

In other words, the vector field $\mathrm{curl}^*(\mathcal{A})$ is divergence-free and its trace on $\partial\Omega$ coincides with $g$.

*Proof.* [6, Lemma 1] ensures that there exists $\mathcal{A} \in C^{2,\alpha}(\bar{\Omega})^{d(d-1)/2}$, with $\mathcal{A} = 0$ on $\partial\Omega$, such that

$$ \nabla \mathcal{A}_{ij} = (g_j n_i - g_i n_j) n \text{ on } \partial \Omega, $$

where $n$ is the outward unit normal vector, $g_i$ and $n_i$ are respectively the $i$th coordinates of $g$ and $n$ in the canonical basis of $\mathbb{R}^d$. It remains to check that $\mathrm{curl}^*(\mathcal{A}) = g$ on $\partial\Omega$. Indeed, we have on $\partial\Omega$:

$$ \begin{aligned} (\mathrm{curl}^* \mathcal{A})_i &= \sum_{j=1}^{i-1} (g_i n_j - g_j n_i) n_j - \sum_{j=i+1}^{d} (g_j n_i - g_i n_j) n_j \\ &= \sum_{j=1, j \neq i}^{d} g_i n_j^2 - \sum_{j=1, j \neq i}^{d} g_j n_i n_j = \sum_{j=1}^{d} g_i n_j^2 = g_i, \end{aligned} $$

where we used the assumption that $g \cdot n = 0$ on $\partial\Omega$. $\square$

The next step is to localize in the vicinity of the boundary $\partial\Omega$ the divergence-free field which lifts a given boundary condition. To do this, we introduce the following classical cut-off function as in [3, 11].

**Definition 4.12 (Cut-off function).** For a small enough constant $c > 0$, define a neighborhood of the boundary $\partial\Omega$, denoted by $O^\epsilon := \{x \in \Omega : \mathrm{dist}(x, \partial\Omega) \le c\epsilon\}$ (band of width $c\epsilon$ near $\partial\Omega$). Then we choose a cut-off function $\theta^\epsilon \in C^\infty(\Omega)$ such that $\theta^\epsilon(x) = 1$ in $O^\epsilon/2$, $\theta^\epsilon(x) = 0$ in $\Omega \setminus O^\epsilon$, and

$$ \|\nabla \theta^\epsilon\|_{L^\infty(\Omega)} \le \frac{C}{\epsilon}, \quad \|\nabla^2 \theta^\epsilon\|_{L^\infty(\Omega)} \le \frac{C}{\epsilon^2}. $$

The main result of this section is the following lemma.

**Lemma 4.13.** For any $g \in C^{1,\alpha}(\bar{\Omega})^d$ with $\mathrm{div} g = 0$ in $\Omega$, $g \cdot n = 0$ on $\partial\Omega$, there exists $v_{b\epsilon} \in H^1(\Omega)^d$ such that $\mathrm{supp} v_{b\epsilon} \subset O^\epsilon := \{x \in \Omega : \mathrm{dist}(x, \partial\Omega) < c\epsilon\}$ and

$$ v_{b\epsilon} = g \text{ on } \partial\Omega, \quad \mathrm{div} \, v_{b\epsilon} = 0 \text{ in } \Omega \text{ and } |v_{b\epsilon}|_{H^1(\Omega)} \le \frac{C}{\sqrt{\epsilon}} \|g\|_{C^{1,\alpha}(\Omega)} $$

where $C > 0$ is a constant independent of $\epsilon$.

**Remark 4.14.** Lemma 4.13 was already proved in [11] in two dimensions. Here, we extend this result to any dimension.
Proof. Let $g \in C^{1,\alpha}(\bar{\Omega})^d$ with $\mathrm{div} g = 0$ in $\Omega$, $g \cdot n = 0$ on $\partial\Omega$. By virtue of Lemma 4.11, there exists a function $\mathcal{A} \in C^{2,\alpha}(\bar{\Omega})^{d(d-1)/2}$, such that $\mathrm{curl}^*(\mathcal{A}) = g$ on $\partial\Omega$. Define:

$$v_{b\epsilon} = \mathrm{curl}^*(\theta^{\epsilon} \mathcal{A}),$$

where $\theta^\epsilon$ is the cut-off function introduced in Definition 4.12. Obviously, $\mathrm{div} v_{b\epsilon} = 0$ in $\Omega$. Then, $v_{b\epsilon} = g$ on $\partial\Omega$ since, in the vicinity $O^\epsilon/2$ of $\partial\Omega$, we recall that $\theta^\epsilon \equiv 1$. Finally, it remains to bound $|v_{b\epsilon}|_{H^1(\Omega)}$. We write

$$\mathrm{curl}^*(\theta^{\epsilon} \mathcal{A}) = \theta^{\epsilon}\mathrm{curl}^*(\mathcal{A}) + \Gamma(\mathcal{A}, \nabla\theta^{\epsilon}),$$

where $\Gamma(\mathcal{A}, \nabla\theta^\epsilon)$ is defined as:

$$\Gamma(\mathcal{A}, \nabla\theta^{\epsilon})_i = \sum_{j=1}^{i-1} \mathcal{A}_{ji} \frac{\partial\theta^{\epsilon}}{\partial x_j} - \sum_{j=i+1}^{d} \mathcal{A}_{ij} \frac{\partial\theta^{\epsilon}}{\partial x_j}.$$

It follows therefore:

$$|v_{b\epsilon}|_{H^1(\Omega)} \le \|\theta^\epsilon \nabla \mathrm{curl}^*(\mathcal{A})\|_{L^2(O^\epsilon)} + \|\nabla\theta^\epsilon \mathrm{curl}^*(\mathcal{A})\|_{L^2(O^\epsilon)} + \|\nabla\Gamma(\mathcal{A}, \nabla\theta^\epsilon)\|_{L^2(O^\epsilon)} \\ \le C \left( \sqrt{\epsilon} \|\nabla\nabla\mathcal{A}\|_{L^\infty(O^\epsilon)} + \frac{1}{\sqrt{\epsilon}} \|\nabla\mathcal{A}\|_{L^\infty(O^\epsilon)} \right) + \|\nabla\Gamma(\mathcal{A}, \nabla\theta^\epsilon)\|_{L^2(O^\epsilon)},$$

since meas($O^\epsilon$) $\le c\epsilon$, and

$$\|\nabla\Gamma(\mathcal{A}, \nabla\theta^{\epsilon})\|_{L^2(O^{\epsilon})} \le C \left( \|\nabla\theta^{\epsilon}\|_{L^{\infty}(O^{\epsilon})} \|\nabla\mathcal{A}\|_{L^2(O^{\epsilon})} + \|\nabla\nabla\theta^{\epsilon}\|_{L^{\infty}(O^{\epsilon})} \|\mathcal{A}\|_{L^2(O^{\epsilon})} \right) \\ \le C \left( \frac{1}{\epsilon} \|\nabla\mathcal{A}\|_{L^2(O^{\epsilon})} + \frac{1}{\epsilon^2} \|\mathcal{A}\|_{L^2(O^{\epsilon})} \right) \\ \le C \left( \frac{1}{\sqrt{\epsilon}} \|\nabla\mathcal{A}\|_{L^{\infty}(\Omega)} + \frac{1}{\epsilon\sqrt{\epsilon}} \|\mathcal{A}\|_{L^{\infty}(O^{\epsilon})} \right).$$

We observe now that any point $x \in O^\epsilon$ can be connected to a point $x' \in \partial\Omega$ by a segment of length not greater than $c\epsilon$ lying in $O^\epsilon$. Recalling that $\mathcal{A}(x') = 0$ and using a Taylor expansion at order 1 with integral remainder gives $|\mathcal{A}|(x) \le c\epsilon \|\nabla\mathcal{A}(z)\|_{L^\infty(O^\epsilon)}$. Consequently, $\|\mathcal{A}\|_{L^\infty(O^\epsilon)} \le c\epsilon \|\nabla\mathcal{A}\|_{L^\infty(O^\epsilon)}$ and

$$\|\nabla\Gamma(\mathcal{A}, \nabla\theta^{\epsilon})\|_{L^2(O^{\epsilon})} \le \frac{C}{\sqrt{\epsilon}} \|\nabla\mathcal{A}\|_{L^{\infty}(O^{\epsilon})},$$

which concludes the proof. $\square$

## 4.4 Representation of divergence-free fields

A key ingredient in the proof of Theorem 4.1 is a representation of periodic divergence-free and zero-mean fields as the image of a differential operator, the range of which is always divergence-free. This representation is given by the following.

**Lemma 4.15.** Let $g \in L^2_\#(Y)^d$ be periodic, satisfying

$$\int_Y g dy = 0, \quad \mathrm{div}_y g = 0 \text{ in } Y.$$

There exists a (non-unique) periodic solution $\varphi \in H^1_\#(Y)^d \cap H^2_\#(Y)^d$ to the problem

$$\begin{cases} -\Delta_y \varphi + \nabla_y \mathrm{div}_y \varphi = g & \text{in } Y, \\ y \to \varphi(y) & \text{Y-periodic,} \end{cases} \qquad (25)$$
Proof. Let $g \in C^{1,\alpha}(\bar{\Omega})^d$ with $\operatorname{div} g = 0$ in $\Omega$, $g \cdot n = 0$ on $\partial\Omega$. By virtue of Lemma 4.11, there exists a function $A \in C^{2,\alpha}(\bar{\Omega})^{d(d-1)/2}$, such that $\operatorname{curl}^*(A) = g$ on $\partial\Omega$. Define:

$$v_{b\epsilon} = \operatorname{curl}^*(\theta^{\epsilon} A),$$

where $\theta^{\epsilon}$ is the cut-off function introduced in Definition 4.12. Obviously, $\operatorname{div} v_{b\epsilon} = 0$ in $\Omega$. Then, $v_{b\epsilon} = g$ on $\partial\Omega$ since, in the vicinity $O^{\epsilon/2}$ of $\partial\Omega$, we recall that $\theta^{\epsilon} \equiv 1$. Finally, it remains to bound $|v_{b\epsilon}|_{H^1(\Omega)}$. We write

$$\operatorname{curl}^*(\theta^{\epsilon} A) = \theta^{\epsilon}\operatorname{curl}^*(A) + \Gamma(A, \nabla\theta^{\epsilon}),$$

where $\Gamma(A, \nabla\theta^{\epsilon})$ is defined as:

$$\Gamma(A, \nabla\theta^{\epsilon})_i = \sum_{j=1}^{i-1} A_{ji} \frac{\partial\theta^{\epsilon}}{\partial x_j} - \sum_{j=i+1}^{d} A_{ij} \frac{\partial\theta^{\epsilon}}{\partial x_j}.$$

It follows therefore:

$$|v_{b\epsilon}|_{H^1(\Omega)} \le \|\theta^{\epsilon} \nabla \operatorname{curl}^*(A)\|_{L^2(O^{\epsilon})} + \|\nabla\theta^{\epsilon} \operatorname{curl}^*(A)\|_{L^2(O^{\epsilon})} + \|\nabla\Gamma(A, \nabla\theta^{\epsilon})\|_{L^2(O^{\epsilon})} \\ \le C \left( \sqrt{\epsilon} \|\nabla\nabla A\|_{L^{\infty}(O^{\epsilon})} + \frac{1}{\sqrt{\epsilon}} \|\nabla A\|_{L^{\infty}(O^{\epsilon})} \right) + \|\nabla\Gamma(A, \nabla\theta^{\epsilon})\|_{L^2(O^{\epsilon})},$$

since meas($O^{\epsilon}$) $\le c\epsilon$, and

$$\|\nabla\Gamma(\mathbf{A}, \nabla\theta^{\epsilon})\|_{L^2(O^{\epsilon})} \le C \left( \|\nabla\theta^{\epsilon}\|_{L^{\infty}(O^{\epsilon})} \|\nabla\mathbf{A}\|_{L^2(O^{\epsilon})} + \|\nabla\nabla\theta^{\epsilon}\|_{L^{\infty}(O^{\epsilon})} \|\mathbf{A}\|_{L^2(O^{\epsilon})} \right) \\ \le C \left( \frac{1}{\epsilon} \|\nabla\mathbf{A}\|_{L^2(O^{\epsilon})} + \frac{1}{\epsilon^2} \|\mathbf{A}\|_{L^2(O^{\epsilon})} \right) \\ \le C \left( \frac{1}{\sqrt{\epsilon}} \|\nabla\mathbf{A}\|_{L^{\infty}(\Omega)} + \frac{1}{\epsilon\sqrt{\epsilon}} \|\mathbf{A}\|_{L^{\infty}(O^{\epsilon})} \right).$$

We observe now that any point $x \in O^{\epsilon}$ can be connected to a point $x' \in \partial\Omega$ by a segment of length not greater than $c\epsilon$ lying in $O^{\epsilon}$. Recalling that $A(x') = 0$ and using a Taylor expansion at order 1 with integral remainder gives $|A|(x) \le c\epsilon \|\nabla A(z)\|_{L^{\infty}(O^{\epsilon})}$. Consequently, $\|A\|_{L^{\infty}(O^{\epsilon})} \le c\epsilon \|\nabla A\|_{L^{\infty}(O^{\epsilon})}$ and

$$\|\nabla\Gamma(\mathbf{A}, \nabla\theta^{\epsilon})\|_{L^2(O^{\epsilon})} \le \frac{C}{\sqrt{\epsilon}} \|\nabla\mathbf{A}\|_{L^{\infty}(O^{\epsilon})},$$

which concludes the proof. $\square$

## 4.4 Representation of divergence-free fields

A key ingredient in the proof of Theorem 4.1 is a representation of periodic divergence-free and zero-mean fields as the image of a differential operator, the range of which is always divergence-free. This representation is given by the following.

**Lemma 4.15.** Let $g \in L^2_\#(Y)^d$ be periodic, satisfying

$$\int_Y g dy = 0, \quad \operatorname{div}_y g = 0 \text{ in } Y.$$

There exists a (non-unique) periodic solution $\varphi \in H^1_\#(Y)^d \cap H^2_\#(Y)^d$ to the problem

$$\begin{cases} -\Delta_y \varphi + \nabla_y \operatorname{div}_y \varphi = g & \text{in } Y, \\ y \to \varphi(y) & \text{Y-periodic,} \end{cases} \qquad (25)$$
which satisfies, for some constant $C > 0$,

$$ \|\varphi\|_{H_{\#}^2(Y)} \leq C \|g\|_{L_{\#}^2(Y)}. $$

Furthermore, for any $\psi \in H_{\#}^2(Y)^d$, it holds that

$$ \operatorname{div}_y (-\Delta_y \psi + \nabla_y \operatorname{div}_y \psi) = 0 \text{ in } Y. $$

**Remark 4.16.** More generally, the identity

$$ \operatorname{div} (-\Delta\psi + \nabla \operatorname{div} \psi) = 0 \text{ in } \mathbb{R}^d, $$

holds true for any vector field $\psi \in (H^2(\mathbb{R}^d))^d$.

*Proof.* Using Fourier series the periodic solution $\varphi$ and the periodic right-hand side $g$ are written as

$$ \varphi(y) = \sum_{\xi \in \mathbb{Z}^d} \hat{\varphi}(\xi) e^{2i\pi \xi \cdot y}, \quad g(y) = \sum_{\xi \in \mathbb{Z}^d} \hat{g}(\xi) e^{2i\pi \xi \cdot y}, $$

where $\hat{g}(0) = 0$ and $\hat{g}(\xi) \cdot \xi = 0$ for any $\xi \in \mathbb{Z}^d$, because of the assumptions on $g$. Then, (25) becomes:

$$ 4\pi^2 |\xi|^2 \hat{\varphi} - 4\pi^2 (\xi \cdot \hat{\varphi}) \xi = \hat{g}. \qquad (26) $$

By writing (26) under matrix form, we have:

$$ A\hat{\varphi} = \hat{g}, \quad A(\xi) = 4\pi^2 (|\xi|^2 I_d - \xi \otimes \xi). $$

We easily see that $\ker(A) = \operatorname{span}(\xi)$. Since $A$ is symmetric, we have $\operatorname{Im}(A) = \ker(A)^\perp = \{u \in \mathbb{R}^d \text{ s.t. } u \cdot \xi = 0\}$ and $\hat{g} \in \ker(A)^\perp$. Now, if we restrict the matrix $A$ to its image, it is reduced to $|\xi|^2 I_d$. Consequently, its inverse is simply given by $|\xi|^{-2} I_d$ which gives the following solution:

$$ \hat{\varphi} = \frac{\hat{g}}{4\pi^2 |\xi|^2} \text{ for any } \xi \neq 0, $$

which is unique, up to the addition of a multiple of $\xi$. In other words, we obtained a solution $\varphi(y)$ of (25) which is unique up to the addition of a constant and of a periodic gradient. Now, recalling that $\xi$ is an integer vector, we note that:

$$ |\hat{\varphi}| = \frac{|\hat{g}|}{4\pi^2 |\xi|^2} \leq C|\hat{g}|, \quad |\xi\hat{\varphi}| \leq C|\hat{g}|, \quad |\xi\xi\hat{\varphi}| \leq C|\hat{g}|, $$

which allows to conclude that $\varphi$ belongs to $H_{\#}^2(Y)^d$. $\square$

In view of the proof of Theorem 4.1, Lemma 4.15 is applied to the solution $\omega_i$ of the cell problem (6). More precisely, introduce

$$ \omega_i' = \omega_i - |Y_F| \langle \omega_i \rangle = \omega_i - \int_Y \omega_i, $$

which satisfies

$$ \left\{ \begin{array}{l} \operatorname{div}_y \omega_i' = 0 \text{ on } Y, \\ \omega_i' \text{ is } \mathbb{Z}^d\text{-periodic,} \\ \omega_i' \text{ has zero mean over } Y. \end{array} \right. $$

According to Lemma 4.15, there exists a periodic function $\psi_i \in H_{\#}^2(Y)^d$ such that:

$$ -\Delta_y \psi_i + \nabla_y \operatorname{div}_y \psi_i = \omega_i' \text{ in } Y. \qquad (27) $$

In particular, since $\omega_i \in H_{\#}^1(Y)^d$, then $\psi_i \in H_{\#}^3(Y)^d$. This regularity cannot be improved since, although $\omega_i$ is smooth inside $Y_F$, it is not in the full cell $Y$. In particular, it implies that, at least for large dimension $d$, $\psi_i$ and $\nabla_y \psi_i$ may not belong to $L^\infty(Y)$. This lack of regularity requires special care in the next subsection.
which satisfies, for some constant $C > 0$,

$$ \|\varphi\|_{H_{\#}^2(Y)} \leq C \|g\|_{L_{\#}^2(Y)}. $$

Furthermore, for any $\psi \in H_{\#}^2(Y)^d$, it holds that

$$ \operatorname{div}_y (-\Delta_y \psi + \nabla_y \operatorname{div}_y \psi) = 0 \text{ in } Y. $$

**Remark 4.16.** More generally, the identity

$$ \operatorname{div} (-\Delta\psi + \nabla\operatorname{div}\psi) = 0 \text{ in } \mathbb{R}^d, $$

holds true for any vector field $\psi \in (H^2(\mathbb{R}^d))^d$.

*Proof.* Using Fourier series the periodic solution $\varphi$ and the periodic right-hand side $g$ are written as

$$ \varphi(y) = \sum_{\xi \in \mathbb{Z}^d} \hat{\varphi}(\xi) e^{2i\pi \xi \cdot y}, \quad g(y) = \sum_{\xi \in \mathbb{Z}^d} \hat{g}(\xi) e^{2i\pi \xi \cdot y}, $$

where $\hat{g}(0) = 0$ and $\hat{g}(\xi) \cdot \xi = 0$ for any $\xi \in \mathbb{Z}^d$, because of the assumptions on $g$. Then, (25) becomes:

$$ 4\pi^2 |\xi|^2 \hat{\varphi} - 4\pi^2 (\xi \cdot \hat{\varphi}) \xi = \hat{g}. \quad (26) $$

By writing (26) under matrix form, we have:

$$ A\hat{\varphi} = \hat{g}, \quad A(\xi) = 4\pi^2 (|\xi|^2 I_d - \xi \otimes \xi). $$

We easily see that $\ker(A) = \operatorname{span}(\xi)$. Since $A$ is symmetric, we have $\operatorname{Im}(A) = \ker(A)^\perp = \{u \in \mathbb{R}^d \text{ s.t. } u \cdot \xi = 0\}$ and $\hat{g} \in \ker(A)^\perp$. Now, if we restrict the matrix $A$ to its image, it is reduced to $|\xi|^2 I_d$. Consequently, its inverse is simply given by $|\xi|^{-2} I_d$ which gives the following solution:

$$ \hat{\varphi} = \frac{\hat{g}}{4\pi^2|\xi|^2} \text{ for any } \xi \neq 0, $$

which is unique, up to the addition of a multiple of $\xi$. In other words, we obtained a solution $\varphi(y)$ of (25) which is unique up to the addition of a constant and of a periodic gradient. Now, recalling that $\xi$ is an integer vector, we note that:

$$ |\hat{\varphi}| = \frac{|\hat{g}|}{4\pi^2|\xi|^2} \leq C|\hat{g}|, \quad |\xi\hat{\varphi}| \leq C|\hat{g}|, \quad |\xi\xi\hat{\varphi}| \leq C|\hat{g}|, $$

which allows to conclude that $\varphi$ belongs to $H_{\#}^2(Y)^d$. $\square$

In view of the proof of Theorem 4.1, Lemma 4.15 is applied to the solution $\omega_i$ of the cell problem (6). More precisely, introduce

$$ \omega_i' = \omega_i - |Y_F| \langle \omega_i \rangle = \omega_i - \int_Y \omega_i, $$

which satisfies

$$ \left\{ \begin{array}{ll} \operatorname{div}_y \omega_i' = 0 & \text{on } Y, \\ \omega_i' \text{ is } \mathbb{Z}^d\text{-periodic,} \\ \omega_i' \text{ has zero mean over } Y. \end{array} \right. $$

According to Lemma 4.15, there exists a periodic function $\psi_i \in H_{\#}^2(Y)^d$ such that:

$$ -\Delta_y \psi_i + \nabla_y \operatorname{div}_y \psi_i = \omega_i' \text{ in } Y. \quad (27) $$

In particular, since $\omega_i \in H_{\#}^1(Y)^d$, then $\psi_i \in H_{\#}^3(Y)^d$. This regularity cannot be improved since, although $\omega_i$ is smooth inside $Y_F$, it is not in the full cell $Y$. In particular, it implies that, at least for large dimension $d$, $\psi_i$ and $\nabla_y \psi_i$ may not belong to $L^\infty(Y)$. This lack of regularity requires special care in the next subsection.
## 4.5 Proof of Theorem 4.1

Inspired by the approximation (13) of the two-scale asymptotic expansion, namely $u_{\varepsilon} \approx \varepsilon^2 u_2 + \varepsilon^3 u_3$, the main idea is to compare the velocity $u_{\varepsilon}$ with $\varepsilon^2 u_2 + \varepsilon^3 u_3$. The difficulty is that the velocity correctors $u_2, u_3$, built in Section 3, do not satisfy the boundary condition, i.e. $u_2, u_3$ do not vanish on $\partial\Omega$, and are not divergence-free. To satisfy this Dirichlet boundary condition for an approximation of $u_{\varepsilon}$, we introduce a cut-off function $\eta^{\varepsilon}$ as in [11].

**Definition 4.17.** For a small enough constant $c > 0$, define a neighborhood of the boundary $\partial\Omega$, denoted by $O^{\varepsilon} := \{x \in \Omega : \text{dist}(x, \partial\Omega) \le c\varepsilon\}$ (band of width $c\varepsilon$ near $\partial\Omega$). Then we choose a cut-off function $\eta^{\varepsilon} \in C^{\infty}(\bar{\Omega})$ such that $\eta^{\varepsilon}$ and all its derivatives vanish on $\partial\Omega$, while $\eta^{\varepsilon}(x) = 1$ on $\Omega \setminus O^{\varepsilon}$ and

$$ \|1 - \eta^{\varepsilon}\|_{L^2(\Omega)} \le C\sqrt{\varepsilon}, \quad \|\eta^{\varepsilon}\|_{L^{\infty}(\Omega)} = 1, \quad \|\nabla\eta^{\varepsilon}\|_{L^{\infty}(\Omega)} \le \frac{C}{\varepsilon}, \quad \|\nabla\nabla\eta^{\varepsilon}\|_{L^{\infty}(\Omega)} \le \frac{C}{\varepsilon^2}. $$

If we simply multiply the approximation $\varepsilon^2 u_2 + \varepsilon^3 u_3$ by the cut-off function $\eta^{\varepsilon}$, of course the Dirichlet boundary condition will be satisfied on $\partial\Omega$ but the divergence of the product is not under control. Therefore, a more clever use of $\eta^{\varepsilon}$ is required. It turns out that only $u_2$ needs a special treatment. Recall from (7) that

$$ u_2(x, y) = \frac{1}{\nu}\omega_i(y)\left(f_i(x) - \frac{\partial p^*}{\partial x_i}(x)\right). $$

The main idea is to decompose $\omega_i$ as

$$ \omega_i = \omega_i' + |Y_F|\langle\omega_i\rangle, $$

and to use (27) to replace $\omega_i'$ by its representation as derivatives of the periodic function $\psi_i$. Then, the cut-off function is inserted inside this differential representation as

$$ (-\Delta + \nabla \text{div}) (\eta^{\varepsilon} \psi_i), $$

which has the advantage of being divergence-free and vanishing on the boundary $\partial\Omega$. We shall not apply the cut-off function to the other (non-oscillating) term

$$ \frac{1}{\nu}|Y_F|\langle\omega_i\rangle\left(f_i(x) - \frac{\partial p^*}{\partial x_i}(x)\right) \equiv |Y_F|u(x), $$

where $u$ is divergence-free (because it is the Darcy velocity) and satisfies $u \cdot n = 0$ on $\partial\Omega$. To correct this boundary condition we apply Lemma 4.13 to $g = |Y_F|u$. The assumptions of Lemma 4.13 are satisfied, by its definition (9), the Darcy velocity $u$ belongs to $C^{1,\alpha}(\bar{\Omega})^d$, satisfies $\text{div} u = 0$ in $\Omega$ and $u \cdot n = 0$ on $\partial\Omega$. Thus Lemma 4.13 ensures that there exists $v_{b\varepsilon} \in H^1(\Omega)^d$ supported in $O^{\varepsilon}$ such that:

$$ \begin{gathered} \text{div } v_{b\varepsilon} = 0 \text{ on } \Omega, \qquad v_{b\varepsilon} = |Y_F|u \text{ on } \partial\Omega, \\ |v_{b\varepsilon}|_{H^1(\Omega^{\varepsilon})} \le \frac{C}{\sqrt{\varepsilon}} \|u\|_{C^{1,\alpha}(\bar{\Omega})} \end{gathered} \qquad (28) $$

The key observation is that, because of our construction of $\Omega^{\varepsilon}$, where obstacles are removed from $\Omega$ only if the entire cell is included in $\Omega$, all solid obstacles are uniformly away of the boundary $\partial\Omega$ by a distance of the order of $\varepsilon$. Therefore, in the choice of the cut-off function $\eta^{\varepsilon}$ (see Definition 4.17), one can choose the constant $c$ small enough so that no obstacles are contained or intersect the neighborhood $O^{\varepsilon}$ of the boundary $\partial\Omega$, where the support of the derivatives of $\eta^{\varepsilon}$ is restricted.
## 4.5 Proof of Theorem 4.1

Inspired by the approximation (13) of the two-scale asymptotic expansion, namely $u_{\varepsilon} \approx \varepsilon^2 u_2 + \varepsilon^3 u_3$, the main idea is to compare the velocity $u_{\varepsilon}$ with $\varepsilon^2 u_2 + \varepsilon^3 u_3$. The difficulty is that the velocity correctors $u_2, u_3$, built in Section 3, do not satisfy the boundary condition, i.e. $u_2, u_3$ do not vanish on $\partial\Omega$, and are not divergence-free. To satisfy this Dirichlet boundary condition for an approximation of $u_{\varepsilon}$, we introduce a cut-off function $\eta^{\varepsilon}$ as in [11].

**Definition 4.17.** For a small enough constant $c > 0$, define a neighborhood of the boundary $\partial\Omega$, denoted by $O^{\varepsilon} := \{x \in \Omega : \text{dist}(x, \partial\Omega) \le c\varepsilon\}$ (band of width $c\varepsilon$ near $\partial\Omega$). Then we choose a cut-off function $\eta^{\varepsilon} \in C^{\infty}(\bar{\Omega})$ such that $\eta^{\varepsilon}$ and all its derivatives vanish on $\partial\Omega$, while $\eta^{\varepsilon}(x) = 1$ on $\Omega \setminus O^{\varepsilon}$ and

$$ \|1 - \eta^{\varepsilon}\|_{L^2(\Omega)} \le C\sqrt{\varepsilon}, \quad \|\eta^{\varepsilon}\|_{L^{\infty}(\Omega)} = 1, \quad \|\nabla\eta^{\varepsilon}\|_{L^{\infty}(\Omega)} \le \frac{C}{\varepsilon}, \quad \|\nabla\nabla\eta^{\varepsilon}\|_{L^{\infty}(\Omega)} \le \frac{C}{\varepsilon^2}. $$

If we simply multiply the approximation $\varepsilon^2 u_2 + \varepsilon^3 u_3$ by the cut-off function $\eta^{\varepsilon}$, of course the Dirichlet boundary condition will be satisfied on $\partial\Omega$ but the divergence of the product is not under control. Therefore, a more clever use of $\eta^{\varepsilon}$ is required. It turns out that only $u_2$ needs a special treatment. Recall from (7) that

$$ u_2(x, y) = \frac{1}{\nu}\omega_i(y)\left(f_i(x) - \frac{\partial p^*}{\partial x_i}(x)\right). $$

The main idea is to decompose $\omega_i$ as

$$ \omega_i = \omega_i' + |Y_F|\langle\omega_i\rangle, $$

and to use (27) to replace $\omega_i'$ by its representation as derivatives of the periodic function $\psi_i$. Then, the cut-off function is inserted inside this differential representation as

$$ (-\Delta + \nabla \text{div}) (\eta^{\varepsilon} \psi_i), $$

which has the advantage of being divergence-free and vanishing on the boundary $\partial\Omega$. We shall not apply the cut-off function to the other (non-oscillating) term

$$ \frac{1}{\nu}|Y_F|\langle\omega_i\rangle\left(f_i(x) - \frac{\partial p^*}{\partial x_i}(x)\right) \equiv |Y_F|\mathbf{u}(x), $$

where $\mathbf{u}$ is divergence-free (because it is the Darcy velocity) and satisfies $\mathbf{u} \cdot \mathbf{n} = 0$ on $\partial\Omega$. To correct this boundary condition we apply Lemma 4.13 to $g = |Y_F|\mathbf{u}$. The assumptions of Lemma 4.13 are satisfied, by its definition (9), the Darcy velocity $\mathbf{u}$ belongs to $C^{1,\alpha}(\bar{\Omega})^d$, satisfies $\text{div} \mathbf{u} = 0$ in $\Omega$ and $\mathbf{u} \cdot \mathbf{n} = 0$ on $\partial\Omega$. Thus Lemma 4.13 ensures that there exists $v_{b\varepsilon} \in H^1(\Omega)^d$ supported in $O^{\varepsilon}$ such that:

$$ \begin{gathered} \text{div } v_{b\varepsilon} = 0 \text{ on } \Omega, \qquad v_{b\varepsilon} = |Y_F|\mathbf{u} \text{ on } \partial\Omega, \\ |v_{b\varepsilon}|_{H^1(\Omega^{\varepsilon})} \le \frac{C}{\sqrt{\varepsilon}} \|\mathbf{u}\|_{C^{1,\alpha}(\bar{\Omega})} \end{gathered} \qquad (28) $$

The key observation is that, because of our construction of $\Omega^{\varepsilon}$, where obstacles are removed from $\Omega$ only if the entire cell is included in $\Omega$, all solid obstacles are uniformly away of the boundary $\partial\Omega$ by a distance of the order of $\varepsilon$. Therefore, in the choice of the cut-off function $\eta^{\varepsilon}$ (see Definition 4.17), one can choose the constant $c$ small enough so that no obstacles are contained or intersect the neighborhood $O^{\varepsilon}$ of the boundary $\partial\Omega$, where the support of the derivatives of $\eta^{\varepsilon}$ is restricted.
Based on the asymptotic expansion of order 3 for the velocity, $\epsilon^2\mathbf{u}_2 + \epsilon^3\mathbf{u}_3$, we define the following velocity estimator:

$$ \begin{aligned} \mathbf{u}_{\varepsilon,3}(x) = & \frac{\varepsilon^2}{\nu} |Y_F| \langle \omega_i \rangle \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) - \varepsilon^2 \mathbf{v}_{b\varepsilon} \\ & + \frac{\varepsilon^4}{\nu} (-\Delta + \nabla \text{div}) \left( \eta^\varepsilon(x) \psi_i \left( \frac{x}{\varepsilon} \right) \right) \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) \\ & + \frac{\varepsilon^3}{\nu} \gamma_{ij} \left( \frac{x}{\varepsilon} \right) \eta^\varepsilon(x) \frac{\partial}{\partial x_j} \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) - \mathbf{v}_{d\varepsilon}, \end{aligned} \quad (29) $$

where $\mathbf{v}_{b\varepsilon}$, defined in (28), corrects the boundary condition of the first term above and $\mathbf{v}_{d\varepsilon} \in H_0^1(\Omega^\varepsilon)^d$ is going to be defined in the proof of Lemma 4.18 to correct the divergence-free condition of $\mathbf{u}_{\varepsilon,3}$. Recall that $\psi_i$ is defined by (27), $\omega_i$ by (6), $\gamma_{ij}$ by (12) and all functions are extended by 0 inside the obstacles so that $\mathbf{u}_{\varepsilon,3}$ is well defined on the whole $\Omega$. Similarly, we define an estimator for the pressure, which is exactly the asymptotic expansion at first order of the pressure (no need of using a cut-off function), as:

$$ p_{\varepsilon,1} = p^* + \varepsilon p_1, \quad (30) $$

where $p^*$ is the homogenized pressure, solution of the Darcy problem (10) and $p_1$ is defined by (7).

**Lemma 4.18.** There exists $\mathbf{v}_{d\varepsilon} \in H_0^1(\Omega^\varepsilon)^d$ such that the velocity estimator $\mathbf{u}_{\varepsilon,3}$, defined by (29), satisfies:

$$ \mathbf{u}_{\varepsilon,3} = \varepsilon^2 \mathbf{u}_2 + \varepsilon^3 \mathbf{u}_3 - \mathbf{v}_{d\varepsilon} \quad \text{in } \Omega \setminus O^\varepsilon, \quad (31) $$

$$ \mathbf{u}_{\varepsilon,3} = \mathbf{0} \quad \text{on } \partial\Omega, \quad (32) $$

$$ \mathbf{u}_{\varepsilon,3} = \mathbf{0} \quad \text{on } \partial B^\varepsilon, \quad (33) $$

$$ \text{div } \mathbf{u}_{\varepsilon,3} = 0 \quad \text{in } \Omega, \quad (34) $$

and

$$ |\mathbf{v}_{d\varepsilon}|_{H^1(\Omega^\varepsilon)} \le C\varepsilon^{\frac{3}{2}} \|\mathbf{f} - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})}, \quad (35) $$

with $C$ independent of $\varepsilon$.

*Proof.* To check (31), simply recall that $\eta^\varepsilon \equiv 1$ in $\Omega \setminus O^\varepsilon$. Thus by virtue of (27) the second line of (29) simplifies and, combined with the first line, yields $\epsilon^2\mathbf{u}_2$, since $\mathbf{v}_{b\varepsilon}$ is supported in $O^\varepsilon$. Similarly, the third line of (29) coincides with $\epsilon^3\mathbf{u}_3$ in $\Omega \setminus O^\varepsilon$.

To check (32), we recognize that the first line of (29) vanishes on $\partial\Omega$ by construction of $\mathbf{v}_{b\varepsilon}$, while the second and third line vanish too because $\eta^\varepsilon$ and all its derivatives vanish on $\partial\Omega$, as well as $\mathbf{v}_{d\varepsilon}$ by definition.

By construction of the cut-off function and Assumption 2.1 on the solid obstacles which stay away from the boundary $\partial\Omega$, we have $\eta^\varepsilon \equiv 1$ in the vicinity of $B^\varepsilon$. Therefore $\mathbf{u}_{\varepsilon,3} = \epsilon^2\mathbf{u}_2 + \epsilon^3\mathbf{u}_3 - \mathbf{v}_{d\varepsilon}$ on $\partial B^\varepsilon$. Furthermore, since $\omega_i$ and $\gamma_{ij}$ vanish on $Y_S$, as well as $\mathbf{v}_{d\varepsilon}$ on $\partial B^\varepsilon$ by definition, we deduce that (33) holds true.

It remains to compute the divergence of the estimator $\mathbf{u}_{\varepsilon,3}$ and to prove (34) and (35). The first line of (29) is divergence-free by construction and recall that $\text{div}(-\Delta + \nabla \text{div}) = 0$. Therefore,

$$ \begin{aligned} \text{div } \mathbf{u}_{\varepsilon,3}(x) = & \frac{\varepsilon^4}{\nu} (-\Delta + \nabla \text{div}) \left( \eta^\varepsilon(x) \psi_i \left( \frac{x}{\varepsilon} \right) \right) \cdot \nabla (f_i - \partial_i p^*)(x) \\ & + \frac{\varepsilon^3}{\nu} \text{div} \left( \gamma_{ij} \left( \frac{x}{\varepsilon} \right) \eta^\varepsilon(x) \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x) \right) - \text{div} \mathbf{v}_{d\varepsilon}, \end{aligned} $$
Based on the asymptotic expansion of order 3 for the velocity, $\epsilon^2\mathbf{u}_2 + \epsilon^3\mathbf{u}_3$, we define the following velocity estimator:

$$ \begin{aligned} \mathbf{u}_{\varepsilon,3}(x) = & \frac{\varepsilon^2}{\nu} |Y_F| \langle \omega_i \rangle \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) - \varepsilon^2 \mathbf{v}_{b\varepsilon} \\ & + \frac{\varepsilon^4}{\nu} (-\Delta + \nabla \mathrm{div}) \left( \eta^\varepsilon(x) \psi_i \left( \frac{x}{\varepsilon} \right) \right) \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) \\ & + \frac{\varepsilon^3}{\nu} \gamma_{ij} \left( \frac{x}{\varepsilon} \right) \eta^\varepsilon(x) \frac{\partial}{\partial x_j} \left( f_i(x) - \frac{\partial p^*}{\partial x_i}(x) \right) - \mathbf{v}_{d\varepsilon}, \end{aligned} \quad (29) $$

where $\mathbf{v}_{b\varepsilon}$, defined in (28), corrects the boundary condition of the first term above and $\mathbf{v}_{d\varepsilon} \in H_0^1(\Omega^\varepsilon)^d$ is going to be defined in the proof of Lemma 4.18 to correct the divergence-free condition of $\mathbf{u}_{\varepsilon,3}$. Recall that $\psi_i$ is defined by (27), $\omega_i$ by (6), $\gamma_{ij}$ by (12) and all functions are extended by 0 inside the obstacles so that $\mathbf{u}_{\varepsilon,3}$ is well defined on the whole $\Omega$. Similarly, we define an estimator for the pressure, which is exactly the asymptotic expansion at first order of the pressure (no need of using a cut-off function), as:

$$ p_{\varepsilon,1} = p^* + \varepsilon p_1, \quad (30) $$

where $p^*$ is the homogenized pressure, solution of the Darcy problem (10) and $p_1$ is defined by (7).

**Lemma 4.18.** There exists $\mathbf{v}_{d\varepsilon} \in H_0^1(\Omega^\varepsilon)^d$ such that the velocity estimator $\mathbf{u}_{\varepsilon,3}$, defined by (29), satisfies:

$$ \mathbf{u}_{\varepsilon,3} = \varepsilon^2 \mathbf{u}_2 + \varepsilon^3 \mathbf{u}_3 - \mathbf{v}_{d\varepsilon} \quad \text{in } \Omega \setminus O^\varepsilon, \quad (31) $$

$$ \mathbf{u}_{\varepsilon,3} = \mathbf{0} \quad \text{on } \partial\Omega, \quad (32) $$

$$ \mathbf{u}_{\varepsilon,3} = \mathbf{0} \quad \text{on } \partial B^\varepsilon, \quad (33) $$

$$ \mathrm{div} \, \mathbf{u}_{\varepsilon,3} = 0 \quad \text{in } \Omega, \quad (34) $$

and

$$ |\mathbf{v}_{d\varepsilon}|_{H^1(\Omega^\varepsilon)} \leq C \varepsilon^{\frac{3}{2}} \|\mathbf{f} - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})}, \quad (35) $$

with $C$ independent of $\varepsilon$.

*Proof.* To check (31), simply recall that $\eta^\varepsilon \equiv 1$ in $\Omega \setminus O^\varepsilon$. Thus by virtue of (27) the second line of (29) simplifies and, combined with the first line, yields $\epsilon^2\mathbf{u}_2$, since $\mathbf{v}_{b\varepsilon}$ is supported in $O^\varepsilon$. Similarly, the third line of (29) coincides with $\epsilon^3\mathbf{u}_3$ in $\Omega \setminus O^\varepsilon$.

To check (32), we recognize that the first line of (29) vanishes on $\partial\Omega$ by construction of $\mathbf{v}_{b\varepsilon}$, while the second and third line vanish too because $\eta^\varepsilon$ and all its derivatives vanish on $\partial\Omega$, as well as $\mathbf{v}_{d\varepsilon}$ by definition.

By construction of the cut-off function and Assumption 2.1 on the solid obstacles which stay away from the boundary $\partial\Omega$, we have $\eta^\varepsilon \equiv 1$ in the vicinity of $B^\varepsilon$. Therefore $\mathbf{u}_{\varepsilon,3} = \epsilon^2\mathbf{u}_2 + \epsilon^3\mathbf{u}_3 - \mathbf{v}_{d\varepsilon}$ on $\partial B^\varepsilon$. Furthermore, since $\omega_i$ and $\gamma_{ij}$ vanish on $Y_S$, as well as $\mathbf{v}_{d\varepsilon}$ on $\partial B^\varepsilon$ by definition, we deduce that (33) holds true.

It remains to compute the divergence of the estimator $\mathbf{u}_{\varepsilon,3}$ and to prove (34) and (35). The first line of (29) is divergence-free by construction and recall that $\mathrm{div}(-\Delta + \nabla \mathrm{div}) = 0$. Therefore,

$$ \mathrm{div} \, \mathbf{u}_{\varepsilon,3}(x) = & \frac{\varepsilon^4}{\nu} (-\Delta + \nabla \mathrm{div}) \left( \eta^\varepsilon(x) \psi_i \left( \frac{x}{\varepsilon} \right) \right) \cdot \nabla (f_i - \partial_i p^*)(x) \\ & + \frac{\varepsilon^3}{\nu} \mathrm{div} \left( \gamma_{ij} \left( \frac{x}{\varepsilon} \right) \eta^\varepsilon(x) \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x) \right) - \mathrm{div} \, \mathbf{v}_{d\varepsilon}, $$
where, for the sake of clarity, the partial derivative with respect to $x_i$ is simply denoted by $\partial_i$.
An easy computation shows that:

$$(-\Delta + \nabla \operatorname{div})(\psi_i \left(\frac{x}{\varepsilon}\right) \eta^\varepsilon(x)) = \frac{1}{\varepsilon^2} \eta^\varepsilon(x)(-\Delta_y + \nabla_y \operatorname{div}_y) \psi_i \left(\frac{x}{\varepsilon}\right) + \frac{2}{\varepsilon} (\nabla \eta^\varepsilon(x) \cdot \nabla_y \psi_i \left(\frac{x}{\varepsilon}\right) \\ - \operatorname{div}_y \psi_i \left(\frac{x}{\varepsilon}\right) \nabla \eta^\varepsilon(x)) + (-\psi_i \left(\frac{x}{\varepsilon}\right) \Delta \eta^\varepsilon(x) + \nabla \nabla \eta^\varepsilon(x) \psi_i \left(\frac{x}{\varepsilon}\right)),$$

which leads to

$$\operatorname{div} \mathbf{u}_{\varepsilon,3} = L_{\varepsilon} + G_{\varepsilon} - \operatorname{div} \mathbf{v}_{d\varepsilon},$$

where $L_\varepsilon$ is made of terms of order 2:

$$L_{\varepsilon} = \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x) (-\Delta_y + \nabla_y \operatorname{div}_y) \psi_i \left(\frac{x}{\varepsilon}\right) \cdot \nabla (f_i - \partial_i p^*)(x) \\ + \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x) \operatorname{div}_y \gamma_{ij} \left(\frac{x}{\varepsilon}\right) \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x),$$

and $G_\varepsilon$ consists of the remaining terms of order 3 and 4:

$$G_{\varepsilon} = \frac{\varepsilon^3}{\nu} \left( \nabla \eta^{\varepsilon}(x) \operatorname{div}_y \psi_i \left(\frac{x}{\varepsilon}\right) - \nabla_y \psi_i \left(\frac{x}{\varepsilon}\right) \nabla \eta^{\varepsilon}(x) \right) \cdot \nabla (f_i - \partial_i p^*)(x) \\ + \frac{\varepsilon^3}{\nu} \gamma_{ij} \left(\frac{x}{\varepsilon}\right) \cdot \nabla \eta^{\varepsilon}(x) \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x) + \frac{\varepsilon^3}{\nu} \gamma_{ij} \left(\frac{x}{\varepsilon}\right) \eta^{\varepsilon}(x) \cdot \nabla \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x) \quad (36) \\ + \frac{\varepsilon^4}{\nu} \left( - \psi_i \left(\frac{x}{\varepsilon}\right) \Delta \eta^{\varepsilon}(x) + \nabla \nabla \eta^{\varepsilon}(x) \psi_i \left(\frac{x}{\varepsilon}\right) \right) \cdot \nabla (f_i - \partial_i p^*)(x).$$

Recalling that $\operatorname{div}_y \gamma_{ij} = -\omega_i \cdot e_j + \langle \omega_i \cdot e_j \rangle$ and $(-\Delta_y + \nabla_y \operatorname{div}_y) \psi_i = \omega_i - |Y_F| \langle \omega_i \rangle$, one can check that $L_\varepsilon$ vanishes since:

$$L_{\varepsilon} = \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x)(\omega_i - |Y_F|\langle \omega_i \rangle) \cdot \nabla(f_i - \partial_i p^*)(x) - \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x)(\omega_i - \langle \omega_i \rangle) \cdot \nabla(f_i - \partial_i p^*)(x) \\ = \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x)(1 - |Y_F|)\langle \omega_i \rangle \cdot \nabla(f_i - \partial_i p^*)(x) = 0,$$

in view of the homogenized equation (10). Therefore, the estimator divergence simplifies as

$$\operatorname{div} \mathbf{u}_{\varepsilon,3} = G_{\varepsilon} - \operatorname{div} \mathbf{v}_{d\varepsilon}$$

and we now prove the following bound

$$\|G_{\varepsilon}\|_{L^2(\Omega^{\varepsilon})} \leq C \varepsilon^{\frac{5}{2}} \|\mathbf{f} - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (37)$$

Indeed, the third term of $G_\varepsilon$ in (36) is bounded by $C\varepsilon^3 \|(f_i - \partial_i p^*)\|_{H^2(\Omega)}$ since $\|\eta^\varepsilon\|_{L^\infty(\Omega)} = 1$ and $\gamma_{ij} \in H^{m+2}(Y_F)^d \subset L^\infty(Y_F)^d$ by Sobolev embeddings (note that this bound is slightly better than (37) and it is the only term in (36) to involve the $H^2$-norm). All other terms in (36) feature derivatives of $\eta^\varepsilon$ and have thus compact supports in $O^\varepsilon$. In all those terms we bound $(f_i - \partial_i p^*)$ by its $C^{1,\alpha}$-norm. Although $\gamma_{ij}$ is bounded in $L^\infty(Y_F)^d$ for the second term of (36), it is not clear that it is the case for $\psi_i$ and $\nabla_y \psi_i$, at least for any space dimension. Therefore, we rely on another argument that we explain on the first term of (36) (the other ones being treated similarly). Since $\|\nabla \eta^\varepsilon\|_{L^\infty(\Omega)} \leq C\varepsilon^{-1}$, this first term is bounded in $L^2(\Omega^\varepsilon)$ by

$$C\varepsilon^2 \|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})} \left( \int_{O^\varepsilon} |\nabla_y \psi_i \left(\frac{x}{\varepsilon}\right)|^2 dx \right)^{1/2} \\ \leq C\varepsilon^2 \|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})} \left( n_{O^\varepsilon} \varepsilon^d \|\nabla_y \psi_i(y)\|_{L^2(Y_F)}^2 \right)^{1/2}, \quad (38)$$
where, for the sake of clarity, the partial derivative with respect to $x_i$ is simply denoted by $\partial_i$. An easy computation shows that:

$$(-\Delta + \nabla \operatorname{div})(\psi_i \left(\frac{x}{\varepsilon}\right) \eta^\varepsilon(x)) = \frac{1}{\varepsilon^2} \eta^\varepsilon(x)(-\Delta_y + \nabla_y \operatorname{div}_y) \psi_i \left(\frac{x}{\varepsilon}\right) + \frac{2}{\varepsilon} (\nabla \eta^\varepsilon(x) \cdot \nabla_y \psi_i \left(\frac{x}{\varepsilon}\right) \\ - \operatorname{div}_y \psi_i \left(\frac{x}{\varepsilon}\right) \nabla \eta^\varepsilon(x)) + (-\psi_i \left(\frac{x}{\varepsilon}\right) \Delta \eta^\varepsilon(x) + \nabla \nabla \eta^\varepsilon(x) \psi_i \left(\frac{x}{\varepsilon}\right)),$$

which leads to

$$\operatorname{div} \mathbf{u}_{\varepsilon,3} = L_{\varepsilon} + G_{\varepsilon} - \operatorname{div} \mathbf{v}_{d\varepsilon},$$

where $L_\varepsilon$ is made of terms of order 2:

$$L_{\varepsilon} = \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x) (-\Delta_y + \nabla_y \operatorname{div}_y) \psi_i \left(\frac{x}{\varepsilon}\right) \cdot \nabla (f_i - \partial_i p^*)(x) \\ + \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x) \operatorname{div}_y \gamma_{ij} \left(\frac{x}{\varepsilon}\right) \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x),$$

and $G_\varepsilon$ consists of the remaining terms of order 3 and 4:

$$G_{\varepsilon} = \frac{\varepsilon^3}{\nu} \left( \nabla \eta^{\varepsilon}(x) \operatorname{div}_y \psi_i \left(\frac{x}{\varepsilon}\right) - \nabla_y \psi_i \left(\frac{x}{\varepsilon}\right) \nabla \eta^{\varepsilon}(x) \right) \cdot \nabla (f_i - \partial_i p^*)(x) \\ + \frac{\varepsilon^3}{\nu} \gamma_{ij} \left(\frac{x}{\varepsilon}\right) \cdot \nabla \eta^{\varepsilon}(x) \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x) + \frac{\varepsilon^3}{\nu} \gamma_{ij} \left(\frac{x}{\varepsilon}\right) \eta^{\varepsilon}(x) \cdot \nabla \frac{\partial}{\partial x_j} (f_i - \partial_i p^*)(x) \quad (36) \\ + \frac{\varepsilon^4}{\nu} \left( - \psi_i \left(\frac{x}{\varepsilon}\right) \Delta \eta^{\varepsilon}(x) + \nabla \nabla \eta^{\varepsilon}(x) \psi_i \left(\frac{x}{\varepsilon}\right) \right) \cdot \nabla (f_i - \partial_i p^*)(x).$$

Recalling that $\operatorname{div}_y \gamma_{ij} = -\omega_i \cdot e_j + \langle \omega_i \cdot e_j \rangle$ and $(-\Delta_y + \nabla_y \operatorname{div}_y) \psi_i = \omega_i - |Y_F| \langle \omega_i \rangle$, one can check that $L_\varepsilon$ vanishes since:

$$L_{\varepsilon} = \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x)(\omega_i - |Y_F|\langle \omega_i \rangle) \cdot \nabla(f_i - \partial_i p^*)(x) - \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x)(\omega_i - \langle \omega_i \rangle) \cdot \nabla(f_i - \partial_i p^*)(x) \\ = \frac{\varepsilon^2}{\nu} \eta^{\varepsilon}(x)(1 - |Y_F|)\langle \omega_i \rangle \cdot \nabla(f_i - \partial_i p^*)(x) = 0,$$

in view of the homogenized equation (10). Therefore, the estimator divergence simplifies as

$$\operatorname{div} \mathbf{u}_{\varepsilon,3} = G_{\varepsilon} - \operatorname{div} \mathbf{v}_{d\varepsilon}$$

and we now prove the following bound

$$\|G_{\varepsilon}\|_{L^2(\Omega^{\varepsilon})} \leq C \varepsilon^{\frac{5}{2}} \|\mathbf{f} - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (37)$$

Indeed, the third term of $G_\varepsilon$ in (36) is bounded by $C\varepsilon^3 \|(f_i - \partial_i p^*)\|_{H^2(\Omega)}$ since $\|\eta^\varepsilon\|_{L^\infty(\Omega)} = 1$ and $\gamma_{ij} \in H^{m+2}(Y_F)^d \subset L^\infty(Y_F)^d$ by Sobolev embeddings (note that this bound is slightly better than (37) and it is the only term in (36) to involve the $H^2$-norm). All other terms in (36) feature derivatives of $\eta^\varepsilon$ and have thus compact supports in $O^\varepsilon$. In all those terms we bound $(f_i - \partial_i p^*)$ by its $C^{1,\alpha}$-norm. Although $\gamma_{ij}$ is bounded in $L^\infty(Y_F)^d$ for the second term of (36), it is not clear that it is the case for $\psi_i$ and $\nabla_y \psi_i$, at least for any space dimension. Therefore, we rely on another argument that we explain on the first term of (36) (the other ones being treated similarly). Since $\|\nabla \eta^\varepsilon\|_{L^\infty(\Omega)} \leq C\varepsilon^{-1}$, this first term is bounded in $L^2(\Omega^\varepsilon)$ by

$$C\varepsilon^2 \|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})} \left( \int_{O^\varepsilon} |\nabla_y \psi_i \left(\frac{x}{\varepsilon}\right)|^2 dx \right)^{1/2} \\ \leq C\varepsilon^2 \|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})} \left( n_{O^\varepsilon} \varepsilon^d \|\nabla_y \psi_i(y)\|_{L^2(Y_F)}^2 \right)^{1/2}, \quad (38)$$
where we used the periodicity of $\nabla_y \psi_i$ and $n_{O^\epsilon} = O(\epsilon^{-d+1})$ is the number of cells $Y_i^\epsilon$ which cover $O^\epsilon$. This implies exactly the bound (37). Finally, remark that $G_\epsilon$ has zero mean on $\Omega^\epsilon$ since

$$ \int_{\Omega^\epsilon} \mathrm{div}(\mathbf{u}_{\epsilon,3} + \mathbf{v}_{d\epsilon}) = \int_{\partial\Omega^\epsilon} (\mathbf{u}_{\epsilon,3} + \mathbf{v}_{d\epsilon}) \cdot \mathbf{n} = 0, $$

because $\mathbf{u}_{\epsilon,3}$ and $\mathbf{v}_{d\epsilon}$ vanish on $\partial B^\epsilon$ and $\partial\Omega$, Consequently, to prove (34) we can rely on Lemma 4.9 which ensures that there exists $\mathbf{v}_{d\epsilon} \in H_0^1(\Omega^\epsilon)^d$ such that

$$ \mathrm{div} \, \mathbf{v}_{d\epsilon} = G_{\epsilon} \quad \text{and} \quad |\mathbf{v}_{d\epsilon}|_{H^1(\Omega^{\epsilon})} \leq \frac{C}{\epsilon} \|G_{\epsilon}\|_{L^2(\Omega^{\epsilon})} \leq C\epsilon^{\frac{3}{2}} \|\mathbf{f} - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})}, $$

which proves (34) and (35).

□

Now that all terms in the velocity estimator (29) are well defined, we turn to the estimates for the residuals of the homogenization process. Define the velocity residual $\mathbf{R}_{\mathbf{u}_\epsilon}$ and the pressure residual $R_{p_\epsilon}$ by:

$$ \mathbf{R}_{\mathbf{u}_\epsilon} = \mathbf{u}_\epsilon - \mathbf{u}_{\epsilon,3}, \quad R_{p_\epsilon} = p_\epsilon - p_{\epsilon,1}, $$

where the estimators $\mathbf{u}_{\epsilon,3}$ and $p_{\epsilon,1}$ are given respectively in (29) and (30).

**Lemma 4.19.** The velocity and pressure homogenization residuals verify:

$$ \left\{ \begin{array}{rcll} -\nu \Delta \mathbf{R}_{\mathbf{u}_\epsilon} + \nabla R_{p_\epsilon} & = & \mathbf{F}_\epsilon + \mathrm{div} \, \mathbf{H}_\epsilon & \text{in } \Omega^\epsilon, \\ \mathrm{div} \, \mathbf{R}_{\mathbf{u}_\epsilon} & = & 0 & \text{in } \Omega^\epsilon, \\ \mathbf{R}_{\mathbf{u}_\epsilon} & = & \mathbf{0} & \text{on } \partial B^\epsilon, \\ \mathbf{R}_{\mathbf{u}_\epsilon} & = & \mathbf{0} & \text{on } \partial \Omega, \end{array} \right. \quad (39) $$

where the source terms are bounded by:

$$ \| \mathbf{F}_{\epsilon} \|_{L^2(\Omega^{\epsilon})} \leq C_{\epsilon} \| \mathbf{f} - \nabla p^{*} \|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})}, \quad (40) $$

$$ \| \mathbf{H}_{\epsilon} \|_{L^2(\Omega^{\epsilon})} \leq C_{\epsilon}^{\frac{3}{2}} \| \mathbf{f} - \nabla p^{*} \|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})}, \quad (41) $$

with C independent of $\epsilon$.

*Proof.* Obviously, $\mathrm{div} \, \mathbf{R}_{\mathbf{u}_\epsilon} = 0$ in $\Omega^\epsilon$ and $\mathbf{R}_{\mathbf{u}_\epsilon} = \mathbf{0}$ on $\partial B^\epsilon \cup \partial\Omega$. We decompose the velocity residual as $\mathbf{R}_{\mathbf{u}_\epsilon} = (\mathbf{u}_\epsilon - \epsilon^2\mathbf{u}_2) + (\epsilon^2\mathbf{u}_2 - \mathbf{u}_{\epsilon,3})$ which leads to the first line of (39), where:

$$ \mathbf{F}_{\epsilon} = -\nu\Delta(\mathbf{u}_{\epsilon} - \epsilon^2\mathbf{u}_2) + \nabla(p_{\epsilon} - p_{\epsilon,1}) \quad \text{and} \quad \mathbf{H}_{\epsilon} = -\nu\nabla(\epsilon^2\mathbf{u}_2 - \mathbf{u}_{\epsilon,3}). $$

To prove estimate (40) for $\mathbf{F}_{\epsilon}$, we compute:

$$ \begin{aligned} \mathbf{F}_{\epsilon} = \mathbf{f} - \nabla p^* &+ (\Delta_y \omega_i - \nabla_y \pi_i) \left(\frac{x}{\epsilon}\right) (f_i - \partial_i p^*)(x) + 2\epsilon \nabla_y \omega_i \left(\frac{x}{\epsilon}\right) \nabla(f_i - \partial_i p^*)(x) \\ &+ \epsilon^2 \omega_i \left(\frac{x}{\epsilon}\right) \Delta(f_i - \partial_i p^*)(x) - \epsilon \pi_i \left(\frac{x}{\epsilon}\right) \nabla(f_i - \partial_i p^*)(x) \end{aligned} $$

and, rearranging the terms,

$$ \begin{aligned} \mathbf{F}_{\epsilon} = & (\Delta_y \omega_i - \nabla_y \pi_i + e_i) \left(\frac{x}{\epsilon}\right) (f_i - \partial_i p^*)(x) \\ & + 2\epsilon \nabla_y \omega_i \left(\frac{x}{\epsilon}\right) \nabla(f_i - \partial_i p^*)(x) - \epsilon \pi_i \left(\frac{x}{\epsilon}\right) \nabla(f_i - \partial_i p^*)(x) \\ & + \epsilon^2 \omega_i \left(\frac{x}{\epsilon}\right) \Delta(f_i - \partial_i p^*)(x). \end{aligned} \quad (42) $$

The first term in (42) actually vanishes because of the cell equation (6). The next three terms above are of order $\epsilon$ or higher. To bound the fourth term requires the $H^2$ regularity of $(f_i - \partial_i p^*)$,
where we used the periodicity of $\nabla_y \psi_i$ and $n_{O^\varepsilon} = O(\varepsilon^{-d+1})$ is the number of cells $Y_i^\varepsilon$ which cover $O^\varepsilon$. This implies exactly the bound (37). Finally, remark that $G_\varepsilon$ has zero mean on $\Omega^\varepsilon$ since

$$ \int_{\Omega^\varepsilon} \mathrm{div}(\mathbf{u}_{\varepsilon,3} + \mathbf{v}_{d\varepsilon}) = \int_{\partial\Omega^\varepsilon} (\mathbf{u}_{\varepsilon,3} + \mathbf{v}_{d\varepsilon}) \cdot \mathbf{n} = 0, $$

because $\mathbf{u}_{\varepsilon,3}$ and $\mathbf{v}_{d\varepsilon}$ vanish on $\partial B^\varepsilon$ and $\partial\Omega$, Consequently, to prove (34) we can rely on Lemma 4.9 which ensures that there exists $\mathbf{v}_{d\varepsilon} \in H_0^1(\Omega^\varepsilon)^d$ such that

$$ \mathrm{div} \, \mathbf{v}_{d\varepsilon} = G_\varepsilon \quad \text{and} \quad |\mathbf{v}_{d\varepsilon}|_{H^1(\Omega^\varepsilon)} \le \frac{C}{\varepsilon} \|G_\varepsilon\|_{L^2(\Omega^\varepsilon)} \le C\varepsilon^{\frac{3}{2}} \|\mathbf{f} - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})}, $$

which proves (34) and (35).

□

Now that all terms in the velocity estimator (29) are well defined, we turn to the estimates for the residuals of the homogenization process. Define the velocity residual $\mathbf{R}_{\mathbf{u}_\varepsilon}$ and the pressure residual $R_{p_\varepsilon}$ by:

$$ \mathbf{R}_{\mathbf{u}_\varepsilon} = \mathbf{u}_\varepsilon - \mathbf{u}_{\varepsilon,3}, \quad R_{p_\varepsilon} = p_\varepsilon - p_{\varepsilon,1}, $$

where the estimators $\mathbf{u}_{\varepsilon,3}$ and $p_{\varepsilon,1}$ are given respectively in (29) and (30).

**Lemma 4.19.** The velocity and pressure homogenization residuals verify:

$$ \left\{ \begin{array}{rcll} -\nu \Delta \mathbf{R}_{\mathbf{u}_\varepsilon} + \nabla R_{p_\varepsilon} & = & \mathbf{F}_\varepsilon + \mathrm{div} \, \mathbf{H}_\varepsilon & \text{in } \Omega^\varepsilon, \\ \mathrm{div} \, \mathbf{R}_{\mathbf{u}_\varepsilon} & = & 0 & \text{in } \Omega^\varepsilon, \\ \mathbf{R}_{\mathbf{u}_\varepsilon} & = & \mathbf{0} & \text{on } \partial B^\varepsilon, \\ \mathbf{R}_{\mathbf{u}_\varepsilon} & = & \mathbf{0} & \text{on } \partial \Omega, \end{array} \right. \qquad (39) $$

where the source terms are bounded by:

$$ \| \mathbf{F}_\varepsilon \|_{L^2(\Omega^\varepsilon)} \le C\varepsilon \| \mathbf{f} - \nabla p^* \|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})}, \qquad (40) $$

$$ \| \mathbf{H}_\varepsilon \|_{L^2(\Omega^\varepsilon)} \le C\varepsilon^{\frac{3}{2}} \| \mathbf{f} - \nabla p^* \|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})}, \qquad (41) $$

with C independent of $\varepsilon$.

*Proof.* Obviously, $\mathrm{div} \, \mathbf{R}_{\mathbf{u}_\varepsilon} = 0$ in $\Omega^\varepsilon$ and $\mathbf{R}_{\mathbf{u}_\varepsilon} = \mathbf{0}$ on $\partial B^\varepsilon \cup \partial\Omega$. We decompose the velocity residual as $\mathbf{R}_{\mathbf{u}_\varepsilon} = (\mathbf{u}_\varepsilon - \varepsilon^2\mathbf{u}_2) + (\varepsilon^2\mathbf{u}_2 - \mathbf{u}_{\varepsilon,3})$ which leads to the first line of (39), where:

$$ \mathbf{F}_\varepsilon = -\nu\Delta(\mathbf{u}_\varepsilon - \varepsilon^2\mathbf{u}_2) + \nabla(p_\varepsilon - p_{\varepsilon,1}) \quad \text{and} \quad \mathbf{H}_\varepsilon = -\nu\nabla(\varepsilon^2\mathbf{u}_2 - \mathbf{u}_{\varepsilon,3}). $$

To prove estimate (40) for $\mathbf{F}_\varepsilon$, we compute:

$$ \mathbf{F}_\varepsilon = \mathbf{f} - \nabla p^* + (\Delta_y \omega_i - \nabla_y \pi_i) \left(\frac{x}{\varepsilon}\right) (f_i - \partial_i p^*)(x) + 2\varepsilon \nabla_y \omega_i \left(\frac{x}{\varepsilon}\right) \nabla (f_i - \partial_i p^*)(x) \\ + \varepsilon^2 \omega_i \left(\frac{x}{\varepsilon}\right) \Delta (f_i - \partial_i p^*)(x) - \varepsilon \pi_i \left(\frac{x}{\varepsilon}\right) \nabla (f_i - \partial_i p^*)(x) $$

and, rearranging the terms,

$$ \begin{aligned} \mathbf{F}_\varepsilon = & (\Delta_y \omega_i - \nabla_y \pi_i + e_i) \left(\frac{x}{\varepsilon}\right) (f_i - \partial_i p^*)(x) \\ & + 2\varepsilon \nabla_y \omega_i \left(\frac{x}{\varepsilon}\right) \nabla (f_i - \partial_i p^*)(x) - \varepsilon \pi_i \left(\frac{x}{\varepsilon}\right) \nabla (f_i - \partial_i p^*)(x) \\ & + \varepsilon^2 \omega_i \left(\frac{x}{\varepsilon}\right) \Delta (f_i - \partial_i p^*)(x). \end{aligned} \qquad (42) $$

The first term in (42) actually vanishes because of the cell equation (6). The next three terms above are of order $\varepsilon$ or higher. To bound the fourth term requires the $H^2$ regularity of $(f_i - \partial_i p^*)$,
while its $C^{1,\alpha}$ regularity is sufficient for all other terms. In the end, we deduce the bound (40) from (42).

Next, to prove estimate (41) for $H_\varepsilon$, we rewrite:

$$ 
\begin{aligned}
\mathbf{u}_{\varepsilon,3} - \varepsilon^2 \mathbf{u}_2 = & -\frac{\varepsilon^2}{\nu} (1 - \eta^\varepsilon(x)) (|Y_F|\langle\omega_i\rangle - \omega_i\left(\frac{x}{\varepsilon}\right))(f_i - \partial_i p^*)(x) - \varepsilon^2 \mathbf{v}_{b\varepsilon} - \mathbf{v}_{d\varepsilon} \\
& + \frac{\varepsilon^3}{\nu} \eta^\varepsilon(x) \gamma_{ij}\left(\frac{x}{\varepsilon}\right) \partial_j (f_i - \partial_i p^*)(x) \\
& + \frac{\varepsilon^3}{\nu} \left(\nabla\eta^\varepsilon(x) \operatorname{div}_y \psi_i\left(\frac{x}{\varepsilon}\right) - \nabla_y \psi_i\left(\frac{x}{\varepsilon}\right) \nabla\eta^\varepsilon(x)\right) (f_i - \partial_i p^*)(x) \\
& + \frac{\varepsilon^4}{\nu} \left(-\psi_i\left(\frac{x}{\varepsilon}\right) \Delta\eta^\varepsilon(x) + \nabla\nabla\eta^\varepsilon(x)\psi_i\left(\frac{x}{\varepsilon}\right)\right) (f_i - \partial_i p^*)(x).
\end{aligned} 
\quad (43) $$

The gradient of (43), equal to $H_\varepsilon$ up to a $\nu$-factor, is:

$$ 
\begin{aligned}
\nabla(\mathbf{u}_{\varepsilon,3} - \varepsilon^2 \mathbf{u}_2) = & -\varepsilon^2 \nabla \mathbf{v}_{b\varepsilon} - \nabla \mathbf{v}_{d\varepsilon} - \frac{\varepsilon^2}{\nu} \nabla \eta^\varepsilon(|Y_F|\langle\omega_i\rangle - \omega_i)(f_i - \partial_i p^*) \\
& - \frac{\varepsilon}{\nu}(1 - \eta^\varepsilon)\nabla_y \omega_i(f_i - \partial_i p^*) + \frac{\varepsilon^2}{\nu}(\nabla\eta^\varepsilon\nabla_y \operatorname{div}_y \psi_i - \nabla_y\nabla_y\psi_i\nabla\eta^\varepsilon)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^3}{\nu}(\nabla\nabla\eta^\varepsilon \operatorname{div}_y \psi_i - \nabla_y\psi_i\nabla\nabla\eta^\varepsilon)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^3}{\nu}(-\nabla_y\psi_i\Delta\eta^\varepsilon + (\nabla\nabla\eta^\varepsilon)\nabla_y\psi_i)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^4}{\nu}(-\psi_i\nabla\Delta\eta^\varepsilon + \nabla(\nabla\nabla\eta^\varepsilon)\psi_i)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^2}{\nu}(1 - \eta^\varepsilon)(|Y_F|\langle\omega_i\rangle - \omega_i)\nabla(f_i - \partial_i p^*) + \frac{\varepsilon^3}{\nu}\nabla\eta^\varepsilon\gamma_{ij}\partial_j(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^3}{\nu}(\nabla\eta^\varepsilon \operatorname{div}_y \psi_i - \nabla_y\psi_i\nabla\eta^\varepsilon)\nabla(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^4}{\nu}(-\psi_i\Delta\eta^\varepsilon + (\nabla\nabla\eta^\varepsilon)\psi_i)\nabla(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^2}{\nu}\eta^\varepsilon\nabla_y\gamma_{ij}\partial_j(f_i - \partial_i p^*) + \frac{\varepsilon^3}{\nu}\eta^\varepsilon\gamma_{ij}\nabla\partial_j(f_i - \partial_i p^*).
\end{aligned} 
\quad (44) $$

We bound (44) in the same way as for the estimate of $G_\varepsilon$. The first two terms of (44) are bounded by (28) and (35). The last term of (44) is the only one requiring the $H^2$-regularity of $(f_i - \partial_i p^*)$: for all other terms the $C^{1,\alpha}$-regularity of $(f_i - \partial_i p^*)$ is sufficient. By Lemma 3.3 the cell solutions $\omega_i$, $\gamma_{ij}$ and the gradient $\nabla_y\omega_i$ belong to $L^\infty(Y)$. This is not the case for $\psi_i$ and its derivatives (at least for any space dimension $d$): therefore we rely on a similar trick as in (38) for all terms involving $\psi_i$. Using the bounds (4.17) on the cut-off function $\eta^\varepsilon$, it follows that the terms in the first five lines of (44) are of order $\varepsilon^{3/2}$. The terms of the three following lines are of order $\varepsilon^{5/2}$ and eventually in the last line the terms are respectively of order $\varepsilon^2$ and $\varepsilon^3$. We do not detail the bounds except for the seemingly higher order term

$$ 
\begin{aligned}
\|\frac{\varepsilon}{\nu}(1 - \eta^\varepsilon)\nabla_y\omega_i(f_i - \partial_i p^*)\|_{L^2(\Omega^\varepsilon)} &\le C\varepsilon\|1 - \eta^\varepsilon\|_{L^2(\Omega^\varepsilon)}\|\nabla_y\omega_i\|_{L^\infty(Y)}\|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})} \\
&\le C\varepsilon^{3/2}\|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})}
\end{aligned} 
$$

because of (4.17). In the end, we deduce estimate (41). □

**Remark 4.20.** We cannot have a better estimate due to the presence of boundary layers. Indeed, the lower order $\sqrt{\varepsilon}$ is due to the boundary condition $\mathbf{u}_\varepsilon = \mathbf{0}$ on $\partial\Omega$, which worsens the approximation near the boundary. Technically, this is taken into account by the introduction of the cut-off function $\eta^\varepsilon$ in the proof. This term satisfies the bound $|\eta^\varepsilon|_{H^1(\Omega)} \le \frac{C}{\sqrt{\varepsilon}}$, which is the origin of the $\sqrt{\varepsilon}$-order in the estimate.
while its $C^{1,\alpha}$ regularity is sufficient for all other terms. In the end, we deduce the bound (40) from (42).

Next, to prove estimate (41) for $H_\varepsilon$, we rewrite:

$$ 
\begin{aligned}
\mathbf{u}_{\varepsilon,3} - \varepsilon^2 \mathbf{u}_2 = & -\frac{\varepsilon^2}{\nu} (1 - \eta^\varepsilon(x)) (|Y_F|\langle\omega_i\rangle - \omega_i\left(\frac{x}{\varepsilon}\right))(f_i - \partial_i p^*)(x) - \varepsilon^2 \mathbf{v}_{b\varepsilon} - \mathbf{v}_{d\varepsilon} \\
& + \frac{\varepsilon^3}{\nu} \eta^\varepsilon(x) \gamma_{ij}\left(\frac{x}{\varepsilon}\right) \partial_j (f_i - \partial_i p^*)(x) \\
& + \frac{\varepsilon^3}{\nu} \left(\nabla\eta^\varepsilon(x) \operatorname{div}_y \psi_i\left(\frac{x}{\varepsilon}\right) - \nabla_y \psi_i\left(\frac{x}{\varepsilon}\right) \nabla\eta^\varepsilon(x)\right) (f_i - \partial_i p^*)(x) \\
& + \frac{\varepsilon^4}{\nu} \left(-\psi_i\left(\frac{x}{\varepsilon}\right) \Delta\eta^\varepsilon(x) + \nabla\nabla\eta^\varepsilon(x)\psi_i\left(\frac{x}{\varepsilon}\right)\right) (f_i - \partial_i p^*)(x).
\end{aligned} 
\quad (43) $$

The gradient of (43), equal to $H_\varepsilon$ up to a $\nu$-factor, is:

$$ 
\begin{aligned}
\nabla(\mathbf{u}_{\varepsilon,3} - \varepsilon^2 \mathbf{u}_2) = & -\varepsilon^2 \nabla \mathbf{v}_{b\varepsilon} - \nabla \mathbf{v}_{d\varepsilon} - \frac{\varepsilon^2}{\nu} \nabla \eta^\varepsilon(|Y_F|\langle\omega_i\rangle - \omega_i)(f_i - \partial_i p^*) \\
& - \frac{\varepsilon}{\nu}(1 - \eta^\varepsilon)\nabla_y \omega_i(f_i - \partial_i p^*) + \frac{\varepsilon^2}{\nu}(\nabla\eta^\varepsilon\nabla_y \operatorname{div}_y \psi_i - \nabla_y\nabla_y\psi_i\nabla\eta^\varepsilon)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^3}{\nu}(\nabla\nabla\eta^\varepsilon \operatorname{div}_y \psi_i - \nabla_y\psi_i\nabla\nabla\eta^\varepsilon)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^3}{\nu}(-\nabla_y\psi_i\Delta\eta^\varepsilon + (\nabla\nabla\eta^\varepsilon)\nabla_y\psi_i)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^4}{\nu}(-\psi_i\nabla\Delta\eta^\varepsilon + \nabla(\nabla\nabla\eta^\varepsilon)\psi_i)(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^2}{\nu}(1 - \eta^\varepsilon)(|Y_F|\langle\omega_i\rangle - \omega_i)\nabla(f_i - \partial_i p^*) + \frac{\varepsilon^3}{\nu}\nabla\eta^\varepsilon\gamma_{ij}\partial_j(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^3}{\nu}(\nabla\eta^\varepsilon \operatorname{div}_y \psi_i - \nabla_y\psi_i\nabla\eta^\varepsilon)\nabla(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^4}{\nu}(-\psi_i\Delta\eta^\varepsilon + (\nabla\nabla\eta^\varepsilon)\psi_i)\nabla(f_i - \partial_i p^*) \\
& + \frac{\varepsilon^2}{\nu}\eta^\varepsilon\nabla_y\gamma_{ij}\partial_j(f_i - \partial_i p^*) + \frac{\varepsilon^3}{\nu}\eta^\varepsilon\gamma_{ij}\nabla\partial_j(f_i - \partial_i p^*).
\end{aligned} 
\quad (44) $$

We bound (44) in the same way as for the estimate of $G_\varepsilon$. The first two terms of (44) are bounded by (28) and (35). The last term of (44) is the only one requiring the $H^2$-regularity of $(f_i - \partial_i p^*)$: for all other terms the $C^{1,\alpha}$-regularity of $(f_i - \partial_i p^*)$ is sufficient. By Lemma 3.3 the cell solutions $\omega_i$, $\gamma_{ij}$ and the gradient $\nabla_y\omega_i$ belong to $L^\infty(Y)$. This is not the case for $\psi_i$ and its derivatives (at least for any space dimension $d$): therefore we rely on a similar trick as in (38) for all terms involving $\psi_i$. Using the bounds (4.17) on the cut-off function $\eta^\varepsilon$, it follows that the terms in the first five lines of (44) are of order $\varepsilon^{3/2}$. The terms of the three following lines are of order $\varepsilon^{5/2}$ and eventually in the last line the terms are respectively of order $\varepsilon^2$ and $\varepsilon^3$. We do not detail the bounds except for the seemingly higher order term

$$ 
\begin{aligned}
\| \frac{\varepsilon}{\nu} (1 - \eta^\varepsilon) \nabla_y \omega_i (f_i - \partial_i p^*) \|_{L^2(\Omega^\varepsilon)} &\le C\varepsilon \|1 - \eta^\varepsilon\|_{L^2(\Omega^\varepsilon)} \| \nabla_y \omega_i \|_{L^\infty(Y)} \|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})} \\
&\le C\varepsilon^{3/2} \|(f_i - \partial_i p^*)\|_{C^{1,\alpha}(\bar{\Omega})}
\end{aligned} 
$$

because of (4.17). In the end, we deduce estimate (41). $\square$

**Remark 4.20.** We cannot have a better estimate due to the presence of boundary layers. Indeed, the lower order $\sqrt{\varepsilon}$ is due to the boundary condition $\mathbf{u}_\varepsilon = \mathbf{0}$ on $\partial\Omega$, which worsens the approximation near the boundary. Technically, this is taken into account by the introduction of the cut-off function $\eta^\varepsilon$ in the proof. This term satisfies the bound $|\eta^\varepsilon|_{H^1(\Omega)} \le \frac{C}{\sqrt{\varepsilon}}$, which is the origin of the $\sqrt{\varepsilon}$-order in the estimate.
Since the residual $R_{u_ε}$ is divergence-free and belongs to $H_0^1(\Omega^\epsilon)^d$, it can be used as test function, leading to

$$
\begin{align*}
\nu \int_{\Omega^\epsilon} \nabla R_{u_\epsilon} \cdot \nabla R_{u_\epsilon} &= \nu \int_{\Omega^\epsilon} \nabla R_{u_\epsilon} \cdot \nabla R_{u_\epsilon} - \int_{\Omega^\epsilon} R_{p_\epsilon} \operatorname{div} R_{u_\epsilon} \\
&= \int_{\Omega^\epsilon} F_\epsilon \cdot R_{u_\epsilon} - \int_{\Omega^\epsilon} H_\epsilon \cdot \nabla R_{u_\epsilon} \\
&\le \|F_\epsilon\|_{L^2(\Omega^\epsilon)} \|R_{u_\epsilon}\|_{L^2(\Omega^\epsilon)} + \|H_\epsilon\|_{L^2(\Omega^\epsilon)} \|R_{u_\epsilon}\|_{H^1(\Omega^\epsilon)}.
\end{align*}
$$

Using estimate (40) and the Poincaré inequality in a perforated domain (Lemma 4.6), it follows:

$$
\| \mathbf{F}_{\varepsilon} \|_{L^2(\Omega^{\varepsilon})} \| \mathbf{R}_{u_{\varepsilon}} \|_{L^2(\Omega^{\varepsilon})} \leq C \varepsilon^2 \| \mathbf{f} - \nabla p^* \|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})} |\mathbf{R}_{u_{\varepsilon}}|_{H^1(\Omega^{\varepsilon})}.
$$

Then, using (41), we conclude that:

$$
|R_{u_\epsilon}|_{H^1(\Omega^\epsilon)} \le C\epsilon^{\frac{3}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (45)
$$

Now, writing $u_\epsilon - \epsilon^2 u_2 = R_{u_\epsilon} + u_{\epsilon,3} - \epsilon^2 u_2$, using (45) and (41) yields:

$$
|u_{\epsilon} - \epsilon^2 u_2|_{H^1(\Omega^{\epsilon})} \leq |R_{u_{\epsilon}}|_{H^1(\Omega^{\epsilon})} + \|H_{\epsilon}\|_{L^2(\Omega^{\epsilon})} \leq C\epsilon^{\frac{3}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})},
$$

which is precisely the $H^1$-estimate (16). The $L^2$-estimate (17) is then a consequence of the Poincaré inequality in perforated domains (Lemma 4.6). To prove the remaining pressure estimate (15), recall that pressure is always defined up to a constant and choose the original pressure $p_\epsilon$, the homogenized pressure $p^*$ and the pressure estimator $p_{\epsilon,1}$ with zero mean in $\Omega^\epsilon$ (if it is not the case, it suffices to consider $p_{\epsilon,1} = p^* + \epsilon(p_1 - \int_{\Omega^\epsilon} p_1)$). Lemma 4.9 ensures the existence of $v \in H_0^1(\Omega^\epsilon)^d$ such that $\operatorname{div} v = p_\epsilon - p_{\epsilon,1}$. Then by integration by parts:

$$
\begin{align*}
\int_{\Omega^{\epsilon}} (p_{\epsilon} - p_{\epsilon,1})^2 &= \int_{\Omega^{\epsilon}} (p_{\epsilon} - p_{\epsilon,1}) \operatorname{div} v = - \int_{\Omega^{\epsilon}} F_{\epsilon} \cdot v + \int_{\Omega^{\epsilon}} \nu \nabla (u_{\epsilon} - \epsilon^2 u_2) \cdot \nabla v \\
&\leq C \epsilon^{\frac{3}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} |v|_{H^1(\Omega^{\epsilon})} \leq C \epsilon^{\frac{1}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \|p_{\epsilon} - p_{\epsilon,1}\|_{L^2(\Omega^{\epsilon})},
\end{align*}
$$

using the estimate in Lemma 4.9. Thus,

$$
\|p_{\epsilon} - p^{*}\|_{L^{2}(\Omega^{\epsilon})} \leq \|p_{\epsilon} - p_{\epsilon,1}\|_{L^{2}(\Omega^{\epsilon})} + \|p_{\epsilon,1} - p^{*}\|_{L^{2}(\Omega^{\epsilon})} \leq C\epsilon^{\frac{1}{2}}\|f - \nabla p^{*}\|_{H^{2}(\Omega) \cap C^{1,\alpha}(\bar{\Omega})},
$$

since $(p_{\epsilon,1} - p^*)$ is of order $\epsilon$ as seen from (30), which concludes the proof of Theorem 4.1.

References

[1] Grégoire Allaire. Homogenization of the Stokes flow in a connected porous medium. *Asymp-totic Analysis*, 2(3):203–222, 1989.

[2] Grégoire Allaire and Micol Amar. Boundary layer tails in periodic homogenization. *ESAIM, Control Optim. Calc. Var.*, 4:209–243, 1999.

[3] Alain Bensoussan, Jacques-Louis Lions, and George Papanicolaou. *Asymptotic analysis for periodic structures*. Providence, RI: AMS Chelsea Publishing, reprint of the 1978 original with corrections and bibliographical additions edition, 2011.

[4] Lamberto Cattabriga. Su un problema al contorno relativo al sistema di equazioni di Stokes. *Rend. Semin. Mat. Univ. Padova*, 31:308–340, 1961.

[5] G. A. Chechkin, A. L. Piatnitski, and A. S. Shamaev. *Homogenization. Methods and applications. Transl. from the Russian by Tamara Rozhkovskaya*, volume 234 of *Transl. Math. Monogr.* Providence, RI: American Mathematical Society (AMS), 2007.
Since the residual $R_{u_ε}$ is divergence-free and belongs to $H_0^1(\Omega^\epsilon)^d$, it can be used as test function, leading to

$$
\begin{align*}
\nu \int_{\Omega^\epsilon} \nabla R_{u_\epsilon} \cdot \nabla R_{u_\epsilon} &= \nu \int_{\Omega^\epsilon} \nabla R_{u_\epsilon} \cdot \nabla R_{u_\epsilon} - \int_{\Omega^\epsilon} R_{p_\epsilon} \operatorname{div} R_{u_\epsilon} \\
&= \int_{\Omega^\epsilon} F_\epsilon \cdot R_{u_\epsilon} - \int_{\Omega^\epsilon} H_\epsilon \cdot \nabla R_{u_\epsilon} \\
&\le \|F_\epsilon\|_{L^2(\Omega^\epsilon)} \|R_{u_\epsilon}\|_{L^2(\Omega^\epsilon)} + \|H_\epsilon\|_{L^2(\Omega^\epsilon)} \|R_{u_\epsilon}\|_{H^1(\Omega^\epsilon)}.
\end{align*}
$$

Using estimate (40) and the Poincaré inequality in a perforated domain (Lemma 4.6), it follows:

$$
\| \mathbf{F}_{\varepsilon} \|_{L^2(\Omega^{\varepsilon})} \| \mathbf{R}_{u_{\varepsilon}} \|_{L^2(\Omega^{\varepsilon})} \leq C \varepsilon^2 \| \mathbf{f} - \nabla p^* \|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})} |\mathbf{R}_{u_{\varepsilon}}|_{H^1(\Omega^{\varepsilon})}.
$$

Then, using (41), we conclude that:

$$
|R_{u_\epsilon}|_{H^1(\Omega^\epsilon)} \leq C\epsilon^{\frac{3}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1,\alpha}(\bar{\Omega})} \quad (45)
$$

Now, writing $u_\epsilon - \epsilon^2 u_2 = R_{u_\epsilon} + u_{\epsilon,3} - \epsilon^2 u_2$, using (45) and (41) yields:

$$
|u_{\epsilon} - \epsilon^2 u_2|_{H^1(\Omega^{\epsilon})} \leq |R_{u_{\epsilon}}|_{H^1(\Omega^{\epsilon})} + \|H_{\epsilon}\|_{L^2(\Omega^{\epsilon})} \leq C \epsilon^{\frac{3}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})},
$$

which is precisely the $H^1$-estimate (16). The $L^2$-estimate (17) is then a consequence of the Poincaré inequality in perforated domains (Lemma 4.6). To prove the remaining pressure estimate (15), recall that pressure is always defined up to a constant and choose the original pressure $p_\epsilon$, the homogenized pressure $p^*$ and the pressure estimator $p_{\epsilon,1}$ with zero mean in $\Omega^\epsilon$ (if it is not the case, it suffices to consider $p_{\epsilon,1} = p^* + \epsilon(p_1 - \int_{\Omega^\epsilon} p_1)$). Lemma 4.9 ensures the existence of $v \in H_0^1(\Omega^\epsilon)^d$ such that $\operatorname{div} v = p_\epsilon - p_{\epsilon,1}$. Then by integration by parts:

$$
\begin{align*}
\int_{\Omega^{\epsilon}} (p_{\epsilon} - p_{\epsilon,1})^2 &= \int_{\Omega^{\epsilon}} (p_{\epsilon} - p_{\epsilon,1}) \operatorname{div} v = - \int_{\Omega^{\epsilon}} F_{\epsilon} \cdot v + \int_{\Omega^{\epsilon}} \nu \nabla (u_{\epsilon} - \epsilon^2 u_2) \cdot \nabla v \\
&\leq C \epsilon^{\frac{3}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})} |v|_{H^1(\Omega^{\epsilon})} \leq C \epsilon^{\frac{1}{2}} \|f - \nabla p^*\|_{H^2(\Omega) \cap C^{1, \alpha}(\bar{\Omega})} \|p_{\epsilon} - p_{\epsilon,1}\|_{L^2(\Omega^{\epsilon})},
\end{align*}
$$

using the estimate in Lemma 4.9. Thus,

$$
\|p_{\epsilon} - p^{*}\|_{L^{2}(\Omega^{\epsilon})} \leq \|p_{\epsilon} - p_{\epsilon,1}\|_{L^{2}(\Omega^{\epsilon})} + \|p_{\epsilon,1} - p^{*}\|_{L^{2}(\Omega^{\epsilon})} \leq C\epsilon^{\frac{1}{2}}\|f - \nabla p^{*}\|_{H^{2}(\Omega) \cap C^{1,\alpha}(\bar{\Omega})},
$$

since $(p_{\epsilon,1} - p^*)$ is of order $\epsilon$ as seen from (30), which concludes the proof of Theorem 4.1.

References

[1] Grégoire Allaire. Homogenization of the Stokes flow in a connected porous medium. *Asymp-totic Analysis*, 2(3):203–222, 1989.

[2] Grégoire Allaire and Micol Amar. Boundary layer tails in periodic homogenization. *ESAIM, Control Optim. Calc. Var.*, 4:209–243, 1999.

[3] Alain Bensoussan, Jacques-Louis Lions, and George Papanicolaou. *Asymptotic analysis for periodic structures*. Providence, RI: AMS Chelsea Publishing, reprint of the 1978 original with corrections and bibliographical additions edition, 2011.

[4] Lamberto Cattabriga. Su un problema al contorno relativo al sistema di equazioni di Stokes. *Rend. Semin. Mat. Univ. Padova*, 31:308–340, 1961.

[5] G. A. Chechkin, A. L. Piatnitski, and A. S. Shamaev. *Homogenization. Methods and applications. Transl. from the Russian by Tamara Rozhkovskaya*, volume 234 of *Transl. Math. Monogr.* Providence, RI: American Mathematical Society (AMS), 2007.
[6] Bernard Dacorogna. Existence and regularity of solutions of $d\omega = f$ with Dirichlet boundary conditions. In *Nonlinear problems in mathematical physics and related topics I. In honor of Professor O. A. Ladyzhenskaya*, pages 67–82. New York, NY: Kluwer Academic/Plenum Publishers, 2002.

[7] Qingqing Feng, Gregoire Allaire, and Pascal Omnes. Enriched nonconforming multiscale finite element method for Stokes flows in heterogeneous media based on high-order weighting functions. *Multiscale Model. Simul.*, 20(1):462–492, 2022.

[8] David Gilbarg and Neil S. Trudinger. *Elliptic Partial Differential Equations of Second Order*, volume 224 of *Classics in Mathematics*. Springer, Berlin, Heidelberg, 2001.

[9] Vivette Girault and Pierre-Arnaud Raviart. *Finite element methods for Navier-Stokes equations. Theory and algorithms.*, volume 5 of *Springer Ser. Comput. Math*. Springer, Cham, (Extended version of the 1979 publ.) edition, 1986.

[10] Ulrich Hornung, editor. *Homogenization and porous media*, volume 6 of *Interdiscip. Appl. Math*. New York, NY: Springer, 1997.

[11] Gaspard Jankowiak and Alexei Lozinski. Non-conforming multiscale finite element method for Stokes flows in heterogeneous media. part II: Error estimates for periodic microstructure. *Discrete and Continuous Dynamical Systems - B*, 29(5):2298–2332, 2024.

[12] Eduard Marusic-Paloka and Andro Mikelić. An Error Estimate for Correctors in the Homogenization of the Stokes and the Navier-Stokes Equations in a Porous Medium. *Boll. Un. Mat. Ital. A (7)*, 10:661–671, January 1996.

[13] Évariste Sanchez-Palencia. *Non-homogeneous media and vibration theory*. Number 127 in Lecture notes in physics. Springer, Berlin, 1980.

[14] Zhongwei Shen. Sharp convergence rates for Darcy’s law. *Communications in Partial Differential Equations*, 47(6):1098–1123, June 2022.

[15] Luc Tartar. Convergence of the homogenization process. *Appendix of [13]*, 1980.
[6] Bernard Dacorogna. Existence and regularity of solutions of $d\omega = f$ with Dirichlet boundary conditions. In *Nonlinear problems in mathematical physics and related topics I. In honor of Professor O. A. Ladyzhenskaya*, pages 67–82. New York, NY: Kluwer Academic/Plenum Publishers, 2002.

[7] Qingqing Feng, Gregoire Allaire, and Pascal Omnes. Enriched nonconforming multiscale finite element method for Stokes flows in heterogeneous media based on high-order weighting functions. *Multiscale Model. Simul.*, 20(1):462–492, 2022.

[8] David Gilbarg and Neil S. Trudinger. *Elliptic Partial Differential Equations of Second Order*, volume 224 of *Classics in Mathematics*. Springer, Berlin, Heidelberg, 2001.

[9] Vivette Girault and Pierre-Arnaud Raviart. *Finite element methods for Navier-Stokes equations. Theory and algorithms.*, volume 5 of *Springer Ser. Comput. Math*. Springer, Cham, (Extended version of the 1979 publ.) edition, 1986.

[10] Ulrich Hornung, editor. *Homogenization and porous media*, volume 6 of *Interdiscip. Appl. Math*. New York, NY: Springer, 1997.

[11] Gaspard Jankowiak and Alexei Lozinski. Non-conforming multiscale finite element method for Stokes flows in heterogeneous media. part II: Error estimates for periodic microstructure. *Discrete and Continuous Dynamical Systems - B*, 29(5):2298–2332, 2024.

[12] Eduard Marusic-Paloka and Andro Mikelić. An Error Estimate for Correctors in the Homogenization of the Stokes and the Navier-Stokes Equations in a Porous Medium. *Boll. Un. Mat. Ital. A (7)*, 10:661–671, January 1996.

[13] Évariste Sanchez-Palencia. *Non-homogeneous media and vibration theory*. Number 127 in Lecture notes in physics. Springer, Berlin, 1980.

[14] Zhongwei Shen. Sharp convergence rates for Darcy’s law. *Communications in Partial Differential Equations*, 47(6):1098–1123, June 2022.

[15] Luc Tartar. Convergence of the homogenization process. *Appendix of [13]*, 1980.