# On the analysis of a mechanically consistent model of fluid-structure-contact interaction

Marguerite Champion, Miguel Angel Fernández, Céline Grandmont, Fabien

Vergnet, Marina Vidrascu

► To cite this version:

Marguerite Champion, Miguel Angel Fernández, Céline Grandmont, Fabien Vergnet, Marina Vidrascu. On the analysis of a mechanically consistent model of fluid-structure-contact interaction. Mathematics in Engineering, 2024, 6 (3), pp.425-467. 10.3934/mine.2024018. hal-04229012

HAL Id: hal-04229012

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

Submitted on 5 Oct 2023

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On the analysis of a mechanically consistent model of fluid-structure-contact interaction

## Marguerite Champion

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France
* marguerite.champion@inria.fr

## Miguel A. Fernández

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France
* miguel.fernandez@inria.fr

## Céline Grandmont

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France
* Département de Mathématique, Université Libre de Bruxelles, Bruxelles, Belgium
* celine.grandmont@inria.fr

## Fabien Vergnet

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France
* fabien.vergnet@sorbonne-universite.fr

## Marina Vidrascu

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France
* marina.vidrascu@inria.fr

Received (Day Month Year)

Revised (Day Month Year)

Communicated by (xxxxxxxxxx)

This paper is devoted to the mathematical analysis of the contact capabilities of the fluid-structure interaction (FSI) model with seepage reported in [*Comput. Methods Appl. Mech., 392:114637, 2022*]. In the case of a rigid disk moving over a fixed horizontal plane, we show that this model encompasses contact and hence removes the non collision paradox of traditional FSI models which rely on Dirichlet or Dirichlet/Navier boundary conditions. Numerical evidence on the theoretical results is also provided.

**Keywords:** Fluid-structure interaction, collision paradox, Darcy flow, Navier boundary conditions.

AMS Subject Classification: 35Q35, 35B44, 74F10, 76D03.

# 1. Introduction

The numerical simulation of systems involving fluid-structure-contact interaction is of fundamental importance in many engineering and biomedical applications. For
2 *Champion, Fernández, Grandmont, Vergnet, Vidrascu*

instance, contact modeling is a crucial ingredient in the simulation of the dynamics
of native or artificial cardiac valves (see, e.g., 19, 20).

Modeling contact between solids within a fluid-structure interaction (FSI) frame-
work raises many modeling, mathematical and numerical issues. First, in the case
of a ball immersed in a viscous incompressible fluid with no-slip boundary condi-
tions and falling over a fixed horizontal plane, the resulting FSI models are unable
to predict contact, both in 2D (see, e.g., 14) and in 3D (see, e.g.,16), which is
known as the *no collision paradox*. One of the most widespread explanation of
this paradox is that one can no longer consider ideal smooth surfaces when solids
come into contact: roughness-induced effects play a fundamental role into enabling
collision. Indeed, contact between smooth solids with no-slip boundary conditions
seem possible only in very specific configurations like, for example, grazing collision
in 3D (see 17). On the other hand, many studies show the *no collision paradox* is
circumvented by taking into account roughness in FSI models, either by enabling
the fluid to slip through Navier-type boundary conditions (see, e.g., 22, 18 in 2D
and 9, 12, 10 in 3D) or by considering rough solids (see, e.g., 8 in 2D and 9
in 3D). A second major difficulty is related to the mechanical consistency of the
model. Indeed, for FSI models which allow for contact, the simple addition of a
non-penetration constraint to the solid can lead to unphysical void creation (at re-
lease from contact) or unbalanced stresses at contact. Recently, these mechanical
inconsistencies have been avoided by considering a poroelastic modeling of the fluid
seepage induced by the roughness of the contacting wall (see, e.g., 1, 4). Yet, very
little is known on the mathematical foundations of these modeling approaches.

In this work, we investigate the capability to encompass contact of the FSI
model with seepage reported in 4. For this purpose, we consider a simplified 2D
setting of a rigid disk immersed in a Stokesian flow and falling over an horizontal
plane (the contacting wall), modeled as porous layer. We provide a well-posedness
analysis for the fluid problem and describe the asymptotics with respect to the
porous layer parameters. By building on the arguments reported in 9, we also derive
an estimate of the fluid drag force, acting on the disk boundary, with respect to the
gap between the disk and the porous layer. A salient feature of this analysis is that
it shows that the considered FSI model with seepage allows for contact between
the disk and the wall. In other words, since the porous layer allows for seepage, the
incompressibility constraint does not create any singularity which prevents contact.
To the best of our knowledge, this is the first time in which contact is allowed in a
FSI model with Dirichlet interface conditions on the falling disk, for the considered
geometrical setting. From the analysis reported in 14, 9, one can indeed show that
the combination of Dirichlet and Navier boundary conditions prevents contact. The
mathematical analysis of the paper is complemented by a comprehensive numerical
study which illustrates the theoretical results obtained.

The rest of the paper is organized as follows. Section 2 presents the considered
2D simplified setting and the fluid-structure-contact interaction model of 4. An ap-
propriate velocity scaling of the problem (see, e.g., 9) reduces the coupled problem
to the evaluation of the drag force operator in terms of the distance to the contacting wall. Section 3 represents the theoretical core of the paper. It provides a mathematical and asymptotic analysis of the resulting scaled fluid system and analytical estimations of the drag force, which allow to describe the contact dynamics of the disk. Numerical evidence on these theoretical results is provided in Section 4. Finally, the main results of the paper are summarized in Section 5 together with some perspectives of future work.

## 2. Mathematical models

In order to investigate the effect of the surface porous medium on the contact dynamics of fluid-structure interaction, we consider a simplified fluid-structure interaction system involving a quasi-steady Stokes flow with an immersed rigid disk. This simplified setting has already been investigated in previous studies (see, e.g., 14, 9, 18), but with a different treatment of the contact wall (notably in terms of boundary conditions).

The Stokesian fluid is assumed to be contained in a rectangular domain $\Omega \stackrel{\text{def}}{=} (-L, L) \times (0, \tilde{L})$ and the current configuration of the immersed rigid disk is denoted by $S(t) \subset \Omega$, for all $t > 0$. We can hence introduce the fluid domain $\Omega(t) \stackrel{\text{def}}{=} \Omega \setminus S(t)$ and its associated non-cylindrical trajectory

$$ \mathcal{T} \stackrel{\text{def}}{=} \bigcup_{t \in \mathbb{R}^+} \Omega(t) \times \{t\}. $$

The boundary of $\Omega$ is partitioned as $\partial\Omega = \Gamma \cup \Sigma$, where $\Sigma$ denotes the bottom contacting wall (see Figure 1). In what follows, the rigid disk is assumed to be of radius one and to move only vertically (without rotation), so that $S(t)$ can be defined as follows:

$$ S(t) \stackrel{\text{def}}{=} B((d(t) + 1)e_2, 1), $$

where $(e_1, e_2)$ denotes the canonical basis of $\mathbb{R}^2$, $B(x, R)$ the disk of radius $R$ centered in $x \in \mathbb{R}^2$, and $d: \mathbb{R}^+ \to \mathbb{R}$ stands for the so so called gap function, viz., the quantity $d(t)$ represents the (signed) distance function between $S(t)$ and the contacting wall $\Sigma$. Owing to the geometrical symmetry of the problem, it is worth noting that this purely vertical motion can be (physically) expected under appropriate initial conditions.

**The fluid.** The state of the fluid can be described in terms of its velocity $u: \mathcal{T} \to \mathbb{R}^2$ and pressure $p: \mathcal{T} \to \mathbb{R}$ fields, which are governed by the following Stokes system

$$ \left\{ \begin{array}{ll} -\text{div } \sigma(u, p) = 0 & \text{in } \Omega(t), \\ \text{div } u = 0 & \text{in } \Omega(t), \\ u = 0 & \text{on } \Gamma, \end{array} \right. \qquad (2.1) $$

where the fluid Cauchy stress tensor $\sigma$ reads

$$ \sigma(u, p) \stackrel{\text{def}}{=} 2D(u) - pI, \quad D(u) \stackrel{\text{def}}{=} \frac{1}{2}(\nabla u + \nabla u^T). $$
4 Champion, Fernández, Grandmont, Vergnet, Vidrascu

<figure><img src="image_2.png" /><figcaption>Figure 1: Geometrical description.</figcaption></figure>

The fluid viscosity is assumed to be one, and no-slip boundary conditions are considered on the exterior boundary $\Gamma$.

**The solid.** Owing to the vertical motion of the rigid disk, the dynamics of the solid is simply given in terms of the momentum conservation law. Assuming that the solid has unit mass, the following relation holds

$$ \ddot{d}(t) + \mathcal{F}(t) = 0, \qquad (2.2) $$

for $t \in \mathbb{R}^+$ and where $\mathcal{F}(t)$ stands for the external force acting on the solid.

**Fluid-solid coupling.** The fluid and solid equations (2.1)-(2.2) have to be coupled with standard geometric, kinematic and dynamic interface conditions on $\partial S(t)$, namely:

$$ \left\{ \begin{array}{l} S(t) = B((d(t) + 1)e_2, 1), \quad \Omega(t) = \Omega \setminus S(t), \\ \mathbf{u} = \dot{d}e_2 \quad \text{on} \quad \partial S(t), \\ \mathcal{F}(t) = \int_{\partial S(t)} \sigma(\mathbf{u}, p)\mathbf{n} \cdot \mathbf{e}_2, \end{array} \right. \qquad (2.3) $$

where $\mathbf{n}$ denotes the exterior unit normal to $\partial \Omega(t)$. Note that the last expression simply relates $\mathcal{F}(t)$ to the fluid drag force exerted on the rigid disk at time $t$.

**Surface roughness model on the contact wall.** Traditionally, fluid-structure interaction models based on (2.1)-(2.3) involve a no-slip boundary condition on the contact wall $\Sigma$, which is known to prevent contact (see, e.g. 14) in contrast to what is physically observed. In order to circumvent this issue, we consider the alternative modeling approach reported in 4, which consists in taking into account the roughness of the contact wall through a surface model of fluid seepage.

The matrix of the porous wall is assumed to be rigid and of thickness $\epsilon_p$. The
normal and tangential conductivities of the porous medium are denoted by $\kappa$ and $\eta$, respectively. We denote by $\boldsymbol{\tau}$ the tangential vector on the surface $\Sigma$; here $\boldsymbol{\tau}$ is simply equal to $e_1$. The porous fluid pressure $\hat{p} : \Sigma \times \mathbb{R}^+ \to \mathbb{R}$ can be described by the following surface Darcy model (obtained after averaging the bulk Darcy problem across the thickness, see 4, 21):

$$ \left\{ \begin{array}{ll} -\operatorname{div}_{\boldsymbol{\tau}} (\varepsilon_p \eta \nabla_{\boldsymbol{\tau}} \hat{p}) = u_{\Sigma} & \text{on } \Sigma, \\ \nabla_{\boldsymbol{\tau}} \hat{p} = 0 & \text{on } \partial\Sigma, \end{array} \right. \qquad (2.4) $$

where $u_{\Sigma}$ denotes the seepage velocity.

**Remark 2.1.** As $\Sigma$ is a 1D line, the symbols $\operatorname{div}_{\boldsymbol{\tau}}$ and $\nabla_{\boldsymbol{\tau}}$, denoting respectively the divergence and gradient along the direction $\boldsymbol{\tau}$, are here both equivalent to $\frac{\partial}{\partial x}$.

The above system has to be coupled with the bulk fluid equations (2.1) on $\Sigma$, via appropriate kinematic and dynamic relations, namely,

$$ \left\{ \begin{array}{ll} u_{\Sigma} = \boldsymbol{u} \cdot \boldsymbol{n} & \text{on } \Sigma, \\ \sigma(\boldsymbol{u}, p)\boldsymbol{n} \cdot \boldsymbol{n} = -\left(\hat{p} + \frac{\varepsilon_p}{4\kappa} \boldsymbol{u} \cdot \boldsymbol{n}\right) & \text{on } \Sigma, \\ \sigma(\boldsymbol{u}, p)\boldsymbol{n} \cdot \boldsymbol{\tau} = 0 & \text{on } \Sigma. \end{array} \right. \qquad (2.5) $$

Note that free slip tangential stress is enforced with the last relation.

**Fully coupled problem.** The considered coupled problem can be summarized as follows: Find the gap $d: \mathbb{R}^+ \to \mathbb{R}$, the fluid velocity $\boldsymbol{u}: \mathcal{T} \to \mathbb{R}^2$, the fluid pressure $p: \mathcal{T} \to \mathbb{R}$ and the porous pressure $\hat{p}: \Sigma \times \mathbb{R} \to \mathbb{R}$ such that:

$$ \left\{ \begin{array}{ll} -\operatorname{div} \boldsymbol{\sigma}(\boldsymbol{u}, p) = \mathbf{0} & \text{in } \Omega(t), \qquad (2.6a) \\ \operatorname{div} \boldsymbol{u} = 0 & \text{in } \Omega(t), \qquad (2.6b) \\ \boldsymbol{u} = \mathbf{0} & \text{on } \Gamma, \qquad (2.6c) \\ \boldsymbol{u} = \dot{e}_2 & \text{on } \partial S(t), \qquad (2.6d) \\ \boldsymbol{\sigma}(\boldsymbol{u}, p)\boldsymbol{n} \cdot \boldsymbol{n} = -\left(\hat{p} + \frac{\varepsilon_p}{4\kappa} \boldsymbol{u} \cdot \boldsymbol{n}\right) & \text{on } \Sigma, \qquad (2.6e) \\ \boldsymbol{\sigma}(\boldsymbol{u}, p)\boldsymbol{n} \cdot \boldsymbol{\tau} = 0 & \text{on } \Sigma, \qquad (2.6f) \\ -\operatorname{div}_{\boldsymbol{\tau}} (\eta \varepsilon_p \nabla_{\boldsymbol{\tau}} \hat{p}) = \boldsymbol{u} \cdot \boldsymbol{n} & \text{on } \Sigma, \qquad (2.6g) \\ \nabla_{\boldsymbol{\tau}} \hat{p} \cdot \boldsymbol{\tau} = 0 & \text{on } \partial \Sigma, \qquad (2.6h) \\ \\ \ddot{d}(t) + \mathcal{F}(t) = 0, \quad \mathcal{F}(t) = \int_{\partial S(t)} \boldsymbol{\sigma}(\boldsymbol{u}, p)\boldsymbol{n} \cdot \boldsymbol{e}_2, \qquad (2.6i) \\ \\ S(t) = B((d(t) + 1)e_2, 1), \quad \Omega(t) = \Omega \setminus S(t), \qquad (2.6j) \end{array} \right. $$

with $d(0) = d_0$ and $\dot{d}(0) = \dot{d}_0$ as initial conditions.
6 Champion, Fernández, Grandmont, Vergnet, Vidrascu

**Velocity scaling.** System (2.6) is a highly non-linear coupled problem, notably due to the fact that the fluid domain $\Omega(t)$ is unknown. In order to mitigate this difficulty, we adopt here a scaling argument, considered for instance in 9, which leverages the linearity of problem (2.6) with respect to $\dot{d}$ and basically consists in performing the following change of variables for $t > 0$:

$$ egin{cases} 
\mathbf{u}(t) = \dot{d}(t)\mathbf{u}_{d(t)} & \text{in } \Omega_{d(t)}, \\ 
p(t) = \dot{d}(t)p_{d(t)} & \text{in } \Omega_{d(t)}, \\ 
\hat{p}(t) = \dot{d}(t)\hat{p}_{d(t)} & \text{in } \Omega_{d(t)}. 
\end{cases} 
\quad (2.7) $$

By inserting the above relations into (2.6) and since the equations are linear for a given location of the rigid disk, we infer that, for a given gap $d \in (0, \tilde{L} - 2)$ and associated solid and fluid domains

$$ S_d \stackrel{\text{def}}{=} B((d+1)e_2, 1), \quad \Omega_d \stackrel{\text{def}}{=} \Omega \setminus S_d, $$

the triplet $\mathbf{u}_d : \Omega_d \to \mathbb{R}^2, p_d : \Omega_d \to \mathbb{R}, \hat{p}_d : \Sigma \to \mathbb{R}$ is solution of the following pure (steady) fluid coupled problem:

$$ egin{cases} 
- \mathrm{div} \, \boldsymbol{\sigma}(\mathbf{u}_d, p_d) = \mathbf{0} & \text{in } \Omega_d, \quad (2.8a) \\ 
\mathrm{div} \, \mathbf{u}_d = 0 & \text{in } \Omega_d, \quad (2.8b) \\ 
\mathbf{u}_d = 0 & \text{on } \Gamma, \quad (2.8c) \\ 
\mathbf{u}_d = e_2 & \text{on } \partial S_d, \quad (2.8d) \\ 
\boldsymbol{\sigma}(\mathbf{u}_d, p_d)\mathbf{n} \cdot \mathbf{n} = -\left(\hat{p}_d + \frac{\varepsilon_p}{4\kappa}\mathbf{u}_d \cdot \mathbf{n}\right) & \text{on } \Sigma, \quad (2.8e) \\ 
\boldsymbol{\sigma}(\mathbf{u}_d, p_d)\mathbf{n} \cdot \boldsymbol{\tau} = 0 & \text{on } \Sigma, \quad (2.8f) \\ 
- \mathrm{div}_{\boldsymbol{\tau}}(\eta\varepsilon_p \nabla_{\boldsymbol{\tau}}\hat{p}_d) = \mathbf{u}_d \cdot \mathbf{n} & \text{on } \Sigma, \quad (2.8g) \\ 
\nabla_{\boldsymbol{\tau}}\hat{p}_d \cdot \boldsymbol{\tau} = 0 & \text{on } \partial\Sigma. \quad (2.8h) 
\end{cases} $$

Note that the velocity scaling (2.7) yields a constant unit velocity as kinematic constraint in (2.8d), and hence independent of $\dot{d}$. The analysis reported in Section 3.1 below shows that the operator $d \in (0, \tilde{L} - 2) \mapsto (\mathbf{u}_d, p_d, \hat{p}_d)$ is well defined, so that the change of variables (2.7) is feasible.

By inserting the scaling (2.7) into (2.6i), it follows by linearity that

$$ \mathcal{F}(t) = \dot{d}(t)\mathcal{F}_{d(t)}, $$

with the notation

$$ \mathcal{F}_d \stackrel{\text{def}}{=} \int_{\partial S_d} \boldsymbol{\sigma}(\mathbf{u}_d, p_d)\mathbf{n} \cdot \mathbf{e}_2 \quad (2.9) $$

for any given $d \in (0, \tilde{L} - 2)$. As a result, the rigid disk dynamics can be rewritten as:

$$ \ddot{d}(t) + \dot{d}(t)\mathcal{F}_{d(t)} = 0. \quad (2.10) $$
We will show in Section 3.3 that $\mathcal{F}_d \ge 0$, so that from (2.10) we can infer that the fluid-solid coupling (2.6) acts as a friction force into the dynamics of the rigid disk. The amount of this friction is nothing but the scaled drag force (2.9), which only depends on the instantaneous gap $d$, trough the solution of the pure fluid problem (2.8). The analysis of the system (2.8) and, more specifically, of the behavior of its associated drag force $\mathcal{F}_d$ with respect to $d$ provides a practical way to estimate the disk dynamics near the wall and, in particular, to conclude on the possibility of contact between the rigid disk and the bottom boundary $\Sigma$.

Previous studies (see, e.g., 14, 9, 10, 22, 16) have shown that Dirichlet or Navier boundary conditions on the contact wall $\Sigma$ lead to a singularity of $\mathcal{F}_d$ when $d \to 0$. In 2D with Dirichlet boundary conditions on both the wall and the disk or a Navier boundary conditions on one side and a Dirichlet boundary conditions on the other side, the singularity of $\mathcal{F}_d$ is

$$ \mathcal{F}_d \sim_{d \to 0} d^{-\frac{3}{2}} $$

and it prevents collision, as proven in Theorem 1 of the seminal paper 14. If Navier boundary conditions are considered on both the wall and the disk, the singularity is weaker (see Proposition 6.1 of 22) and allows collision:

$$ \mathcal{F}_d \sim_{d \to 0} d^{-\frac{1}{2}}. $$

In this paper, we investigate the case in which Dirichlet boundary conditions are applied on the disk and a porous layer is considered on the bottom wall. To the best of our knowledge, this new setting has not yet been investigated from a mathematical analysis perspective. In the next section, we prove existence and uniqueness of weak solution to the coupled system (2.8). Then we study the asymptotic behavior of the system with respect to the porous layer parameters $\kappa$ and $\eta$. Finally, we provide an estimate of the drag force when $d \to 0$ and show that

$$ 0 \le \mathcal{F}_d \le C(\kappa, \eta), $$

from which we conclude that the porous layer enables contact in the 2D case.

### 3. Mathematical analysis

This section gathers the main theoretical results of the paper. We provide a thorough mathematical analysis of system (2.8), whose main distinctive feature lies in the coupling with the porous layer on the bottom wall $\Sigma$. As, to our knowledge, this system has not yet been studied, we address its well-posedness in Section 3.1, using standard arguments from saddle point problems theory. Then, in order to highlight the role of the normal and tangential conductivity parameters $\kappa$ and $\eta$, we provide an asymptotic analysis of (2.8) with respect to these parameters in Section 3.2. Finally, in Section 3.3 we investigate the influence of the porous layer on the disk contact dynamics, by estimating analytically the drag force acting on the disk.
### 3.1. Existence and uniqueness

In this section, the existence and uniqueness of weak solution for system (2.8) is investigated. We also derive a priori estimates satisfied by this weak solution.

Let $\omega$ be a given bounded domain of $\mathbb{R}^2$ and let $\Upsilon$ denote a subset of its boundary. We denote by $L_0^2(\omega)$ the space of $L^2(\omega)$ functions with zero mean value in $\omega$ and by $H_{0,\Upsilon}^1(\omega)$ the space of functions in $H^1(\omega)$ with zero trace on $\Upsilon$. The scalar product in $L^2(\omega)$ is denoted by $(\cdot, \cdot)_{0,\omega}$, the $L^2(\omega)$ norm by $\|\cdot\|_{0,\omega}$ and the $H^1(\omega)$ norm by $\|\cdot\|_{1,\omega}$. We recall that the $H^1$ semi-norm $|\cdot|_{1,\omega}$ is a norm in $H_{0,\Upsilon}^1(\omega)$, since the Poincaré inequality holds true for functions in $H_{0,\Upsilon}^1(\omega)$. On $\Upsilon$, we also consider the $H^s$ norms, denoted by $\|\cdot\|_{s,\Upsilon}$, for every $0 < s \le 1$, and the $H^1$ semi-norm, denoted by $|\cdot|_{1,\Upsilon}$.

Since only the gradient of $p_d$ and of $\hat{p}_d$ is involved in (2.8a) and (2.8g), respectively, both $p_d$ and $\hat{p}_d$ are in principle defined up to a constant. Nonetheless, owing to (2.8e), once one of the constants is fixed, the other is uniquely defined, therefore, only one of them needs to be fixed to guarantee uniqueness. In the analysis below, we choose to fix the constant of the Darcy pressure $\hat{p}_d$ by considering the space

$$ \mathcal{D} \stackrel{\text{def}}{=} H^1(\Sigma) \cap L_0^2(\Sigma), \quad (3.1) $$

so that the fluid pressure belongs to $L^2(\Omega_d)$. As regards the fluid velocity, we denote by $\mathcal{U}_d$ and $\mathcal{U}_d^0$ the following spaces,

$$ \mathcal{U}_d \stackrel{\text{def}}{=} \{ \mathbf{u} \in H_\Gamma^1(\Omega_d)^2 : \mathbf{u} = e_2 \text{ on } \partial S_d \}, \quad \mathcal{U}_d^0 \stackrel{\text{def}}{=} H_{\Gamma \cup \partial S_d}^1(\Omega_d)^2. $$

In order to derive a weak formulation for problem (2.8), we test the fluid momentum equation (2.8a) with $\mathbf{v} \in \mathcal{U}_d^0$, the incompressibility condition (2.8b) with $q \in L^2(\Omega_d)$ and the surface Darcy equation (2.8g) with $\hat{q} \in \mathcal{D}$. The weak formulation then follows by integrating by parts and using the different boundary and coupling conditions: Find $(\mathbf{u}_d, \hat{p}_d, p_d) \in \mathcal{U}_d \times \mathcal{D} \times L^2(\Omega_d)$ such that

$$ \begin{aligned} & 2(\mathcal{D}(\mathbf{u}_d), \mathcal{D}(\mathbf{v}))_{0,\Omega_d} + \frac{\varepsilon_p}{4\kappa} (\mathbf{u}_d \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} + \varepsilon_p \eta (\nabla_{\tau} \hat{p}_d, \nabla_{\tau} \hat{q})_{0,\Sigma} \\ & + (\hat{p}_d, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} - (\mathbf{u}_d \cdot \mathbf{n}, \hat{q})_{0,\Sigma} - (p_d, \text{div} \, \mathbf{v})_{0,\Omega_d} + (\text{div} \, \mathbf{u}_d, q)_{0,\Omega_d} = 0 \quad (3.2) \end{aligned} $$

for all $(\mathbf{v}, \hat{q}, q) \in \mathcal{U}_d^0 \times \mathcal{D} \times L^2(\Omega_d)$.

The main results of this section are stated in the following theorem.

**Theorem 3.1.** For any $d \in (0, \tilde{L} - 2)$, problem (3.2) admits a unique solution $(\mathbf{u}_d, p_d, \hat{p}_d) \in \mathcal{U}_d \times L^2(\Omega_d) \times \mathcal{D}$. Moreover, there exist three positive constants $C_1, C_2, C_3 > 0$, which only depend on $\Omega_d$, such that $(\mathbf{u}_d, p_d, \hat{p}_d)$ satisfies the fol-
10 Champion, Fernández, Grandmont, Vergnet, Vidrascu

the standard arguments to prove existence and uniqueness. We consider the Hilbert spaces

$$X \stackrel{\text{def}}{=} \mathcal{U}_d^0 \times \mathcal{D}, \quad M \stackrel{\text{def}}{=} L^2(\Omega_d).$$

Let $a: X \times X \to \mathbb{R}$, $b: X \times M \to \mathbb{R}$ and $l \in X'$ be the bi-linear and linear forms defined by

$$ 
\begin{align} 
a((\mathbf{u}, \hat{\mathbf{p}}), (\mathbf{v}, \hat{\mathbf{q}})) &\stackrel{\text{def}}{=} 2 (\mathcal{D}(\mathbf{u}), \mathcal{D}(\mathbf{v}))_{0,\Omega_d} + \frac{\varepsilon_p}{4\kappa} (\mathbf{u} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} \nonumber \\ 
&\quad + \varepsilon_p \eta (\nabla_{\tau} \hat{\mathbf{p}}, \nabla_{\tau} \hat{\mathbf{q}})_{0,\Sigma} + (\hat{\mathbf{p}}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} - (\mathbf{u} \cdot \mathbf{n}, \hat{\mathbf{q}})_{0,\Sigma}, \tag{3.10} \\ 
b((\mathbf{v}, \hat{\mathbf{q}}), q) &\stackrel{\text{def}}{=} -(q, \text{div } \mathbf{v})_{0,\Omega_d}, \nonumber \\ 
\langle l, (\mathbf{v}, \hat{\mathbf{q}}) \rangle_{X',X} &\stackrel{\text{def}}{=} (\mathcal{D}(\mathbf{w}), \mathcal{D}(\mathbf{v}))_{0,\Omega_d}. \nonumber 
\end{align} 
$$

Problem (3.9) can therefore be reformulated in an abstract fashion as:

$$ \text{Find } ((\bar{\mathbf{u}}, \hat{\mathbf{p}}), p) \in X \times M \text{ such that} \quad \left\{ \begin{array}{ll} a((\bar{\mathbf{u}}, \hat{\mathbf{p}}), (\mathbf{v}, \hat{\mathbf{q}})) + b((\mathbf{v}, \hat{\mathbf{q}}), p_d) = \langle l, (\mathbf{v}, \hat{\mathbf{q}}) \rangle_{X',X} & \forall (\mathbf{v}, \hat{\mathbf{q}}) \in X, \\ b((\bar{\mathbf{u}}, \hat{\mathbf{p}}), q) = 0 & \forall q \in M. \end{array} \right. \quad (3.11) $$

The form $l$ is continuous on $X$, the form $b$ is continuous on $X \times M$, by continuity of the divergence operator and, using the continuity of the trace operator, we show that $a$ is continuous on $X \times X$. Moreover, using Körn's inequality (see, e.g., 7), it follows that form $a$ is coercive on $X \times X$.

**Inf-sup condition.** To prove that problem (3.11) is well-posed, it remains to show that the continuous bilinear form $b$ satisfies the usual inf-sup condition. We first observe that

$$ \sup_{(\mathbf{v}, \hat{\mathbf{q}}) \in \mathcal{U}_d^0 \times \mathcal{D} \setminus \{(0,0)\}} \frac{|(\text{div } \mathbf{v}, q)|}{\|(\mathbf{v}, \hat{\mathbf{q}})\|_{\mathcal{U}_d^0 \times \mathcal{D}} \|q\|_{0,\Omega_d}} \ge \sup_{\mathbf{v} \in \mathcal{U}_d^0 \setminus \{0\}} \frac{|(\text{div } \mathbf{v}, q)|}{\|\mathbf{v}\|_{1,\Omega_d} \|q\|_{0,\Omega_d}}, $$

for all $q \in L^2(\Omega_d) \setminus \{0\}$. From Theorem 3.3 of 6 and by adapting the Bogovskii's lemma stated in 2 to the case of homogeneous boundary condition on part of the boundary, there exists a constant $C > 0$ such that for any $q \in L^2(\Omega_d)$ there exists $\mathbf{w}_q \in \mathcal{U}_d^0$ such that $\text{div } \mathbf{w}_q = q$ and $\|\mathbf{w}_q\|_{1,\Omega_d} \le C \|q\|_{0,\Omega_d}$. As a result

$$ \sup_{\mathbf{v} \in \mathcal{U}_d^0 \setminus \{0\}} \frac{|(\text{div } \mathbf{v}, q)|}{\|\mathbf{v}\|_{1,\Omega_d}} \ge \frac{|(\text{div } \mathbf{w}_q, q)|}{\|\mathbf{w}_q\|_{1,\Omega_d}} = \frac{\|q\|_{0,\Omega_d}^2}{\|\mathbf{w}_q\|_{1,\Omega_d}} \ge \frac{1}{C} \|q\|_{0,\Omega_d}, $$

for all $q \in L^2(\Omega_d) \setminus \{0\}$, which yields the inf-sup condition

$$ \inf_{q \in L^2(\Omega_d) \setminus \{0\}} \sup_{(\mathbf{v}, \hat{\mathbf{q}}) \in \mathcal{U}_d^0 \times \mathcal{D} \setminus \{(0,0)\}} \frac{|(\text{div } \mathbf{v}, q)|}{\|(\mathbf{v}, \hat{\mathbf{q}})\|_{\mathcal{U}_d^0 \times \mathcal{D}} \|q\|_{0,\Omega_d}} \ge \frac{1}{C}. \quad (3.12) $$

Therefore, since the bilinear form $a$ is continuous and coercive on $X \times X$, that $l$ is continuous on $X$ and that $b$ is continuous on $X \times M$ and satisfies the inf-sup condition (3.12), we conclude that problem (3.11) admits a unique solution in $X \times M$ (see, e.g., Corollary 4.1 of 11). Moreover, it also proves that problem (3.2) admits a unique solution in $\mathcal{U}_d \times \mathcal{D} \times L^2(\Omega_d)$.
**A priori estimates.** To obtain a priori estimates on $\mathbf{u}$ and $\hat{\mathbf{p}}$, we first consider $(\bar{\mathbf{u}}, \hat{\mathbf{p}}) \in X$ as a test function in the weak formulation (3.9), which yields

$$2\|\mathbf{D}(\bar{\mathbf{u}})\|_{0,\Omega_d}^2 + \frac{\varepsilon_p}{4\kappa}\|\bar{\mathbf{u}} \cdot \mathbf{n}\|_{0,\Sigma}^2 + \varepsilon_p \eta \|\nabla_{\tau} \hat{\mathbf{p}}\|_{0,\Sigma}^2 = -2(\mathbf{D}(\mathbf{w}), \mathbf{D}(\bar{\mathbf{u}}))_{\Omega_d}.$$

Thus, from (3.6) and (3.8), we obtain

$$2\|\mathbf{D}(\mathbf{u})\|_{0,\Omega_d}^2 + \frac{\varepsilon_p}{4\kappa}\|\mathbf{u} \cdot \mathbf{n}\|_{0,\Sigma}^2 + \varepsilon_p \eta \|\nabla_{\tau} \hat{\mathbf{p}}\|_{0,\Sigma}^2 = 2(\mathbf{D}(\mathbf{w}), \mathbf{D}(\mathbf{u}))_{0,\Omega_d}.$$

Finally, by using the Cauchy-Schwarz and Young inequalities, we get

$$2\|\mathbf{D}(\mathbf{u})\|_{0,\Omega_d}^2 + \frac{\varepsilon_p}{4\kappa}\|\mathbf{u} \cdot \mathbf{n}\|_{0,\Sigma}^2 + \varepsilon_p \eta \|\nabla_{\tau} \hat{\mathbf{p}}\|_{0,\Sigma}^2 \le \|\mathbf{D}(\mathbf{w})\|_{0,\Omega_d}^2 + \|\mathbf{D}(\mathbf{u})\|_{0,\Omega_d}^2. \quad (3.13)$$

The a priori estimate (3.3) then simply follows from (3.13) and (3.7).

Thanks to (3.3) on $\mathbf{u}_d \cdot \mathbf{n}$, we can derive a sharper estimate on $\hat{\mathbf{p}}_d$. Indeed using the variational formulation (3.2) with $\mathbf{v} = \mathbf{0}$ and $\hat{\mathbf{q}} = \hat{\mathbf{p}}_d$, we obtain

$$\varepsilon_p \eta \|\nabla_{\tau} \hat{\mathbf{p}}_d\|_{0,\Sigma}^2 = (\mathbf{u}_d \cdot \mathbf{n}, \hat{\mathbf{p}}_d)_{0,\Sigma},$$

so that, by using Cauchy-Schwarz and Poincaré-Wirtinger inequalities we get

$$\varepsilon_p \eta \|\hat{\mathbf{p}}_d\|_{1,\Sigma} \le C \|\mathbf{u}_d \cdot \mathbf{n}\|_{0,\Sigma},$$

with $C > 0$. Finally, from (3.3) on $\mathbf{u}_d \cdot \mathbf{n}$ we get estimate (3.4) with $C_2 = C\sqrt{C_1}$.

Then, to estimate the fluid pressure $p$, we test the weak formulation (3.9) with $\hat{\mathbf{q}} = \mathbf{0}$, which yields

$$(p, \operatorname{div} \mathbf{v})_{0,\Omega_d} = 2(\mathbf{D}(\bar{\mathbf{u}}), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} + 2(\mathbf{D}(\mathbf{w}), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} + \frac{\varepsilon_p}{4\kappa} (\bar{\mathbf{u}} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} \\ + (\hat{\mathbf{p}}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma},$$

for all $\mathbf{v} \in \mathcal{U}_d^0$. Therefore, owing to (3.8), we have

$$(p, \operatorname{div} \mathbf{v})_{0,\Omega_d} = 2(\mathbf{D}(\mathbf{u}), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} + \frac{\varepsilon_p}{4\kappa} (\mathbf{u} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} + (\hat{\mathbf{p}}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma}, \quad (3.14)$$

for all $\mathbf{v} \in \mathcal{U}_d^0$. Using once more Theorem 3.3 of 6, as $p \in L^2(\Omega_d)$ there exists a constant $C > 0$ such that there exists $\tilde{\mathbf{v}} \in \mathcal{U}_d^0$ that satisfies both $\operatorname{div} \tilde{\mathbf{v}} = p$ and

$$|\tilde{\mathbf{v}}|_{1,\Omega_d} \le C \|p\|_{0,\Omega_d}. \quad (3.15)$$

Taking $\tilde{\mathbf{v}}$ as a test function in (3.14) gives

$$\|p\|_{0,\Omega_d}^2 = 2(\mathbf{D}(\mathbf{u}), \mathbf{D}(\tilde{\mathbf{v}}))_{0,\Omega_d} + \frac{\varepsilon_p}{4\kappa} (\mathbf{u} \cdot \mathbf{n}, \tilde{\mathbf{v}} \cdot \mathbf{n})_{0,\Sigma} + (\hat{\mathbf{p}}, \tilde{\mathbf{v}} \cdot \mathbf{n})_{0,\Sigma},$$

so that, using Cauchy-Schwarz inequality, the continuity of the trace operator and estimate (3.15) on $\tilde{\mathbf{v}}$, we get

$$\|p\|_{0,\Omega_d}^2 \le C_3 \left( |\mathbf{u}|_{1,\Omega_d} + \frac{\varepsilon_p}{4\kappa} \|\mathbf{u} \cdot \mathbf{n}\|_{0,\Sigma} + \|\hat{\mathbf{p}}\|_{0,\Sigma} \right) \|p\|_{0,\Omega_d},$$

with $C_3 > 0$ a positive constant. Estimate (3.5) then follows after division by $\|p\|_{0,\Omega_d}$ in the previous bound and by using (3.3), which concludes the proof of Theorem 3.□
12 Champion, Fernández, Grandmont, Vergnet, Vidrascu

Estimates (3.3)-(3.5) will be essential for the asymptotic analysis (with respect to $\kappa$ and $\eta$) conducted in Section 3.2. On one hand, we can infer from (3.3) that the fluid velocity $\mathbf{u}_d$ is bounded in $H^1(\Omega_d)$, independently of $\kappa$ and $\eta$. On the other hand, we lose control on the fluid pressure in (3.5) when either $\kappa$ or $\eta$ vanish. The purpose of the next proposition is to show that the zero-mean part of the fluid pressure can be bounded independently of $\kappa$ and $\eta$. We hence decompose the fluid pressure $p_d = p_d^* + c_d$, using the direct sum $L^2 = L_0^2 \oplus \mathbb{R}$, so that

$$c_d = \frac{1}{|\Omega_d|} \int_{\Omega_d} p_d, \quad p_d^* = p_d - c_d \in L_0^2(\Omega_d).$$

**Proposition 3.1.** For any $d \in (0, \tilde{L} - 2)$, there exists a positive constant $C_4 > 0$, independent of the conductivity parameters, such that the zero-mean part $p_d^*$ of the fluid pressure $p_d$, solution of system (2.8), satisfies the a priori estimate

$$\|p_d^*\|_{0,\Omega_d} \le C_4. \qquad (3.16)$$

**Proof.** The proof is very similar to the proof of estimate (3.5). We use the notations introduced in the proof of Theorem 3.1 and, as before, we drop the subscript $d$ of all variables for readability.

We consider a test function $\mathbf{v} \in \mathcal{U}_d^0$ which also satisfies $\mathbf{v} \cdot \mathbf{n} = 0$ on the bottom wall $\Sigma$. Thus, the constant $c$ does not play any role in the weak formulation (3.9) since

$$(c, \operatorname{div} \mathbf{v})_{0,\Omega_d} = c \int_{\Omega_d} \operatorname{div} \mathbf{v} = c \int_{\partial \Omega_d} \mathbf{v} \cdot \mathbf{n} = 0.$$

The weak formulation (3.2) with $\hat{q} = 0$ therefore gives

$$(p^*, \operatorname{div} \mathbf{v})_{0,\Omega_d} = 2(D(\tilde{u}), D(\mathbf{v}))_{0,\Omega_d} \qquad (3.17)$$

for all $\mathbf{v} \in \mathcal{U}_d^0 \cap \{\mathbf{v} \in H^1(\Omega_d)^2, \mathbf{v} \cdot \mathbf{n} = 0 \text{ on } \Sigma\}$. According to Bogovskii's Lemma 2, there exists $C_4 > 0$ such that for all $p^* \in L_0^2(\Omega_d)$ there exists $\tilde{v} \in H_0^1(\Omega_d)$ such that $\operatorname{div} \tilde{v} = p^*$ and

$$|\tilde{v}|_{1,\Omega_d} \le C_4 \|p^*\|_{0,\Omega_d}.$$

Taking $\mathbf{v} = \tilde{v}$ in (3.17), applying the Cauchy-Schwarz inequality and using estimate (3.3) yields

$$\|p^*\|_{0,\Omega_d}^2 = 2(D(\tilde{u}), D(\tilde{v}))_{0,\Omega_d} \le |\tilde{u}|_{1,\Omega_d} |\tilde{v}|_{1,\Omega_d} \le C_4 \|p^*\|_{0,\Omega_d}.$$

This completes the proof. $\square$

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

**Remark 3.1.** Now, to bring insight on the fluid pressure constant $c_d$, let us note that the weak formulation (3.2) gives the following link between $\hat{p}_d$, $p_d^*$ and $c_d$

$$\langle \sigma(\mathbf{u}_d, p_d^* + c_d)\mathbf{n} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n} \rangle_{(H_{00}^{1/2}(\Sigma))', H_{00}^{1/2}(\Sigma)} = - \left( \hat{p}_d + \frac{\varepsilon_P}{4\kappa} \mathbf{u}_d \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n} \right)_{0,\Sigma}$$
for all quad $v \in \mathcal{U}_d^0$. Thanks to the $H^1$ regularity of $\hat{p}_d$ and the $H^{1/2}$ regularity of the trace of $u_d \cdot n$ and using the Hahn-Banach theorem and density arguments, we infer that

$$ \sigma(u_d, p_d^* + c_d)n \cdot n = - \left( \hat{p}_d + \frac{\varepsilon_p}{4\kappa} u_d \cdot n \right) \quad \text{in } L^2(\Sigma). \quad (3.18) $$

We can get rid of the Darcy pressure term by integrating this relation over $\Sigma$ and by using that $\hat{p}_d \in L_0^2(\Omega_d)$ and $\int_{\Sigma} u_d \cdot n = 0$, which yields

$$ \int_{\Sigma} \sigma(u_d, p_d^*) n \cdot n - 2Lc_d = - \int_{\Sigma} \hat{p}_d - \frac{\varepsilon_p}{4\kappa} \int_{\Sigma} u_d \cdot n = 0, $$

so that the pressure constant is given by

$$ c_d = \frac{1}{2L} \int_{\Sigma} \sigma(u_d, p_d^*) n \cdot n. \quad (3.19) $$

Note that it is not obvious to obtain an a priori estimate on $c_d$ with respect to the data of the problem.

## 3.2. Asymptotic analysis

In this section, we use the a priori estimates of Theorem 3.1 to study the asymptotic behavior of system (2.8) with respect to $\kappa$ or $\eta$. These results are stated in Theorem 3.2 and summarized in Table 1 below.

### 3.2.1. Preliminary definitions

We first introduce the different boundary conditions on $\Sigma$ that can be obtained as a limit of problem (2.8). We look for the fluid velocity $u_d: \Omega_d \to \mathbb{R}^2$ and the fluid pressure $p_d: \Omega_d \to \mathbb{R}$ solution of (2.8a)-(2.8d), that we recall here,

$$ \left\{ \begin{array}{ll} -\mathrm{div} \, \sigma(u_d, p_d) = 0 & \text{in } \Omega_d, \\ \mathrm{div} \, u_d = 0 & \text{in } \Omega_d, \\ u_d = 0 & \text{on } \Gamma, \\ u_d = e_2 & \text{on } \partial S, \end{array} \right. $$

supplemented with either one of the following boundary conditions on the bottom wall $\Sigma$:

* Navier boundary conditions:

$$ \left\{ \begin{array}{ll} u_d \cdot n = 0 & \text{on } \Sigma, \\ \sigma(u_d, p_d)n \cdot \tau = 0 & \text{on } \Sigma. \end{array} \right. \quad (3.20) $$

* Neumann boundary conditions:

$$ \sigma(u_d, p_d)n = 0 \quad \text{on } \Sigma. \quad (3.21) $$
14 Champion, Fernández, Grandmont, Vergnet, Vidrascu

• Robin boundary conditions:

$$ \left\{ \begin{array}{l} \sigma(\mathbf{u}_d, p_d)\mathbf{n} \cdot \mathbf{n} = -\frac{\varepsilon_p}{4\kappa}\mathbf{u}_d \cdot \mathbf{n} \quad \text{on } \Sigma, \\ \sigma(\mathbf{u}_d, p_d) \cdot \mathbf{n} \cdot \boldsymbol{\tau} = 0 \quad \text{on } \Sigma. \end{array} \right. \qquad (3.22) $$

• Darcy boundary conditions without correction term

$$ \left\{ \begin{array}{ll} \sigma(\mathbf{u}_d, p_d)\mathbf{n} \cdot \mathbf{n} = -\hat{p}_d & \text{on } \Sigma, \\ \sigma(\mathbf{u}_d, p_d)\mathbf{n} \cdot \boldsymbol{\tau} = 0 & \text{on } \Sigma, \\ -\operatorname{div}_{\boldsymbol{\tau}} (\varepsilon_p \eta \nabla_{\boldsymbol{\tau}} \hat{p}_d) = \mathbf{u}_d \cdot \mathbf{n} & \text{on } \Sigma, \\ \nabla_{\boldsymbol{\tau}} \hat{p}_d \cdot \boldsymbol{\tau} = 0 & \text{on } \partial\Sigma. \end{array} \right. \qquad (3.23) $$

In what follows, we shall denote the weak solutions of system (2.8a)-(2.8d) with boundary conditions of type Navier (3.20), Neumann (3.21), Robin (3.22) or Darcy without correction (3.23) by

$$ \begin{aligned} (\mathbf{u}_d^{Na}, p_d^{Na}) &\in \mathcal{U}_d^0 \times L_0^2(\Omega_d), \\ (\mathbf{u}_d^{Ne}, p_d^{Ne}) &\in \mathcal{U}_d^0 \times L^2(\Omega_d), \\ (\mathbf{u}_d^R, p_d^R) &\in \mathcal{U}_d^0 \times L^2(\Omega_d), \\ (\mathbf{u}_d^{Dwc}, \hat{p}_d^{Dwc}, p_d^{Dwc}) &\in \mathcal{U}_d^0 \times \mathcal{D} \times L^2(\Omega_d), \end{aligned} $$

respectively. The Navier, Neumann and Robin boundary conditions are classical in different fluid-structure interaction settings. Existence and uniqueness of the associated weak solutions can be established using standard arguments. For the coupling with the Darcy without correction (3.23), the well-posedness of the solution can be proved by using the same arguments as for system (2.8) in Section 3.1.

In what follows, we shall make use of the following functional spaces involving Navier boundary conditions

$$ \mathcal{U}_d^{Na} \stackrel{\text{def}}{=} \{ \mathbf{u} \in H^1(\Omega_d)^2 : \mathbf{u} = \mathbf{e}_2 \text{ on } \partial S_d, \mathbf{u} = \mathbf{0} \text{ on } \Gamma, \mathbf{u} \cdot \mathbf{n} = \mathbf{0} \text{ on } \Sigma \}, $$

$$ \mathcal{U}_d^{0,Na} \stackrel{\text{def}}{=} \{ \mathbf{u} \in H^1(\Omega_d)^2 : \mathbf{u} = \mathbf{0} \text{ on } \partial S_d \cup \Gamma, \mathbf{u} \cdot \mathbf{n} = \mathbf{0} \text{ on } \Sigma \} $$

$$ \text{and the divergence-free subspace } V_d \stackrel{\text{def}}{=} \{ \mathbf{u} \in H^1(\Omega_d)^2 : \operatorname{div} \mathbf{u} = \mathbf{0} \}. $$

### 3.2.2. Main result

The following result provides the asymptotic behavior of system (2.8) with respect to the parameters of the Darcy layer.

**Theorem 3.2 (Asymptotic analysis).** Let $d > 0$ be given and let $(\mathbf{u}_d, \hat{p}_d, p_d) \in \mathcal{U}_d \times \mathcal{D} \times L^2(\Omega_d)$ be the solution of (2.8) which depends on the conductivity parameters $\kappa$ and $\eta$. The following propositions hold true:

(i) When either $\kappa$ or $\eta$ goes to zero, $(\mathbf{u}_d, p_d - \frac{1}{|\Omega_d|} \int_{\Omega_d} p_d)$ weakly converges towards $(\mathbf{u}_d^{Na}, p_d^{Na})$ in $H^1(\Omega_d) \times L^2(\Omega_g)$ solution of (2.8a)-(2.8d) with Navier boundary conditions (3.20) on $\Sigma$;
(ii) When both $\kappa$ and $\eta$ go to infinity, $(u_d, p_d)$ weakly converges towards $(u_d^{Ne}, p_d^{Ne})$ in $H^1(\Omega_d)^2 \times L^2(\Omega_g)$ solution of (2.8a)-(2.8d) with Neumann boundary conditions (3.21) on $\Sigma$;

(iii) When $\kappa$ fixed and $\eta \to \infty$, $(u_d, p_d)$ weakly converges towards $(u_d^R, p_d^R)$ in $H^1(\Omega_d)^2 \times L^2(\Omega_g)$ solution of (2.8a)-(2.8d) with Robin boundary conditions (3.22) on $\Sigma$;

(iv) When $\kappa \to \infty$ and $\eta$ fixed, $(u_d, p_d, \hat{p}_d)$ weakly converges towards $(u_d^{Dwc}, p_d^{Dwc}, \hat{p}_d^{Dwc})$ in $H^1(\Omega_d)^2 \times L^2(\Omega_g) \times H^1(\Sigma)$ solution of (2.8a)-(2.8d) coupled with Darcy without correction system (3.23) on $\Sigma$.

**Remark 3.2.** Theorem 3.2 mainly focuses on the asymptotic behavior of the pure fluid part of (2.8). Some asymptotic results for $\hat{p}$ in the limit cases mentioned above are discussed in Remark 3.3. In particular, it is worth mentioning that, in the different cases studied in Theorem 3.2 (except in case (iv)), the fluid unknowns $(u, p)$ and the porous pressure $\hat{p}$ are not coupled anymore.

**Proof.** We first gather the cases which converge towards the Navier boundary conditions ($\kappa$ or $\eta$ goes to 0) and then derive the other cases. This distinction is motivated by estimates (3.3) and (3.5), in which some control on the solution is lost when either $\kappa$ or $\eta$ goes to 0, so that the passage to the limit is more involved than in the other cases.

In this proof, we assume the gap distance $d$ is fixed and we drop the $d$ subscript in all variables. Instead, we introduce the notations $u_{\kappa,\eta}$, $p_{\kappa,\eta}$ and $\hat{p}_{\kappa,\eta}$, that makes the dependence on the conductivity parameters explicit. Thanks to a priori estimate (3.3) on $u_{\kappa,\eta}$, we know there exists a subsequence that converges weakly in $H^1(\Omega_d)^2$. By abuse of notation, we also denote by $u_{\kappa,\eta}$ the subsequence and $u_l$ its limit:

$$u_{\kappa,\eta} \xrightarrow[\kappa \text{ or } \eta \to 0]{H^1} u_l.$$

By continuity of the trace and of the divergence operators, we have $u_l \in V_d \cap U_d$. We now consider the different cases.

*Proof of (i).* When $\kappa$ or $\eta$ goes to 0, to show that $u_l \cdot n = 0$, we have to consider two cases. When $\kappa \to 0$, a priori estimate (3.3) gives straightforwardly $u_l \cdot n = 0$ in $L^2(\Sigma)$. When $\eta \to 0$, taking the weak formulation (3.2) with $v = 0$ yields

$$(u_{\kappa,\eta} \cdot n, \hat{q})_{0,\Sigma} = \sqrt{\varepsilon_p \eta} (\sqrt{\varepsilon_p \eta} \nabla_{\tau} \hat{p}_{\kappa,\eta}, \nabla_{\tau} \hat{q})_{0,\Sigma} \quad \forall \hat{q} \in \mathcal{D}.$$

Owing to (3.3), we have $\sqrt{\varepsilon_p \eta} \|\hat{p}_{\kappa,\eta}\|_{1,\Sigma} \le C$, so that we can pass to limit in the previous identity and obtain

$$(u_l \cdot n, \hat{q})_{0,\Sigma} = 0 \quad \forall \hat{q} \in \mathcal{D}$$

and, by the density of $H^1(\Sigma)$ in $L^2(\Sigma)$, we have

$$(u_l \cdot n, \hat{q})_{0,\Sigma} = 0 \quad \forall \hat{q} \in L_0^2(\Sigma),$$
16 Champion, Fernández, Grandmont, Vergnet, Vidrascu

which implies that $\mathbf{u}_l \cdot \mathbf{n}$ is constant on $\Sigma$. Moreover, since $\mathbf{u}_l \cdot \mathbf{n} \in H_{00}^{\frac{1}{2}}(\Sigma)$, we recover $\mathbf{u}_l \cdot \mathbf{n} = 0$ also in this case. We therefore have that $\mathbf{u}_l \in V_d \cap \mathcal{U}_d^{Na}$ either when $\kappa$ or $\eta$ goes to zero.

Testing (3.2) with $\mathbf{v} \in \mathcal{U}_d^{0,Na} \cap V_d$ and $\hat{q} = 0$, yields

$$ (2\mathbf{D}(\mathbf{u}_{\kappa,\eta}), \mathbf{D}(\mathbf{v})) = 0 \quad \forall \mathbf{v} \in \mathcal{U}_d^{0,Na} \cap V_d, $$

so that, passing to the limit,

$$ (2\mathbf{D}(\mathbf{u}_l), \mathbf{D}(\mathbf{v})) = 0 \quad \forall \mathbf{v} \in \mathcal{U}_d^{0,Na} \cap V_d. $$

By uniqueness of the weak solution associated to the Navier problem (2.8a)-(2.8d) with boundary conditions (3.20), we obtain that $\mathbf{u}_l = \mathbf{u}^{Na}$. Furthermore, thanks to the sequential characterization of the limit, we conclude that

$$ \mathbf{u}_{\kappa, \eta} \xrightarrow[\kappa \text{ or } \eta \to 0]{H^1} \mathbf{u}^{Na}. \quad (3.24) $$

For the pressure, we consider its the zero mean part

$$ p_{\kappa, \eta}^* \stackrel{\text{def}}{=} p_{\kappa, \eta} - \frac{1}{|\Omega_d|} \int_{\Omega_d} p_{\kappa, \eta}, $$

which, according to Proposition 3.1, is bounded independently of $\kappa$ and $\eta$. We then have a subsequence $(\mathbf{u}_{\kappa,\eta}, p_{\kappa,\eta}^*)$ in $\mathcal{U}_d \times L_0^2(\Omega_g)$ which weakly converges towards $(\mathbf{u}^{Na}, p_l^*)$. Testing (3.2) with $\mathbf{v} \in \mathcal{U}_d^{0,Na}$ and $\hat{q} = 0$, yields

$$ 2(\mathbf{D}(\mathbf{u}_{\kappa,\eta}), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} - (p_{\kappa,\eta}^*, \operatorname{div} \mathbf{v})_{0,\Omega_d} = 0 \quad \forall \mathbf{v} \in \mathcal{U}_d^{0,Na}, $$

which, passing to the limit, gives

$$ 2(\mathbf{D}(\mathbf{u}^{Na}), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} - (p_l^*, \operatorname{div} \mathbf{v})_{0,\Omega_d} = 0 \quad \forall \mathbf{v} \in \mathcal{U}_d^{0,Na}. $$

We hence have $p_l^* = p^{Na}$ and

$$ p_{\kappa, \eta}^* \xrightarrow[\kappa \text{ or } \eta \to 0]{L^2} p^{Na}. \quad (3.25) $$

This concludes the proof of (i).

*Proof of (ii) and (iii).* When neither $\eta$ nor $\kappa$ vanish, we have control on $\mathbf{u}_{\kappa,\eta}$, $p_{\kappa,\eta}$ and $\hat{p}_{\kappa,\eta}$ thanks to the estimates (3.3) and (3.5). When $\eta \to \infty$, from (3.3) we get that

$$ \hat{p}_{\kappa, \eta} \xrightarrow[\eta \to \infty]{H^1} 0. $$

So, passing to the limit in the weak formulation (3.2) with $\hat{q} = 0$, gives

$$ (2\mathbf{D}(\mathbf{u}_{\kappa,\eta}), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} + \frac{\varepsilon_D}{4\kappa} (\mathbf{u}_{\kappa,\eta} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} - (p_{\kappa,\eta}, \operatorname{div} \mathbf{v})_{0,\Omega_d} = 0 \quad \forall \mathbf{v} \in \mathcal{U}_d^0. $$

For $\kappa$ fixed, we recover Robin by uniqueness

$$ \mathbf{u}_{\kappa, \eta} \xrightarrow[\kappa \text{ fixed, } \eta \to \infty]{H^1} \mathbf{u}^R, \quad p_{\kappa, \eta} \xrightarrow[\kappa \text{ fixed, } \eta \to \infty]{L^2} p^R. $$
Whenever $\kappa \to \infty$, the term $\frac{\epsilon_p}{4\kappa}(\mathbf{u}_{\kappa,\eta} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma}$ disappears as the bound on $\mathbf{u}_{\kappa,\eta}$ in $H^1(\Omega_d)$ implies that the trace $\mathbf{u}_{\kappa,\eta} \cdot \mathbf{n}$ on $\Sigma$ is bounded in $L^2(\Sigma)$, so that we retrieve the Neumann boundary conditions at the limit and have

$$ \mathbf{u}_{\kappa,\eta} \xrightarrow[\kappa \to \infty, \eta \to \infty]{H^1} \mathbf{u}^{Ne}, \quad p_{\kappa,\eta} \xrightarrow[\kappa \to \infty, \eta \to \infty]{L^2} p^{Ne}. $$

*Proof of (iv).* For $\eta$ fixed and $\kappa \to \infty$, we can pass to the limit in the weak formulation (3.2) and note that only the correction term $\frac{\epsilon_p}{4\kappa}(\mathbf{u}_{\kappa,\eta} \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma}$ vanishes. One can easily prove that this limit problem admits a unique solution corresponding to the Darcy without correction, using the same arguments than for the proof of Theorem 3.1. This concludes the proof. $\square$

<table><caption>Table 1: Asymptotic behavior of the solution $(\mathbf{u}_d, p_d, \hat{p}_d)$ of system (2.8), with respect to the normal and tangential conductivity parameters $\kappa$ and $\eta$. We denote by $p_d^* \stackrel{\text{def}}{=} p_d - \frac{1}{|\Omega_d|} \int_{\Omega_d} p_d$ the zero mean value of $p_d$.</caption><thead><tr><th>η<br>κ</th><th>0</th><th>cst.</th><th>→ ∞</th></tr></thead><tbody><tr><td rowspan="3">0</td><td>Navier<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sup>*</sup><sub>d</sub> <u><sup>L</sup><sub>d</sub></u></td><td>Navier<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sup>*</sup><sub>d</sub> <u><sup>L</sup><sub>d</sub></u><br>p̂<sub>d</sub> <u><sup>H</sup><sub>d</sub></u></td><td>Navier<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sup>*</sup><sub>d</sub> <u><sup>L</sup><sub>d</sub></u><br>p̂<sub>d</sub> <u><sup>H</sup><sub>d</sub></u></td></tr><tr><td>Navier<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sup>*</sup><sub>d</sub> <u><sup>L</sup><sub>d</sub></u></td><td>Darcy</td><td>Robin<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sub>d</sub> <u><sup>L</sup><sub>d</sub></u><br>p̂<sub>d</sub> <u><sup>H</sup><sub>d</sub></u></td></tr><tr><td>Navier<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sup>*</sup><sub>d</sub> <u><sup>L</sup><sub>d</sub></u><br>p̂<sub>d</sub> <u><sup>H</sup><sub>d</sub></u></td><td>Darcy no cor.<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sub>d</sub> <u><sup>L</sup><sub>d</sub></u><br>p̂<sub>d</sub> <u><sup>H</sup><sub>d</sub></u></td><td>Neumann<br><b>u</b><sub>d</sub> <u><sup>H</sup><sub>d</sub></u><br>p<sub>d</sub> <u><sup>L</sup><sub>d</sub></u><br>p̂<sub>d</sub> <u><sup>H</sup><sub>d</sub></u></td></tr></tbody></table>

**Remark 3.3.** The following asymptotic behaviors for the Darcy pressure can also
18 Champion, Fernández, Grandmont, Vergnet, Vidrascu

be proved:

$$ \hat{p}_d \xrightarrow[\eta \to \infty]{H^1} 0, \quad \hat{p}_d \xrightarrow[\kappa \to 0, \eta \text{ fixed}]{H^1} 0, \quad \hat{p}_d \xrightarrow[\frac{\sqrt{\kappa}}{\eta} \to 0]{H^1} 0, $$

$$ \hat{p}_d - \frac{1}{|\Omega_d|} \int_{\Omega_d} p_d \xrightarrow[\kappa \text{ fixed or } \eta \to \infty]{(H_{00}^{1/2})'} -\sigma(\mathbf{u}_d^{Na}, p_d^{Na})\mathbf{n} \cdot \mathbf{n}. $$

Indeed, from (3.3), we have $|\hat{p}_d|_{1,\Sigma}^2 \le \frac{C_1}{\epsilon_p \eta}$, so that $\hat{p}_d$ strongly converges to 0 in $H^1(\Sigma)$ when $\eta \to \infty$. When $\kappa \to 0$ and $\eta$ is fixed, from estimate (3.4) we straightforwardly deduce that

$$ \hat{p}_d \xrightarrow[\kappa \to 0, \eta \text{ fixed}]{H^1} 0. $$

When $\eta \to 0$ and $\kappa$ is fixed or goes to $\infty$, we lose the control on $\hat{p}_d$ provided by (3.3). Nevertheless, from (3.2) with $\hat{q} = 0$, we have

$$ (\hat{p}_d, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} = -2(\mathbf{D}(\mathbf{u}_d), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} - \frac{\epsilon_p}{4\kappa}(\mathbf{u}_d \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Omega_d} + (p_d, \operatorname{div} \mathbf{v})_{0,\Omega_d} $$

for all $v \in \mathcal{U}_d^0$. By splitting the fluid pressure in terms of its zero mean part, $p_d = p_d^* + c_d$, we have

$$ (\hat{p}_d - c_d, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} = -2(\mathbf{D}(\mathbf{u}_d), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} - \frac{\epsilon_p}{4\kappa}(\mathbf{u}_d \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Omega_d} + (p_d^*, \operatorname{div}(\mathbf{v}))_{0,\Omega_d} \quad (3.26) $$

for all $v \in \mathcal{U}_d^0$. As $\kappa$ does not vanish, we can pass to the limit in this expression using the weak convergences (3.24) and (3.25), so that

$$ \lim_{\eta \to 0} (\hat{p}_d - c_d, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} = -2(\mathbf{D}(\mathbf{u}_d^{Na}), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} + (p_d^{Na}, \operatorname{div}(\mathbf{v}))_{0,\Omega_d} \quad \forall \mathbf{v} \in \mathcal{U}_d^0. $$

By density arguments, we conclude that

$$ \hat{p}_d - c_d \xrightarrow[\eta \to 0, \kappa \text{ fixed or } \eta \to \infty]{(H_{00}^{1/2})'} -\sigma(\mathbf{u}_d^{Na}, p_d^{Na})\mathbf{n} \cdot \mathbf{n}. $$

Finally, when $\kappa$ and $\eta$ both go to 0, we cannot pass to the limit in (3.26), but we can conclude in some particular cases. Indeed, from (3.4), we have

$$ \|\hat{p}_d\|_{1,\Sigma} \le C \frac{\sqrt{\kappa}}{\eta}. $$

When $\sqrt{\kappa}$ converges faster than $\eta$ towards 0, the Darcy pressure $\hat{p}_d$ converges strongly in $H^1$ towards 0. If $\eta$ and $\sqrt{\kappa}$ converges at the same speed, $\hat{p}_d$ is bounded in $H^1$ and therefore admits a subsequence weakly convergent in $H^1$ towards an indeterminate limit. Nothing can be said when $\eta$ converges faster than $\sqrt{\kappa}$ towards 0.

We conclude this section by summarizing all the obtained convergence results in Table 1. It should be noted that, in principle, the sole cases which are physically relevant are those where the conductivity parameters are small. It is interesting to
note also that we always recover the solution of the Stokes equations with Navier boundary conditions on the bottom wall when either $\kappa$ or $\eta$ goes to 0. This is in agreement with the following physical intuition: reducing $\kappa$ or $\eta$ makes the fluid seepage through the porous medium more difficult, as $u \cdot n$ converges weakly towards zero. In terms of the contact dynamics of the disk moving towards the bottom wall, this means that in the limit cases where either $\kappa$ or $\eta$ goes to 0, the drag force acting on it is expected to behave like the one with the Navier boundary conditions on the bottom wall. The next section is devoted to the estimation of this drag force.

### 3.3. Drag force

This section is devoted to the main result of the paper: modeling seepage with a porous layer on the bottom wall, as in (2.6), enables contact even with Dirichlet boundary conditions on the disk. The proof of this result is based on an estimate (Theorem 3.3 below) of the asymptotic behavior of the drag force (2.9) acting on the disk as it gets closer to the wall (i.e., $d \to 0$). To this purpose, we build on the approach reported in 14, 9. For Navier or Dirichlet boundary conditions, similar arguments have previously shown that the drag force becomes singular as $d \to 0$ (see, e.g., 14, 22 in 2D and 16, 10 in 3D). In some cases, such as Dirichlet boundary conditions on both the disk and the wall or the combination of Navier and Dirichlet, the singularity in the drag force scales as $d^{-3/2}$ in 2D. This prevents the disk to reach the wall in finite time, which is known as the no-collision paradox.

#### 3.3.1. Main result - Collision result

The next result states that the drag force associated to the scaled fluid problem (2.8) is bounded as $d \to 0$.

**Theorem 3.3 (Drag force estimate).** The drag force $F_d$ given by (2.9) is non-negative and bounded from above independently of $d$. More precisely, we have the following estimate

$$|\mathcal{F}_d| = \mathcal{O}_{d \to 0} \left( 1 + \sum_{i=1}^{6} \frac{1}{(\epsilon_p \eta)^i} + \frac{\epsilon_p}{4\kappa} \sum_{i=2}^{4} \frac{1}{(\epsilon_p \eta)^i} \right). \quad (3.27)$$

Before proceeding with the proof of Theorem 3.3, we can already combine estimate (3.27) with the disk dynamics equation (2.10) to conclude that system (2.6) enables contact between the disk and the bottom wall. This is the purpose of the next corollary.

**Corollary 3.1 (Collision result).** Assume that equation (2.10), with initial conditions $\dot{d}(0) = \dot{d}_0 < 0$ and $d(0) = d_0 > 0$, admits a solution $d \in C^2(0, t_c)$, where $t_c \in \mathbb{R}_+ \cup \{+\infty\}$ denotes the time at which the disk hits the bottom wall. There exists a constant $C > 0$, depending only on $\kappa$ and $\eta$, such that for $\dot{d}_0 < -Cd_0$ we have $t_c < +\infty$ with $d(t_c) = 0$ and $\dot{d}(t_c) < 0$.
20 Champion, Fernández, Grandmont, Vergnet, Vidrascu

**Proof.** We first recall that the disk dynamics associated to the FSI system (2.6) are given by the ordinary differential equation (2.10), namely,

$$ \ddot{d} + \mathcal{F}_d \dot{d} = 0, $$

as long as $d$ remains positive. Owing to Theorem 3.3, there exists a constant $C(\kappa, \eta) > 0$ such that, for all $d > 0$, the drag force satisfies

$$ \mathcal{F}_d \le C(\kappa, \eta). \qquad (3.28) $$

The drag force is also positive $\mathcal{F}_d \ge 0$. This is a direct consequence of the result, shown in the next section (see (3.35) and (3.36)), that the drag force is the minimum of an energy.

Assume that, at time $t = 0$, the disk is located at position $x_0 = d_0 e_2$ and has the initial velocity $\dot{x}_G = \dot{d}_0 e_2$ with $\dot{d}_0 < 0$. In other words, it falls towards the bottom wall. By continuity, the velocity $\dot{d}$ is still negative for a short time interval $0 \le t \le t_1$, so that from (2.10) and (3.28) we get

$$ \ddot{d} = \mathcal{F}_d(-\dot{d}) \quad \implies \quad \ddot{d} \le -C(\kappa, \eta)\dot{d} \qquad (3.29) $$

for all $0 \le t \le t_1$, and thus

$$ \dot{d}(t) \le \dot{d}_0 e^{-Ct}. \qquad (3.30) $$

This implies that, for all $t > 0$ such that $d(t) > 0$, the velocity $\dot{d}(t)$ remains negative and (3.30) holds true. Integrating over time gives for any $t > 0$ such that $d(t) > 0$, we have

$$ d(t) \le d_0 + \dot{d}_0 \int_0^t e^{-Ct} ds \quad \implies \quad d(t) \le d_0 + \underbrace{\frac{\dot{d}_0}{C} (1 - e^{-Ct})}_{\stackrel{\text{def}}{=} f(t)}. $$

The right-hand side function $f$ is obviously decreasing and its limit when $t \to \infty$ is given by

$$ \lim_{t \to \infty} f(t) = d_0 + \frac{\dot{d}_0}{C}. $$

If $\dot{d}_0 < -Cd_0$, there exists a time $t_c > 0$ such that $f(t_c) = 0$, namely,

$$ t_c = - \frac{\ln \left( 1 + C(\kappa, \eta) \frac{d_0}{\dot{d}_0} \right)}{C(\kappa, \eta)}. \qquad (3.31) $$

By continuity, we show that the disk touches the wall at a time $t^* \le t_c$, with $\dot{d}(t^*) < 0$. This completes the proof. □

**Remark 3.4.** Note that in the current model nothing prevents the disk of penetrating the bottom wall. Indeed, if the initial velocity $|\dot{d}_0|$ is large enough, the distance $d$ can vanish with a negative velocity. To avoid this, one often adds a non-penetration constraint in the fluid-structure interaction model (2.6). Usually, this
constraint is imposed by duality with a Lagrange multiplier $\lambda$ that represents the upward force acting on the disk at contact. The resulting problem writes

$$ \begin{cases} \ddot{d} + \mathcal{F}_d \dot{d} = \lambda, \\ d \ge 0, \quad \lambda \ge 0, \quad \lambda d = 0, \end{cases} $$

where the last three conditions represent the so-called non-penetration, compressibility and complementary conditions.

The next sections are devoted to the proof of Theorem 3.3. The first idea consists in expressing the drag force $\mathcal{F}_d$ as the minimum of the energy of the system. This is the purpose of Section 3.3.2. Then, in Section 3.3.3, we construct an admissible minimizer which catches the contact dynamics. Finally, in Section 3.3.4, we estimate the energy associated to the minimizer which yields the bound (3.27). Since the proof builds on 14, we put particular emphasis on the main difficulties arising in the present model with the porous layer.

### 3.3.2. Energy minimization problem

Let $\mathcal{A}_d$ be the admissible function space for the fluid velocity defined by

$$ \mathcal{A}_d \stackrel{\text{def}}{=} \mathcal{U}_d \cap V_d = \{ \mathbf{u} \in H_\Gamma^1 (\Omega_d)^2 : \ \text{div} \, \mathbf{u} = 0 \text{ in } \Omega_d, \ \mathbf{u} = \mathbf{e}_2 \text{ on } \partial S_d \}. \quad (3.32) $$

In this section, we express the Darcy pressure $\hat{p}$ in terms of the fluid velocity via the following operator.

**Definition 3.1 (Operator A).** For any $f \in L_0^2(\Sigma)$ the problem

$$ \text{Find } \hat{p} \in \mathcal{D} \text{ such that } \quad \varepsilon_p \eta \int_{\Sigma} \nabla_{\tau} \hat{p} \cdot \nabla_{\tau} \hat{q} = \int_{\Sigma} f \hat{q}, \quad \forall \hat{q} \in \mathcal{D}, \quad (3.33) $$

is well-posed. Therefore, we denote by $A: L_0^2(\Sigma) \to \mathcal{D}$ the operator that associates, to any $f$ in $L_0^2(\Sigma)$, the unique solution $\hat{p}$ of problem (3.33): $A f = \hat{p}$.

In particular, for all $v \in \mathcal{A}_d$, $v \cdot n|_{\Sigma} \in L_0^2(\Sigma)$ as $v$ is divergence free. We have that $A(v \cdot n)$ is the solution of the following strong form of problem (3.33) with $f = v \cdot n$:

$$ \begin{cases} -\text{div}_{\tau}(\varepsilon_p \eta \nabla_{\tau} \hat{p}) = v \cdot n \quad \text{on } \Sigma, & (3.34a) \\ \nabla_{\tau} \hat{p}(-L) = 0, \quad \nabla_{\tau} \hat{p}(L) = 0. & (3.34b) \end{cases} $$

The following proposition is the cornerstone of our approach to prove Theorem 3.3.

**Proposition 3.2 (Energy minimization problem).** We introduce the energy functional $\mathcal{E}_d: \mathcal{A}_d \to \mathbb{R}^+$ defined by

$$ \mathcal{E}_d(\mathbf{u}) \stackrel{\text{def}}{=} 2 \int_{\Omega_d} \|\mathbf{D}(\mathbf{u})\|^2 + \varepsilon_p \eta \int_{\Sigma} |\nabla_{\tau} A(\mathbf{u} \cdot \mathbf{n})|^2 + \frac{\varepsilon_p}{4\kappa} \int_{\Sigma} |\mathbf{u} \cdot \mathbf{n}|^2 \quad (3.35) $$
22 Champion, Fernández, Grandmont, Vergnet, Vidrascu

for all $\mathbf{u} \in \mathcal{A}_d$. The drag force of (2.8) satisfies

$$ \mathcal{F}_d = \min_{\mathbf{u} \in \mathcal{A}_d} \mathcal{E}_d(\mathbf{u}). \qquad (3.36) $$

Moreover, the energy minimization problem (3.36) admits a unique solution in $\mathcal{A}_d$ which is precisely $\mathbf{u}_d$, the solution of (3.2).

**Proof.** We first prove that

$$ \mathcal{F}_d = \mathcal{E}_d(\mathbf{u}_d). \qquad (3.37) $$

Owing to the weak formulation (3.2), we have

$$ \mathrm{div}(\sigma(\mathbf{u}_d, p_d)) = 0 \quad \text{in} \quad \mathcal{D}'(\Omega_d), $$

which by a standard density argument and thanks to Hahn-Banach Theorem leads to

$$ \mathrm{div}(\sigma(\mathbf{u}_d, p_d)) = 0 \quad \text{in} \quad L^2(\Omega_d). $$

Since $\mathrm{div}(\sigma(\mathbf{u}_d, p_d)) \in L^2(\Omega_d)$, the trace of $\sigma(\mathbf{u}_d, p_d)\mathbf{n}$ can be defined by duality as follows:

$$ \langle \sigma(\mathbf{u}_d, p_d)\mathbf{n}, \xi \rangle_{H^{-\frac{1}{2}}(\partial\Omega_d), H^{\frac{1}{2}}(\partial\Omega_d)} = 2(\mathbf{D}(\mathbf{u}_d), \mathbf{D}(L\xi))_{0,\Omega_d} - (p_d, \mathrm{div}(L\xi))_{0,\Omega_d} $$

for all $\xi \in H^{\frac{1}{2}}(\partial\Omega_d)$, where the symbol L denotes a standard lifting operator from $H^{\frac{1}{2}}(\Sigma)^2$ to $H^1(\Omega_d)^2$. Taking $\xi = \mathbf{u}_d|_{\Sigma}$ and $L\xi = \mathbf{u}_d$, and using that $\mathbf{u}_d \in \mathcal{U}_d \cap V_d$ yields

$$ \langle \sigma(\mathbf{u}_d, p_d)\mathbf{n}, \mathbf{u}_d \rangle_{H^{-\frac{1}{2}}(\Sigma \cup \Gamma), H^{\frac{1}{2}}(\Sigma \cup \Gamma)} + \langle \sigma(\mathbf{u}_d, p_d)\mathbf{n}, e_2 \rangle_{H^{-\frac{1}{2}}(\partial S_d), H^{\frac{1}{2}}(\partial S_d)} = 2 \| \mathbf{D}(\mathbf{u}_d) \|_{0, \Omega_d}^2, \quad (3.38) $$

where drag force appears in the term on $\partial S_d$. It is standard to check that $\sigma(\mathbf{u}_d, p)\mathbf{n} \in L^2(\partial S_d)$. Indeed, the $C^\infty$ regularity of the disk and the Dirichlet boundary condition on $\partial S_d$ guarantee $H^2$ regularity for the fluid velocity $\mathbf{u}_d$ and $H^1$ regularity for the pressure $p_d$ near the ball (see, e.g., 3), so that we have

$$ \mathcal{F}_d = \langle \sigma(\mathbf{u}_d, p_d)\mathbf{n}, e_2 \rangle_{H^{-\frac{1}{2}}(\partial S_d), H^{\frac{1}{2}}(\partial S_d)} = \int_{\partial S_d} \sigma(\mathbf{u}_d, p_d)\mathbf{n} \cdot \mathbf{e}_2. \quad (3.39) $$

As regards the first term (3.38), we use the relation (3.18) and the fact that $\mathbf{u}_d|_{\Gamma} = \mathbf{0}$ to get

$$ \langle \sigma(\mathbf{u}_d, p_d)\mathbf{n}, \mathbf{u}_d \rangle_{H^{-\frac{1}{2}}(\Sigma \cup \Gamma), H^{\frac{1}{2}}(\Sigma \cup \Gamma)} = \langle \sigma(\mathbf{u}_d, p_d)\mathbf{n}, \mathbf{u}_d \rangle_{(H_{00}^{\frac{1}{2}}(\Sigma))', H_{00}^{\frac{1}{2}}(\Sigma)} \\ = -(\hat{p}_d, \mathbf{u}_d \cdot \mathbf{n})_{0,\Sigma} - \frac{\varepsilon_p}{4\kappa} (\mathbf{u}_d \cdot \mathbf{n}, \mathbf{u}_d \cdot \mathbf{n})_{0,\Sigma}. $$

Since by construction $\hat{p}_d = A(\mathbf{u}_d \cdot \mathbf{n})$, taking $\hat{q} = A(\mathbf{u}_d \cdot \mathbf{n})$ in (3.33) yields

$$ (\hat{p}_d, \mathbf{u}_d \cdot \mathbf{n})_{0,\Sigma} = (A(\mathbf{u}_d \cdot \mathbf{n}), \mathbf{u}_d \cdot \mathbf{n})_{0,\Sigma} = \varepsilon_p \eta \int_{\Sigma} |\nabla_{\tau} A(\mathbf{u}_d \cdot \mathbf{n})|^2, $$
so that

$$ \langle \sigma(\mathbf{u}_d, p_d)\mathbf{n}, \mathbf{u}_d \rangle_{H^{-\frac{1}{2}}(\Gamma \cup \Sigma), H^{\frac{1}{2}}(\Gamma \cup \Sigma)} = -\varepsilon_p \eta \int_{\Sigma} |\nabla_{\tau} A(\mathbf{u}_d \cdot \mathbf{n})|^2 - \frac{\varepsilon_p}{4\kappa} \int_{\Sigma} |\mathbf{u}_d \cdot \mathbf{n}|^2. \quad (3.40) $$

Finally, by inserting (3.39) and (3.40) into (3.38), we obtain

$$ \begin{aligned} \mathcal{F}_d &= \int_{\partial S_d} \sigma(\mathbf{u}_d, p_d)\mathbf{n} \cdot e_2 \\ &= \int_{\Omega_d} 2\|\mathbf{D}(\mathbf{u}_d)\|^2 + \eta \varepsilon_p \int_{\Sigma} |\nabla_{\tau} A(\mathbf{u}_d \cdot \mathbf{n})|^2 + \frac{\varepsilon_p}{4\kappa} \int_{\Sigma} |\mathbf{u}_d \cdot \mathbf{n}|^2 \\ &= \mathcal{E}_d(\mathbf{u}_d), \end{aligned} $$

which yields (3.37).

In order to show that $\mathbf{u}_d$ is a minimizer of the functional $\mathcal{E}_d$, we introduce the tangent Hilbert space to $\mathcal{A}_d$:

$$ \mathcal{A}_d^0 \stackrel{\text{def}}{=} \{ \mathbf{u} \in H^1(\Omega_d)^2, \ \text{div} \, \mathbf{u} = 0 \ \text{in} \ \Omega_d, \ \mathbf{u} = 0 \ \text{on} \ \partial S_d \cup \Gamma \}. \quad (3.41) $$

Straightforward computations give

$$ \begin{aligned} \mathcal{E}_d(\mathbf{u}_d + \mathbf{v}) &= \mathcal{E}_d(\mathbf{u}_d) + 2(\mathbf{D}(\mathbf{u}_d), \mathbf{D}(\mathbf{v}))_{0,\Omega_d} + \frac{\varepsilon_p}{4\kappa} (\mathbf{u}_d \cdot \mathbf{n}, \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} \\ &\quad + \varepsilon_p \eta (\nabla_{\tau} A(\mathbf{u}_d \cdot \mathbf{n}), \nabla_{\tau} A(\mathbf{v} \cdot \mathbf{n}))_{0,\Sigma} + (A(\mathbf{u}_d \cdot \mathbf{n}), \mathbf{v} \cdot \mathbf{n})_{0,\Sigma} \\ &\quad - (\mathbf{u}_d \cdot \mathbf{n}, A(\mathbf{v} \cdot \mathbf{n}))_{0,\Sigma} + \mathcal{E}_d(\mathbf{v}). \end{aligned} $$

for all $\mathbf{v} \in \mathcal{A}_d^0$. Testing (3.2) with $\mathbf{v} \in \mathcal{A}_d^0 \subset \mathcal{U}_d^0$, $q = 0$ and $\hat{q} = A(\mathbf{v} \cdot \mathbf{n})$ and since $\hat{p}_d = A(\mathbf{u}_d \cdot \mathbf{n})$, we finally have

$$ \mathcal{E}_d(\mathbf{u}_d + \mathbf{v}) = \mathcal{E}_d(\mathbf{u}_d) + \mathcal{E}_d(\mathbf{v}) \quad \forall \mathbf{v} \in \mathcal{A}_d^0, $$

which, since $\mathcal{E}_d(\mathbf{v}) \ge 0$, yields (3.36).

Finally, to prove the uniqueness of the minimizer, we can use the equivalence

$$ \min_{\mathbf{v} \in \mathcal{A}_d} \mathcal{E}_d(\mathbf{v}) \Leftrightarrow \min_{\mathbf{v} \in \mathcal{A}_d^0} J_d(\mathbf{v}), \quad J_d(\mathbf{v}) \stackrel{\text{def}}{=} \mathcal{E}_d(\mathbf{u}_d + \mathbf{v}) $$

and evaluate the second differential of $J_d$ to show that

$$ D^2 J_d(\mathbf{w})(\mathbf{v}, \mathbf{v}) = a((\mathbf{v}, A(\mathbf{v} \cdot \mathbf{n})), (\mathbf{v}, A(\mathbf{v} \cdot \mathbf{n}))) > \alpha |\mathbf{v}|_{1,\Omega_d}^2, $$

for all non zero $\mathbf{v} \in \mathcal{A}_d^0$ and $\mathbf{w} \in \mathcal{A}_d^0$, where the last inequality follows from the coercivity of the bilinear form $a$ given in (3.10). This guarantees the strict convexity of $J_d$, which implies the uniqueness of the solution of the minimization problem (3.36) and, hence, completes the proof. □

In order to ease the computations, we rewrite the minimization problem, in an equivalent form, in terms of stream functions. This is the purpose of the next result.
24 Champion, Fernández, Grandmont, Vergnet, Vidrascu

**Lemma 3.1 (Stream function).** The admissible space $\mathcal{A}_d$, defined in (3.32), is also given by

$$ \mathcal{A}_d = \left\{ \mathbf{u} \in H^1(\Omega_d) : \exists! \varphi \in H^2(\Omega_d), \partial_1 \varphi|_{\partial S} = 1, \partial_2 \varphi|_{\partial S} = 0, \partial_1 \varphi|_{\Gamma} = 0, \partial_2 \varphi|_{\Gamma} = 0, \varphi(-L, 0) = 0 : \mathbf{u} = -\partial_2 \varphi e_1 + \partial_1 \varphi e_2 \right\}. \quad (3.42) $$

Moreover, we have

$$ \mathcal{E}_d(\mathbf{u}) = \mathcal{E}_d(\mathbf{u}(\varphi)) = \int_{\Omega_d} \left( 4|\partial_1 \partial_2 \varphi|^2 + |\partial_1^2 \varphi|^2 + |\partial_2^2 \varphi|^2 - 2|\partial_2^2 \varphi \partial_1^2 \varphi| \right) + \frac{1}{\varepsilon_p \eta} \int_{\Sigma} |\varphi|^2 + \frac{\varepsilon_p}{4\kappa} \int_{\Sigma} |\partial_1 \varphi|^2 \quad (3.43) $$

for all $\mathbf{u} \in \mathcal{A}_d$.

**Proof.** We start by proving (3.42). Given $\mathbf{u} \in \mathcal{A}_d$, $\mathbf{u}$ is a 2D divergence-free vector field so there exist a streamfunction $\varphi \in H^2(\Omega_d)$, unique up to an additive constant, such that $\mathbf{u} = -\partial_2 \varphi e_1 + \partial_1 \varphi e_2$ (see, e.g., Theorem 3.1 of 11). In (3.42), this constant is fixed by imposing $\varphi(-L, 0) = 0$. The reasons of this choice will be made clear later in the proof. One can easily check that

$$ \mathbf{u}|_{\partial S_d} = e_2, \quad \mathbf{u}|_{\Gamma} = 0, $$

enforce

$$ \partial_1 \varphi|_{\partial S_d} = 1, \quad \partial_2 \varphi|_{\partial S_d} = 0, \quad \partial_1 \varphi|_{\Gamma} = 0, \quad \partial_2 \varphi|_{\Gamma} = 0. \quad (3.44) $$

Conversely, given $\varphi \in H^2(\Omega_d)$, if $\varphi$ satisfies (3.44) then obviously $\mathbf{u} = -\partial_2 \varphi e_1 + \partial_1 \varphi e_2$ belongs to $\mathcal{A}_d$.

In order to get (3.43), we simply replace $\mathbf{u}$ by its associated stream-function $\varphi$ in (3.35). Computations are straightforward for the bulk term $\mathbf{D}(\mathbf{u})$ and the last boundary term. For the term involving $\nabla_{\tau} \mathbf{A}(\mathbf{u} \cdot \mathbf{n})$ we proceed as follows. The Darcy pressure $\hat{\mathbf{p}} = \mathbf{A}(\mathbf{u} \cdot \mathbf{n})$ satisfies (3.34) which gives

$$ \left\{ \begin{array}{l} \mathrm{div}_{\tau} (\varepsilon_p \eta \nabla_{\tau} \hat{\mathbf{p}}(x)) = \mathbf{u}(x, 0) \cdot \mathbf{e}_2 = \partial_1 \varphi(x, 0) \quad \text{in } (-L, L), \\ \nabla_{\tau} \hat{\mathbf{p}}(-L) = 0. \end{array} \right. \quad (3.45) $$

Integrating the first relation yields

$$ \varepsilon_p \eta (\nabla_{\tau} \hat{\mathbf{p}}(x) - \nabla_{\tau} \hat{\mathbf{p}}(-L)) = \int_{-L}^{x} \partial_1 \varphi(s, 0) ds = \varphi(x, 0) - \varphi(-L, 0) \quad \text{in } (-L, L). $$

On the other hand, using (3.45)$_2$, we have

$$ \varepsilon_p \eta \nabla_{\tau} \hat{\mathbf{p}} = \varphi - \varphi(-L, 0) \quad \text{in } \Sigma. \quad (3.46) $$

As a result, since $\varphi(-L, 0) = 0$, we get

$$ \varepsilon_p \eta \int_{\Sigma} |\nabla_{\tau} \mathbf{A}(\mathbf{u} \cdot \mathbf{n})|^2 = \varepsilon_p \eta \int_{\Sigma} \frac{|\varphi|^2}{(\varepsilon_p \eta)^2} = \frac{1}{\varepsilon_p \eta} \int_{\Sigma} |\varphi|^2, $$

which completes the proof. □
### 3.3.3. Focus on contact dynamics - Relaxed energy minimization problem

In the section, the energy functional is split into two contributions: one part leads to a relaxed problem, that we can solve analytically, while the remaining terms will be bounded in Section 3.3.4. In order to both catch the contact dynamics and to be able to simplify the computations, two main choices are made to relax problem (3.36).

First, we limit the analysis to the following region under the disk

$$ \Omega_{d,\delta}^{\text{def}} = \{(x, y) \in \Omega_d : |x| < \delta, \ 0 < y < d + \gamma(x)\}, \quad (3.47) $$

for all $d > 0$ and $0 < \delta < 1$, and where $\gamma: x \in [-\delta, \delta] \mapsto 1 - \sqrt{1 - x^2}$ is a parametrization of a subset of the boundary of the disk (see Figure 2). We also consider the notations

$$ \partial S_{d,\delta}^{\text{def}} = \partial S_d \cap \bar{\Omega}_{d,\delta} = \{(x, \gamma(x)), x \in (-\delta, \delta)\}, \quad \Sigma_{\delta}^{\text{def}} = \Sigma \cap \bar{\Omega}_{d,\delta} = (-\delta, \delta). $$

Second, we introduce the relaxed energy functional $\tilde{\mathcal{E}}_d$ given by

<figure><img src="image_4.png" /><figcaption>Figure 2: Contact zone $\Omega_{d,\delta}$.</figcaption></figure>

$$ \tilde{\mathcal{E}}_d(\mathbf{u}) \stackrel{\text{def}}{=} \int_{\Omega_{d,\delta}} |\partial_2 u_1|^2 + \varepsilon_p \eta \int_{-\delta}^{\delta} |\nabla_{\tau}(A u_2)|^2 = \int_{\Omega_{d,\delta}} |\partial_2^2 \varphi|^2 + \frac{1}{\varepsilon_p \eta} \int_{-\delta}^{\delta} |\varphi(x, 0)|^2. \quad (3.48) $$

Note that only the contribution $\partial_2 u_1 = \partial_2^2 \varphi$ has been kept in the volumetric energy term. This is motivated by physical intuition and previous studies (see, e.g, 14, 22, 9) that show that, in the case of Dirichlet or Navier boundary conditions, the term $\partial_2 u_1 = \partial_2^2 \varphi$ causes the no-collision paradox. Indeed, high velocity gradients arise when the fluid has to escape tangentially in the narrow vertical space between the disk and the bottom wall. The boundary term in $u_2 = \partial_1 \varphi$ has been removed to avoid having to deal with terms on the side boundaries $\partial \Omega_{d,\delta} \setminus (\partial S_{d,\delta} \cup \Sigma_{\delta})$. It should be noted that, all the left out terms will be estimated in Section 3.3.
Similarly, we introduce the relaxed admissible function space $\tilde{\mathcal{A}}_d$ as the restriction of functions of $\mathcal{A}_d$ to the contact zone $\Omega_{d,\delta}$, viz.,

$$ \tilde{\mathcal{A}}_d \stackrel{\text{def}}{=} \{ \mathbf{u}|_{\Omega_{d,\delta}} : \mathbf{u} \in \mathcal{A}_d \}. \qquad (3.49) $$

In order to give a precise characterization of $\tilde{\mathcal{A}}_d$, we need the following two technical lemmas.

<figure><img src="image_5.png" /><figcaption>Figure 3: Symmetries of the problem.</figcaption></figure>

**Lemma 3.2 (Symmetry of the solution).** *Solution of (2.8) satisfies the following symmetries:*

$$ \begin{aligned} u_1(-x, y) &= -u_1(x, y), \\ u_2(-x, y) &= u_2(x, y), \\ p(-x, y) &= p(x, y), \\ \hat{p}_d(-x) &= \hat{p}_d(x). \end{aligned} $$

**Proof.** The proof is straightforward thanks to the symmetries of the problem (see Figure 3) and the uniqueness of the solution of (2.8). □

The next lemma gives the characterization of $\tilde{\mathcal{A}}_d$.

**Lemma 3.3.** Let $\tilde{\mathcal{A}}_d$ be the relaxed space given by (3.49). There holds

$$ \tilde{\mathcal{A}}_d = \left\{ \mathbf{u} \in H^1(\Omega_{d,\delta}) : \exists! \varphi \in H^2(\Omega_{d,\delta}), \quad \varphi|_{\partial S_{d,\delta}} = x, \quad \partial_2 \varphi|_{\partial S_{d,\delta}} = 0 : \mathbf{u} = -\partial_2 \varphi e_1 + \partial_1 \varphi e_2, \right\}. \quad (3.50) $$

**Proof.** Let $\mathbf{u} = -\partial_2 \varphi e_1 + \partial_1 \varphi e_2 \in \mathcal{A}_d$. The boundary conditions $\mathbf{u}|_{\Sigma} = e_2$ implies

$$ \partial_1 \varphi|_{\partial S_{d,\delta}} = 1, \quad \partial_2 \varphi|_{\partial S_{d,\delta}} = 0, $$
from which we deduce

$$ \frac{d}{dx}\varphi(x, \gamma(x)) = \partial_1\varphi(x, \gamma(x)) + \gamma'(x)\partial_2\varphi(x, \gamma(x)) = 1 \quad \text{in} \quad (-\delta, \delta), \\ \partial_2\varphi|_{\partial S_{d,\delta}} = 0, $$

so that

$$ \varphi|_{\partial S_{d,\delta}} = x + c, \quad \partial_2\varphi|_{\partial S_{d,\delta}} = 0. $$

As we already fixed the constant by imposing $\varphi(0, -L) = 0$ in (3.42), the constant $c$ cannot be arbitrarily chosen. We have $c = \varphi(0, d)$. Owing to Lemma 3.2, by a symmetry argument on $u$, we have $u_1(0, y) = 0$ in $(0, d)$, so that

$$ 0 = \int_{0}^{d} \partial_{2}\varphi(0, y)dy \Rightarrow \varphi(0, d) = \varphi(0, 0) = c. $$

As $\hat{p}$ is also symmetric and thanks to (3.46), we have

$$ 0 = \nabla_{\tau} \hat{p}(0) = \frac{1}{\varepsilon_p \eta} (\varphi(0, 0) - \varphi(0, -L)) = \frac{1}{\varepsilon_p \eta} \varphi(0, 0), $$

which gives $c = 0$. As a result, we finally recover the following boundary conditions on $\partial S_{d,\delta}$

$$ \varphi|_{\partial S_{d,\delta}} = x, \quad \partial_2\varphi|_{\partial S_{d,\delta}} = 0. $$

For the boundary conditions on $\partial\Omega_{d,\delta} \setminus (\Sigma_{\delta} \cup \partial S_{d,\delta})$, we remark that for all $\delta < 1$ the following restriction function on $\Omega_{d,\delta}$

$$ \{\varphi \in H^2(\Omega_d), \partial_1\varphi|_{\Gamma} = \partial_2\varphi|_{\Gamma} = 0\} \rightarrow H^2(\Omega_{d,\delta}) \\ \varphi \mapsto \varphi|_{\Omega_{d,\delta}} $$

is surjective (see for example the extension $\hat{\varphi}$ in Section 3.3.4) that is why the boundary conditions on $\partial\Omega_{d,\delta} \setminus (\Sigma_{\delta} \cup \partial S_{d,\delta})$ do not appear in the relaxed space definition (3.50). This completes the proof. $\square$

The relaxed minimization problem reads therefore as follows: Find $\tilde{u} = -\partial_2\tilde{\varphi}e_1 + \partial_1\tilde{\varphi}e_2 \in \tilde{\mathcal{A}}_d$ such that

$$ \tilde{\mathcal{E}}_d(\tilde{u}) = \min_{\tilde{v} \in \tilde{\mathcal{A}}_d} \tilde{\mathcal{E}}_d(\tilde{v}). \qquad (3.51) $$

We finally define the relaxed drag force as the energy minimum of this relaxed energy, namely,

$$ \tilde{\mathcal{F}}_d \stackrel{\text{def}}{=} \tilde{\mathcal{E}}_d(\tilde{u}). \qquad (3.52) $$

The next proposition solves explicitly the relaxed minimization problem (3.51).

**Proposition 3.3 (Explicit expression of $\tilde{\varphi}$).** The relaxed minimization problem (3.51) admits a unique solution $\tilde{u} = -\partial_2\tilde{\varphi}e_1 + \partial_1\tilde{\varphi}e_2 \in \tilde{\mathcal{A}}_d$, given by

$$ \tilde{\varphi}(x, y) = x\phi\left(x, \frac{y}{d + \gamma(x)}\right), \qquad (3.53) $$
28 Champion, Fernández, Grandmont, Vergnet, Vidrascu

with

$$ \alpha_d(x) \stackrel{\text{def}}{=} \frac{(d + \gamma(x))^3}{\varepsilon_p \eta}, \quad \phi(x, t) \stackrel{\text{def}}{=} -\frac{\alpha_d(x)t^3}{6 + 2\alpha_d(x)} + \frac{3\alpha_d(x)t}{6 + 2\alpha_d(x)} + \frac{6}{6 + 2\alpha_d(x)}. \quad (3.54) $$

**Proof.** As in the proof of Proposition 3.2, we can rewrite (3.51) as a minimization problem in the tangent Hilbert space $\tilde{\mathcal{A}}_d^0$:

$$ \tilde{\mathcal{A}}_d^0 \stackrel{\text{def}}{=} \left\{ \mathbf{u} \in H^1(\Omega_{d,\delta}), \quad \mathbf{u} = -\partial_2 \varphi \mathbf{e}_1 + \partial_1 \varphi \mathbf{e}_2, \right. \\ \left. \exists! \varphi \in H^2(\Omega_d), \quad \varphi|_{\partial S_{d,\delta}} = 0, \quad \partial_2 \varphi|_{\partial S_{d,\delta}} = 0 \right\}. $$

So, let $\mathbf{u}_0 \in \tilde{\mathcal{A}}_d$ be given, instead of (3.51), we consider the auxiliary equivalent problem

$$ \min_{\mathbf{v} \in \tilde{\mathcal{A}}_d^0} \tilde{J}_d(\mathbf{v}), \quad \tilde{J}_d(\mathbf{v}) \stackrel{\text{def}}{=} \tilde{\mathcal{E}}_d(\mathbf{v} + \mathbf{u}_0). \quad (3.55) $$

For all $\mathbf{u} \in \tilde{\mathcal{A}}_d^0$, we have

$$ D^2 \tilde{J}_d(\mathbf{u})(\mathbf{v}, \mathbf{v}) = 2 \left( \int_{\Omega_{d,\delta}} |\partial_2 v_1|^2 + \varepsilon_p \eta \int_{-\delta}^{\delta} |\nabla_{\tau}(A v_2)|^2 \right) \ge 0 \quad \forall \mathbf{v} \in \tilde{\mathcal{A}}_d^0. $$

Let $\mathbf{v} \in \tilde{\mathcal{A}}_d^0$ and $\psi$ its associated streamfunction. The relation $D^2 \tilde{J}_d(\mathbf{u})(\mathbf{v}, \mathbf{v}) = 0$ implies that

$$ \left\{ \begin{array}{ll} \partial_2^2 \psi = 0 & \text{in } \Omega_{d,\delta}, \\ \psi = 0 & \text{on } \Sigma_{\delta}. \end{array} \right. $$

Therefore, there exists a function $g$ such that $\psi(x, y) = g(x)y$ in $\Omega_{d,\delta}$. Since $\psi(x, d+\gamma(x)) = 0$, we conclude that $\psi = 0$, so that $\mathbf{v} = \mathbf{0}$. As a result, $D^2 \tilde{J}_d(\mathbf{u})(\mathbf{v}, \mathbf{v}) > 0$ for all $\mathbf{v} \in \tilde{\mathcal{A}}_d^0 \setminus \{\mathbf{0}\}$. This guarantees the strict convexity of $\tilde{J}_d$, so that, if a minimum exists for problem (3.55), it is unique.

The minimum $\tilde{\mathbf{u}} \in \tilde{\mathcal{A}}_d$ of (3.51) can be characterized by the Euler equations associated to the auxiliary minimization problem (3.55), which writes $J'_d(\tilde{\mathbf{u}} - \mathbf{u}_0)\mathbf{v} = 0$ for all $\mathbf{v} \in \tilde{\mathcal{A}}_d^0$. After integration by parts, we obtain

$$ \int_{\Omega_{d,\delta}} \partial_2^4 \tilde{\varphi} \cdot \psi + \int_{-\delta}^{\delta} \left( \partial_2^3 \tilde{\varphi} + \frac{1}{\varepsilon_p \eta} \tilde{\varphi} \right) \cdot \psi - \int_{-\delta}^{\delta} \partial_2^2 \tilde{\varphi} \cdot \partial_2 \psi = 0 \quad \forall \psi \in \mathcal{B}, \quad (3.56) $$

where $\mathcal{B}$ simply denotes the stream-function space associated to $\tilde{\mathcal{A}}_d^0$, namely,

$$ \mathcal{B} = \left\{ \varphi \in H^2(\Omega_{d,\delta}), \quad \varphi|_{\partial S_{d,\delta}} = 0, \quad \partial_2 \varphi|_{\partial S_{d,\delta}} = 0 \right\}. $$
From (3.50) and (3.56) we obtain, using standard density and trace arguments, the
following PDE system for $\tilde{\varphi}$

$$
\left\{
\begin{array}{ll}
\partial_2^4 \tilde{\varphi} = 0 & \text{in } \Omega_{d,\delta}, \\
\partial_2^3 \tilde{\varphi} + \frac{1}{\epsilon_p \eta} \tilde{\varphi} = 0 & \text{on } \Sigma_{\delta}, \\
\partial_2^2 \tilde{\varphi} = 0 & \text{on } \Sigma_{\delta}, \\
\partial_2 \tilde{\varphi} = 0 & \text{on } \partial S_{d,\delta}, \\
\tilde{\varphi} = x & \text{on } \partial S_{d,\delta}.
\end{array}
\right.
$$

whose solution is given by (3.53), which characterizes the unique minimizer $\tilde{u}$ of (3.51). This completes the proof. □

The next theorem shows that the relaxed drag force $\tilde{\mathcal{F}}_d$ is bounded irrespectively of $d$.

**Theorem 3.4.** The relaxed drag force (3.52) is given by

$$
\tilde{\mathcal{F}}_d = 12 \int_{-\delta}^{\delta} \frac{\alpha_d(x)^2 + 3\alpha_d(x)}{(6 + 2\alpha_d(x))^2} \frac{x^2}{(d + \gamma(x))^3} dx. \quad (3.57)
$$

In particular, we have that

$$
\tilde{\mathcal{F}}_d = O_{d \to 0} \left( \delta^3 \left( \frac{1}{(\epsilon_p \eta)^2} + \frac{1}{\epsilon_p \eta} \right) \right). \quad (3.58)
$$

**Proof.** By inserting $\alpha_d(x)$ into (3.57), we have

$$
\tilde{\mathcal{F}}_d = \frac{12}{\epsilon_p \eta} \int_{-\delta}^{\delta} f_d(x) x^2 dx, \quad f_d(x) \stackrel{\text{def}}{=} \frac{3 + \alpha_d(x)}{(6 + 2\alpha_d(x))^2}.
$$

Since $\lim_{d \to 0} f_d(x) = \frac{3+\frac{\gamma(x)^3}{\epsilon_p \eta}}{(6+2\frac{\gamma(x)^3}{\epsilon_p \eta})^2}$, the dominated convergence theorem gives

$$
\lim_{d \to 0} \tilde{\mathcal{F}}_d = \frac{12}{\epsilon_p \eta} \int_{-\delta}^{\delta} \frac{3 + \frac{\gamma(x)^3}{\epsilon_p \eta}}{(6 + 2 \frac{\gamma(x)^3}{\epsilon_p \eta})^2} \cdot x^2 dx \le \frac{1}{36 \epsilon_p \eta} \left( 3 + \frac{1}{\epsilon_p \eta} \right) \cdot \frac{\delta^3}{3},
$$

which yields (3.58) and completes the proof. □

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

We conclude this section with a series of remarks.

**Remark 3.5 (Darcy versus Navier or Dirichlet).** A salient difference of the of the porous surface model considered in this paper with respect to Navier or Dirichlet boundary conditions lies in the fact that the fluid can escape vertically, i.e. **u**·**n** can be non-zero on Σ, so that ∂₁φ can play an important role in the estimation of the drag force. On the contrary, the work reported in 9 only focuses on the horizontal velocity of the fluid. Note that we neglected the impact of κ in the relaxed energy
30 Champion, Fernández, Grandmont, Vergnet, Vidrascu

(3.48) in order to facilitate the solution of (3.51) (Proposition 3.3). However, we will retrieve the effect of $\kappa$ in the bound of the full drag force $\mathcal{F}_d$, provided in the next section. It should be noted we do not have a natural lower bound on the drag force because the cross derivatives $\partial_2^2\varphi\partial_1^2\varphi$ cannot be simplified in the energy functional (3.43), so that we do not have $\tilde{\mathcal{F}}_d \le \mathcal{F}_d$ as in the Navier or Dirichlet case (see 9).

**Remark 3.6 (Asymptotic of $\tilde{\mathcal{F}}_d$ for $\varepsilon_p\eta = \mathcal{O}(d^p)$, $p \ge 3$).** If $\varepsilon_p\eta = \mathcal{O}(d^p)$ with $p \ge 3$, we obtain an asymptotic behavior for the relaxed drag force (3.4) similar to the case with Dirichlet boundary conditions on both the disk and the wall:

$$ \tilde{\mathcal{F}}_d \sim d^{-\frac{3}{2}}. \qquad (3.59) $$

Indeed, by using $\varepsilon_p\eta = \mathcal{O}(d^p)$, $p \ge 3$, in (3.54) we infer that $\alpha_d \ge 1$. Since $f: x \mapsto \frac{x^2+x}{(6+2x)^2}$ is continuous in $[1, +\infty)$ and $\lim_{x \to \infty} f(x) = \frac{1}{4}$, we have

$$ c \le f(\alpha) = \frac{\alpha^2 + 3\alpha}{(6 + 2\alpha)^2} \le C, $$

from which we deduce

$$ c \int_{-\delta}^{\delta} \frac{x^2}{(d + \gamma(x))^3} dx \le \tilde{\mathcal{F}}_d = \tilde{\mathcal{E}}_d(\tilde{\varphi}) \le 12C \int_{-\delta}^{\delta} \frac{x^2}{(d + \gamma(x))^3} dx. $$

We therefore recover Dirichlet asymptotic behavior as $\int_{-\delta}^{\delta} \frac{x^2}{(d+\gamma(x))^3} dx \sim d^{-\frac{3}{2}}$. We recall that the idea to retrieve this behavior consists in combining the expansion of $\gamma(x) = \frac{x^2}{2} + O_{x \to 0}(x^4)$ with the change of variables $u = \frac{v}{\sqrt{d}}$ to exhibit $d^{-\frac{3}{2}}$ out of the integral (see 9 for the no-slip case).

Note that the obtained asymptotic (3.59) of the relaxed drag force is consistent with the asymptotic analysis carried out at the Section 3.2. Indeed, as stated in Theorem 3.2, when $\eta$ tends to zero the solution of problem (2.8) converges weakly towards solution of a problem with Navier boundary conditions on the wall and Dirichlet boundary conditions on the disk. So, for $\varepsilon_p\eta = \mathcal{O}(d^p)$, $p \ge 3$ we expect to recover (3.59) as the asymptotic of the drag force associated to the Navier-Dirichlet boundary conditions, which is known to behave like the Dirichlet-Dirichlet case (see, e.g., 10). Yet, the proved asymptotic behavior (3.59) only concerns the relaxed drag force $\tilde{\mathcal{F}}_d$. Nothing can be said of the asymptotic of the whole drag force due to the lack of lower bound on $\mathcal{F}_d$ (see Remark 3.5).

**Remark 3.7 (Asymptotic of $\tilde{\mathcal{F}}_d$ for $\varepsilon_p\eta = \mathcal{O}(d^p)$, $0 < p < 3$).** In this case, it can be shown that

$$ \tilde{\mathcal{F}}_d \sim d^{-\frac{p}{2}}. $$

by combining the Taylor expansion $\gamma(x) = \frac{x^2}{2} + O_{x \to 0}(x^4)$ with the change of variable $t = \frac{x}{d^{\frac{1}{6}}}$.
### 3.3.4. Proof of Theorem 3.3

The main idea to complete the proof consists in extending $\tilde{\varphi}$ to the whole domain $\Omega_d$ in order to obtain an admissible velocity $\tilde{u}_d \in \mathcal{A}_d$. To obtain a bound on $\mathcal{F}_d$, we have to estimate the whole energy $\mathcal{E}_d(\tilde{u}_d)$ and, in particular, all the terms left out above in the introduction of the relaxed energy functional $\tilde{\mathcal{E}}_d$. We consider the same extension as in 14 and use similar notations for the auxiliary functions.

**Proof.** We consider the following extension of $\tilde{\varphi}_d$ (see 14):

$$ \tilde{\varphi}_d(x, y) \stackrel{\text{def}}{=} \begin{cases} \chi_{2\delta}(1 + |x|)\tilde{\varphi}_d(x, y) + [1 - \chi_{2\delta}(1 + |x|)]\hat{\varphi}_d(x, y) & (x, y) \in \Omega_{d,2\delta}, \\ \hat{\varphi}_d(x, y) & \text{otherwise,} \end{cases} \qquad (3.60) $$

where the auxiliary functions $\chi_{2\delta}$ and $\hat{\varphi}_d$ are given by

$$ \chi_{\delta}(x) \stackrel{\text{def}}{=} \begin{cases} 1 & \text{if } x \le 1 + \delta/2 \\ 0 & \text{if } x \ge 1 + \delta \end{cases} \qquad \hat{\varphi}_d(x, y) \stackrel{\text{def}}{=} x\chi_{\delta}(\|\mathbf{x} - \mathbf{x}_d\|), \qquad (3.61) $$

with $\mathbf{x}_d \stackrel{\text{def}}{=} (0, 1 + d)$ denoting the center of the disk. In between $1 + \delta/2$ and $1 + \delta$, $\chi_{\delta}$ has a smooth $\mathcal{C}^2$ decreasing transition. The different zones have been depicted in Figure 4 for illustration purposes.

<figure><img src="image_7.png" /><figcaption>Figure 4: Definition of the extension $\tilde{\varphi}_d$.</figcaption></figure>
We denote $\tilde{u}_d = -\partial_2\tilde{\varphi}_d e_1 + \partial_2\tilde{\varphi}_d e_2$. We first check $\tilde{u}_d$ is admissible. We skip the tedious computations, letting them to the reader. According to Proposition 3.2, we have the following bound on $\mathcal{F}_d$

$$ \mathcal{F}_d \le \mathcal{E}_d(\tilde{u}_d). \qquad (3.62) $$

Thus in order to derive (3.27), we estimate the following energy,

$$ \begin{aligned} \mathcal{E}_d(\tilde{u}_d(\tilde{\varphi}_d)) &= \int_{\Omega_d} 2\mathbf{D}(\tilde{u}_d(\tilde{\varphi}_d)) : \mathbf{D}(\tilde{u}_d(\tilde{\varphi}_d)) + \frac{1}{\varepsilon_p \eta} \int_{\Sigma} |\tilde{\varphi}_d|^2 + \frac{\varepsilon_p}{4\kappa} \int_{\Sigma} |\partial_1 \tilde{\varphi}_d|^2 \\ &= \tilde{\mathcal{E}}_d(\tilde{\varphi}_d) + \int_{\Omega_{d,\delta}} \left[ 2\mathbf{D}(\tilde{u}_d(\tilde{\varphi}_d)) : \mathbf{D}(\tilde{u}_d(\tilde{\varphi}_d)) - |\partial_2^2 \tilde{\varphi}_d|^2 \right] \\ &\quad + \frac{\varepsilon_p}{4\kappa} \int_{\Sigma \cap \Omega_{d,\delta}} |\partial_1 \tilde{\varphi}_d|^2 + \int_{\Omega_{d,2\delta} \setminus \Omega_{d,\delta}} 2\mathbf{D}(\tilde{u}_d(\tilde{\varphi}_d)) : \mathbf{D}(\tilde{u}_d(\tilde{\varphi}_d)) \\ &\quad + \frac{1}{\varepsilon_p \eta} \int_{\Sigma \cap \Omega_{d,2\delta} \setminus \Omega_{d,\delta}} |\tilde{\varphi}_d|^2 + \frac{\varepsilon_p}{4\kappa} \int_{\Sigma \cap \Omega_{d,2\delta} \setminus \Omega_{d,\delta}} |\partial_1 \tilde{\varphi}_d|^2 \\ &\quad + \int_{\Omega_d \setminus \Omega_{d,2\delta}} 2\mathbf{D}(\hat{u}_d(\hat{\varphi}_d)) : \mathbf{D}(\hat{u}_d(\hat{\varphi}_d)). \end{aligned} $$

We notice that the first term corresponds to the relaxed drag force, already estimated in Theorem (3.4). We proceed by bounding all the remaining terms.

Under the disk. We consider the following energy

$$ \begin{aligned} E_d(\tilde{\varphi}_d) &\stackrel{\text{def}}{=} \int_{\Omega_{d,\delta}} (2\mathbf{D}(\tilde{\varphi}_d) : \mathbf{D}(\tilde{\varphi}_d) - |\partial_2^2 \tilde{\varphi}_d^2|) + \frac{\varepsilon_p}{4\kappa} \int_{-\delta}^{\delta} |\partial_1 \tilde{\varphi}_d|^2 \\ &= \int_{\Omega_{d,\delta}} \left( |\partial_1^2 \tilde{\varphi}_d|^2 + 4|\partial_1 \partial_2 \tilde{\varphi}_d|^2 - 2\partial_1^2 \tilde{\varphi}_d \partial_2^2 \tilde{\varphi}_d \right) + \frac{\varepsilon_p}{4\kappa} \int_{-\delta}^{\delta} |\partial_1 \tilde{\varphi}_d|^2. \end{aligned} $$

We then estimate each term separately. For the first term, we have

$$ \partial_1^2 \tilde{\varphi}_d(x, y) = -\frac{1}{\varepsilon_p \eta} f''(x)y^3 + \frac{3}{\varepsilon_p \eta} g''(x)y + f''(x), $$

with $\alpha_d$ defined in (3.54) and the auxiliary functions

$$ f(x) \stackrel{\text{def}}{=} \frac{x}{6 + 2\alpha_d(x)}, \quad g(x) \stackrel{\text{def}}{=} \frac{x(d + \gamma(x))^2}{6 + 2\alpha_d(x)}. $$

Since $d < \tilde{L}$, the derivatives of $\alpha_d$ can be bounded for all $x \in [-\delta, \delta]$ as

$$ |\alpha_d(x)| \le \frac{(\tilde{L} + 1)^3}{\varepsilon_p \eta}, \quad |\alpha'_d(x)| \le \frac{3(\tilde{L} + 1)^2}{\varepsilon_p \eta} \gamma'(\delta), \quad |\alpha''_d(x)| \le \frac{C(\tilde{L})}{\varepsilon_p \eta} (\gamma''(\delta) + \gamma'(\delta)). $$

We can show that if $\delta < 1$ there exist strictly positive constants $C_1$ and $C_2$ depending only on $\tilde{L}$ and $\delta$ such that

$$ |f''(x)| \le C_1(\tilde{L}, \delta) \left( \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} \right), \quad |g''(x)| \le C_2(\tilde{L}, \delta) \left( 1 + \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} \right). $$
We finally obtain the estimate of the first term

$$ \int_{\Omega_{d,\delta}} |\partial_1^2 \tilde{\varphi}(x,y)|^2 \le \frac{C_3(\tilde{L}, \delta)}{(\varepsilon_p \eta)^2} \left( 1 + \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} + \frac{1}{(\varepsilon_p \eta)^3} + \frac{1}{(\varepsilon_p \eta)^4} \right). $$

Similarly, we can obtain the following estimates for the other terms

$$ \int_{\Omega_{d,\delta}} |\partial_1 \partial_2 \tilde{\varphi}(x,y)|^2 \le \frac{\delta C_4(\tilde{L}, \delta)}{(\varepsilon_p \eta)^2} \left( 1 + \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} \right), $$

$$ \int_{\Omega_{d,\delta}} |\partial_1 \partial_2 \tilde{\varphi}(x,y)|^2 \le \frac{\delta C_5(\tilde{L}, \delta)}{(\varepsilon_p \eta)^2} \left( 1 + \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} \right), $$

$$ \frac{\varepsilon_p}{4\kappa} \int_{\Sigma \cap \Omega_{d,\delta}} |\partial_1 \tilde{\varphi}(x,0)|^2 \le \frac{\varepsilon_p}{4\kappa} \frac{\delta C_6(\tilde{L}, \delta)}{(\varepsilon_p \eta)^2} \left( 1 + \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} \right). $$

Finally, we obtain the following estimate

$$ |E_d(\tilde{\varphi})| \le \frac{C_7(\tilde{L}, \delta)}{(\varepsilon_p \eta)^2} \left( 1 + \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} + \frac{1}{(\varepsilon_p \eta)^3} + \frac{1}{(\varepsilon_p \eta)^4} \right) \\ + \frac{\varepsilon_p}{4\kappa} \frac{C_8(\tilde{L}, \delta)}{(\varepsilon_p \eta)^2} \left( 1 + \frac{1}{\varepsilon_p \eta} + \frac{1}{(\varepsilon_p \eta)^2} \right). \quad (3.63) $$

It is worth noting that if $\kappa$ or $\eta$ vanishes we lose the upper bound.

Away from the disk $\Omega_d \setminus \Omega_{d,2\delta}$. Outside of $\Omega_{d,2\delta}$, $\tilde{\varphi}_d = \hat{\varphi}_d$ and $\hat{\varphi}_d$ depends on the gap distance $d$ only through a translation (see proof of Lemma 8 of 14). Thus, energy in $\Omega \setminus \Omega_{d,2\delta}$ of $\tilde{u}_d = \nabla^\perp \tilde{\varphi}_d$ does not depend on $d$ and is therefore bounded. So, there exists a constant $C_9(\delta) > 0$ such that

$$ \int_{\Omega_d \setminus \Omega_{d,2\delta}} 2\mathbf{D}(\tilde{u}_d) : \mathbf{D}(\tilde{u}_d) \le C_9(\delta). \quad (3.64) $$

for $d < \tilde{L}$.

Transitional zone $\Omega_{d,2\delta} \setminus \Omega_{d,\delta}$. Since $|\chi_\delta| \le 1$, using Young's inequality, we have
$$ \tilde{\varphi}_d = \chi_\delta \tilde{\varphi}_d + (1-\chi_\delta) \hat{\varphi}_d \Rightarrow |\tilde{\varphi}_d| \le |\tilde{\varphi}_d| + |\hat{\varphi}_d| \Rightarrow |\tilde{\varphi}_d|^2 \le C_{10} (|\tilde{\varphi}_d|^2 + |\hat{\varphi}_d|^2). $$

So, the estimate of the transitional term can simply be obtained by combining the bounds (3.63) and (3.64).

Finally, (3.27) then follows by combining the previous bounds with (3.58), which concludes the proof of Theorem 3.3. □

## 4. Numerical results

In this section we provide numerical evidence of the theoretical results obtained in the previous sections by performing a series of numerical experiments. To this purpose, we consider finite element approximations of the solution to the static system (2.8), for different values of $d$ (viz., the distance between the disk and the
34 Champion, Fernández, Grandmont, Vergnet, Vidrascu

bottom wall), and of $\kappa$ and $\eta$ (the normal and tangential conductivity parameters, respectively). We recover numerically the theoretical upper bound of the drag force of Theorem 3.3, and the weak convergence results stated in Theorem 3.2 when either $\kappa$ or $\eta$ vanishes. Then, we numerically investigate the ability of the porous layer to enable collision as stated Corollary 3.1. Finally, we provide numerical insight on the validity of the fluid-reduced Darcy model (2.8) by comparing it against the full model.

## 4.1. Experiment setup

In what follows, we consider a disk of radius $r = 1$ centered at position $(40, d+r)$, immersed in a rectangular box of width $L = 80$ and of height $\tilde{L} = 40$, whose bottom-left corner lies at position $(0, 0)$. We numerically solve the system (2.8) for different values of the gap $d$ between the disk and the bottom wall, decreasing from 1 to $10^{-3}$. We do the same with system (2.8a)-(2.8d) supplemented with Navier boundary conditions (3.20). The values for $d$ are chosen to be equidistant in log scale

$$d \in \{0.001, 0.002, 0.004, 0.008, 0.016, 0.032, 0.064, 0.126, 0.252, 0.502, 1\}. \quad (4.1)$$

This experiment setup is very similar to the one considered in 15. Therein the authors retrieve the $d^{-3/2}$ asymptotic behavior of the drag force for the problem with Dirichlet boundary conditions (on both the disk and the bottom wall).

For a more intuitive visualization and without loss of generality, we choose the disk to move down with unit velocity, instead of going up as originally in system (2.8). One can straightforwardly check that this only changes the sign of the computed drag force, so it does not affect the expected theoretical results.

For a given gap distance $d$, let $\mathcal{T}_d^f$ denote a triangulation of the fluid domain $\Omega_d$ and let $\mathcal{T}_d^p$ the 1D triangulation of the porous mid-surface $\Sigma$, given simply as the edges of $\mathcal{T}_d^f$ lying on $\Sigma$. For the Navier problem, we consider the standard space of continuous piecewise polynomial functions of degree two for the fluid velocity, namely,

$$\mathcal{U}_{d,h}^{Na} \stackrel{\text{def}}{=} \left\{ \mathbf{v}_h \in C_0(\bar{\Omega}_d)^2 : \mathbf{v}_h|_K \in \mathbb{P}_2(K)^2 \forall K \in \mathcal{T}_d^f, \mathbf{v}_h = \mathbf{e}_2 \text{ on } \partial S_d, \mathbf{v}_h = \mathbf{0} \text{ on } \Gamma, \mathbf{v}_h \cdot \mathbf{n} = 0 \text{ on } \Sigma \right\},$$

and of continuous piecewise linear functions for the fluid pressure

$$Q_{d,h} \stackrel{\text{def}}{=} \{ q_h \in C_0(\bar{\Omega}_d) : q_h|_K \in \mathbb{P}_1(K) \ \forall K \in \mathcal{T}_d^f \}.$$
We also introduce the homogeneous counterpart of $\mathcal{U}_{d,h}^{Na}$

$$ \mathcal{U}_{d,h}^{0,Na} \stackrel{\text{def}}{=} \left\{ \mathbf{v}_h \in C_0(\overline{\Omega_d})^2 : \mathbf{v}_h|_K \in \mathbb{P}_2(K)^2 \ \forall K \in \mathcal{T}_d^f, \mathbf{v}_h = \mathbf{0} \text{ on } \Gamma \cup \partial S_d, \mathbf{v}_h \cdot \mathbf{n} = 0 \text{ on } \Sigma \right\}. $$

For any $d > 0$, the corresponding finite element approximation of (2.8a)-(2.8d) with the Navier boundary conditions writes

$$ \begin{array}{l} \text{Find } (\mathbf{u}_h, p_h) \in \mathcal{U}_{d,h}^{Na} \times Q_h \text{ such that} \\ \quad 2(\mathbf{D}(\mathbf{u}_h), \mathbf{D}(\mathbf{v}_h))_{\Omega_d} + (\mathbf{div} \, \mathbf{u}_h, q_h)_{\Omega_d} - (p_h, \mathbf{div} \, \mathbf{v}_h)_{\Omega_d} + \varepsilon_s (p_h, q_h)_{\Omega_d} = 0 \quad (4.2) \\ \forall (\mathbf{v}_h, q_h) \in \mathcal{U}_{d,h}^{0,Na} \times Q_{d,h}. \end{array} $$

Here, the small perturbation $\varepsilon_s (p_h, q_h)_{\Omega_d}$, with $\varepsilon_s = 10^{-10}$, serves to guarantee that the pressure has zero mean.

To formulate the finite element approximation of (3.2), we consider a slightly different trial space for the fluid velocity, without any constraint on $\Sigma$,

$$ \mathcal{U}_{d,h} \stackrel{\text{def}}{=} \left\{ \mathbf{v}_h \in C_0(\overline{\Omega_d})^2 : \mathbf{v}_h|_K \in \mathbb{P}_2(K)^2 \ \forall K \in \mathcal{T}_d^f, \mathbf{v}_h = \mathbf{e}_2 \text{ on } \partial S_d, \mathbf{v}_h = 0 \text{ on } \Gamma \right\}, $$

with its homogeneous counterpart

$$ \mathcal{U}_{d,h}^0 \stackrel{\text{def}}{=} \left\{ \mathbf{v}_h \in C_0(\overline{\Omega_d})^2 : \mathbf{v}_h|_K \in \mathbb{P}_2(K)^2 \ \forall K \in \mathcal{T}_d^f, \mathbf{v}_h = 0 \text{ on } \Gamma \cup \partial S_d \right\}. $$

The finite element approximation of (3.2) finally writes

$$ \begin{array}{l} \text{Find } (\mathbf{u}_h, \hat{p}_h, p_h) \in \mathcal{U}_{d,h} \times \mathcal{D}_{d,h} \times Q_{d,h} \text{ such that} \\ \quad 2(\mathbf{D}(\mathbf{u}_{d,h}), \mathbf{D}(\mathbf{v}_{d,h}))_{\Omega_d} + \frac{\varepsilon_p}{4\kappa} (\mathbf{u}_{d,h} \cdot \mathbf{n}, \mathbf{v}_h \cdot \mathbf{n})_{\Sigma} \\ \quad + (\mathbf{div} \, \mathbf{u}_{d,h}, q_h)_{\Omega_d} - (p_{d,h}, \mathbf{div} \, \mathbf{v}_h)_{\Omega_d} + \varepsilon_p \eta (\nabla_{\tau} \hat{p}_{d,h}, \nabla_{\tau} \hat{q}_h)_{\Sigma} \quad (4.3) \\ \quad + (\hat{p}_h, \mathbf{v}_h \cdot \mathbf{n})_{\Sigma} - (\mathbf{u}_h \cdot \mathbf{n}, \hat{q}_h)_{\Sigma} + \varepsilon_s^{-1} \int_{\Sigma} \hat{p}_h \int_{\Sigma} \hat{q}_h = 0, \\ \forall (\mathbf{v}_h, \hat{q}_h, q_h) \in \mathcal{U}_{d,h}^0 \times \mathcal{D}_{d,h} \times Q_{d,h}. \end{array} $$

The last term enforces the mean value of the Darcy pressure $\hat{p}_h$ to be 0, with again $\varepsilon_s = 10^{-10}$.

For each value of the gap (4.1), a mesh of $\Omega_d$ has been generated. The mesh size between the disk and the bottom wall depends on the gap $d$: we have at least ten triangles between the disk and the bottom wall in order to properly catch the fluid dynamics in the contact zone. The mesh size below the disk is defined by

$$ h_{min} = \min \left\{ \frac{d}{10}, 0.01 \right\}. $$
The meshes are not structured uniform?, to avoid unnecessary refinement outside
of the contact zone. Away from the disk, the mesh size is $h_{max}$ = 0.1. In order to
ensure a smooth transition between the small mesh elements below the disk and
the bigger mesh elements away from the gap, we first build the meshes using the
FreeFem++ (see 13) integrated mesh generator then we optimize them with the
Mmg remeshing software (see 5). The resulting meshes have between $10^4$ and $10^5$
nodes.

<figure><img src="image_8.png" /><figcaption>Figure 5: Mesh for the smaller gap distance $d = 10^{-3}$. Zooms on the disk and on the contact zone.</figcaption></figure>

**4.2. Results**

First, we numerically validate the porous layer model. We then investigate that we can numerically recover the main theoretical results of the paper, that is, the asymptotic analysis in $\kappa$ and $\eta$ of Theorem 3.2 and the estimate of the drag force stated in Theorem 3.3. Finally, we use the numerical simulations to bring insight on the way the porous layer model influences the fluid and solid dynamics, by reporting the fluid velocity and pressure fields and also by simulating the disk fall towards the wall.

**Validation of the porous layer model.** The porous layer model is derived in 4 by averaging the original Darcy system across the thickness of the bottom wall under suitable assumptions. In order to investigate the validity of the surface Darcy model, we simulate the fall of the disk on a 2D porous layer of the same thickness $\varepsilon_p = 0.01$. Figures 6 and 7 report the velocity and pressure field obtained in both cases and indicate that the behavior of 2D Darcy is very similar to the Darcy surface model considered in this paper. In Figures 6b and 7b, the porous layer thickness has been amplified by a factor 25 to better visualize the Darcy pressure.

**Drag force asymptotics.** In Figure 8, we report the drag force $F_d$ for $d \in [10^{-3}, 1]$ arising from system (2.8) where the fluid is coupled to the porous layer model.
<figure><img src="image_9.png" /><figcaption>Figure 6: Comparison of velocity field with the porous layer model versus Darcy 2D.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Comparison of pressure field with the porous layer model versus Darcy 2D.</figcaption></figure>

For the sake of simplicity, we call it the Darcy drag force. We also plot the drag force asymptotic arising from system (2.8a)-(2.8d) with Navier boundary conditions (3.20), that we similarly call Navier drag force asymptotic.

We test various set of values for the conductivity parameters $κ$ and $η$, notably to investigate the behavior when the conductivity parameters gets smaller. Thus, we can compare with the asymptotic results of Theorem 3.2 and the estimate of the Darcy drag force (3.27) stated in Theorem 3.3.

Away from the disk at $d = 1$, the wall does not have any influence on the disk, so there is practically no difference between Navier and Darcy drag force. When the disk gets closer of the wall, the Darcy drag force asymptotic reaches a plateau which is consistent with the result of Theorem 3.3, whereas for Navier the drag force explodes. The very good recovery of the theoretical slope of -0.5 for Navier drag force is clearly visible.

As expected from the estimate on the drag force (3.27), the plateau value of the Darcy asymptotic increases when $κ$ and $η$ goes to 0. We can also note that the conductivity parameters $κ$ and $η$ play a role in the switch from the Navier
<figure><img src="image_11.png" /><figcaption>Figure 8: Drag force asymptotic when $d \to 0$ and $(\kappa, \eta) \to 0$.</figcaption></figure>

asymptotic to a plateau, smaller conductivity values leading to a switch to a plateau closer to the wall. The combination of both effects make that the Darcy drag force asymptotic converges towards the Navier one, when $\eta$ and/or $\kappa$ goes to 0, which is consistent with the asymptotic analysis results of Theorem 3.2. It is worth noting
that the convergence is faster when the tangential conductivity $\eta$ vanishes than
when the normal one $\kappa$ vanishes.

It is interesting to see that the porous layer simply removes the singularity of
the drag force arising in the case of Navier boundary conditions. The introduction
of a porous modeling of the asperities solves the no-collision paradox by enabling
contact, but it also changes the contact dynamics with respect to Navier.

**Navier versus Darcy velocity and pressure fields.** In order to better understand the fundamental difference between the Navier and Darcy drag force asymptotics, Figures 9 and 10 show the velocity and pressure fields, respectively obtained with the Navier and Darcy models, when the disk is very close to the wall ($d = 0.032$).

<figure><img src="image_12.png" /><figcaption>Figure 9: Velocity and pressure fields for a gap distance $d = 0.032$ obtained with Navier boundary conditions on the below wall.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: Velocity and pressure fields for a gap distance $d = 0.032$ and a porous layer model with conductivity parameters $\kappa = 1$ and $\eta = 1$ on the below wall.</figcaption></figure>

In the case of the coupling with a porous layer, Figure 10, we can clearly see that
the fluid velocity on the wall right below the disk is vertical and of magnitude 1,
which yields a plateau of fluid entering into the porous layer below the disk. In this
new regime where the disk is in the close vicinity of the wall, all the fluid escapes
vertically through the porous layer without any tangential leak. Thus, reducing the
gap does not affect the flow of the fluid. On the contrary, with Navier boundary
conditions on the wall, Figure 9 indicates that size of the escaping jet can only
increase whenever $d$ gets smaller, as only tangential escape is allowed.
**Porous layer dynamics when $\eta \to 0$ or $\kappa \to 0$.** We now focus on the influence of $\kappa$ and $\eta$ on the porous layer dynamics. For that purpose, we fix the disk at a tiny distance $d = 4.10^{-3}$ from the wall and we look at the evolution of the incoming fluid flow $u_h \cdot n$ on $\Sigma$ and of the Darcy pressure $\hat{p}_h$ when $\kappa$ or $\eta$ vanishes, see Figures 11 and 12. We zoom around the contact point at $x = 40$.

<figure><img src="image_14.png" /><figcaption>Figure 11: Zoom below the disk of $u_d \cdot n$ and $\hat{p}_d$ when $\kappa \to 0$.</figcaption></figure>

For $\kappa = \eta = 1$, a plateau that stagnates at $-1$ is observed, which corresponds to the velocity of the disk. When the disk is in the close vicinity of the wall, all the fluid leaks into the porous medium without any tangential escape, as we have already seen in Figure 10. Reducing the conductivity hinders the fluid escape: the plateau first narrows then it shrinks to 0. We recover in both cases $\kappa \to 0$ and $\eta \to 0$ the convergence of $u_h \cdot n$ to 0, which is consistent with weak convergence results of Darcy solutions towards Navier when either $\kappa$ or $\eta$ tend to 0, see Theorem 3.2.

It is interesting to highlight the differences between the two cases. Let us look at the behavior of the Darcy pressure. When $\eta \to 0$, it becomes more and more singular below the disk, converging towards the Navier normal stresses accordingly to the proved weak convergence of Theorem 3.2, with a huge peak of the Darcy pressure up to $8.10^3$ at $\eta = 10^{-6}$, see Figure 12. This singularity makes it difficult to run the simulation for smaller values of $\eta$ because of numerical instabilities. It
<figure><img src="image_15.png" /><figcaption>Figure 12: Zoom below the disk of $u_d \cdot n$ and $\hat{p}_d - \frac{1}{|\Omega_d|} \int_{\Omega_d} p_d$ when $\eta \to 0$.</figcaption></figure>

also explains why the Darcy drag force converges much faster towards the Navier one for $\eta \to 0$ than $\kappa \to 0$, as we can see on Figure 8. On the contrary, reducing $\kappa$ makes $\hat{p}_h$ flatter, consistently with the $H^1$ strong convergence of $\hat{p}_h$ towards 0 when $\kappa \to 0$ proved in Theorem 3.2.

The difference between the two asymptotics $\kappa \to 0$ and $\eta \to 0$ is also noticeable if we look at the way the fluid escapes in Figures 11 and 12: the two flat bosses on both side of the disk when the normal conductivity $\eta$ goes to 0 contrast with the two narrow peaks of outling fluid appearing when the tangential conductivity $\eta$ goes to 0.

**Disk trajectories with Darcy layer on the wall versus with Navier boundary conditions.** The numerical approximations obtained for the drag force, can be combined with the ODE (2.2) to simulate the evolution of the gap distance $d(t)$. We recall that the disk dynamics are given by the relation

$$ \ddot{d}(t) + \dot{d}(t) \mathcal{F}_{d(t)} = 0, \quad (4.4) $$

where the drag force depends on the roughness model considered for the bottom wall (here, Navier boundary conditions or the porous layer model).

In order to approximate $d$, we first consider a linear interpolation of $\log(\mathcal{F}_d)$ that we then insert in an explicit time-stepping of (4.4). In the following, we fix
<figure><img src="image_16.png" /><figcaption>Figure 13: Evolution of the gap distance for a ball starting at distance $d_0 = 1$ with reduced Darcy on the wall and for various initial conditions.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 14: Comparison of the evolution of the gap distance between Navier and Darcy porous layer for a ball starting at distance $d_0 = 1$.</figcaption></figure>

the conductivity parameters $\kappa$ and $\eta$ to be equal to 1. In Figure 8, we see that the bound of the Darcy drag force is around $C = 100$. We first consider $d_0 = 1$ and then start a downward motion with $|\dot{d}_0|$ going from 1 to 100, see Figure 13. As expected
from Corollary 3.1, we have collision for $|\dot{d}_0| \le Cd_0$.

We now compare the effect of Darcy and Navier on the ball dynamics in Figure 14, which provides the trajectories obtained in both cases for different initial velocities. For $\dot{d}_0 = -30$, we see that the disk reaches the wall in the case of the Darcy layer, whereas with Navier boundary conditions the wall stops the disk at a distance $d = 0.1$.

## 5. Conclusion

In this paper, we have analyzed the contact capabilities of a fluid-structure interaction model with seepage reported in 4. The key feature of this model lies in taking into account the surface roughness of the contacting wall in terms of a reduced Darcy model. The analysis shows that this modeling approach removes the no-collision paradox. The contact dynamics are also modified with respect to more standard boundary conditions, such as Navier, since contact is allowed with a non-zero velocity. A non penetration condition must hence be added in order to prevent the solid to go through the contacting wall.

Extensions of this work can explore several directions. From the mathematical analysis point of view, one interesting question would be the study of the complete fluid-structure-contact interaction model, with appropriate non penetration conditions. Another interesting question would be the formulation of the model in the case of multiple solids getting into contact.

## Acknowledgment

MF and MV were partially supported by the French National Research Agency (ANR), through the SIMR project (ANR-19-CE45-0020).

## Bibliography

1. C. Ager, B. Schott, A.-T. Vuong, A. Popp and W. A. Wall, A consistent approach for fluid-structure-contact interaction based on a porous flow model for rough surface contact, *International Journal for Numerical Methods in Engineering* **119** (2019) 1345–1378.

2. M. Bogovskii, Solution of the first boundary value problem for an equation of continuity of an incompressible medium, *Soviet Math. Dokl.* **20** (1979) 1094–1098.

3. F. Boyer and P. Fabrie, *Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models*, volume 183 (Springer Science & Business Media, 2012).

4. E. Burman, M. A. Fernández, S. Frei and F. M. Gerosa, A mechanically consistent model for fluid-structure interactions with contact including seepage, *Computer Methods in Applied Mechanics and Engineering* **392** (2022) 114637.

5. C. Dobrzynski and P. Frey, Anisotropic delaunay mesh adaptation for unsteady simulations, in *Proceedings of the 17th international Meshing Roundtable* (Springer, 2008), pp. 177–194.

6. Q. Du, M. D. Gunzburger, L. S. Hou and J. Lee, Analysis of a linear fluid-structure interaction problem, *Discrete and continuous dynamical systems* **9** (2003) 633–650.
44 Champion, Fernández, Grandmont, Vergnet, Vidrascu

7. A. Ern and J.-L. Guermond, *Theory and practice of finite elements*, volume 159 (Springer, 2004).

8. D. Gérard-Varet and M. Hillairet, Regularity issues in the problem of fluid structure interaction, *Archive for rational mechanics and analysis* **195** (2010) 375–407.

9. D. Gérard-Varet and M. Hillairet, Computation of the drag force on a sphere close to a wall: the roughness issue, *ESAIM: Mathematical Modelling and Numerical Analysis* **46** (2012) 1201–1224.

10. D. Gérard-Varet, M. Hillairet and C. Wang, The influence of boundary conditions on the contact problem in a 3d navier–stokes flow, *Journal de Mathématiques Pures et Appliquées* **103** (2015) 1–38.

11. V. Girault and P.-A. Raviart, *Finite element approximation of the Navier-Stokes equations*, volume 749 (Springer Berlin, 1979).

12. D. Gérard-Varet and M. Hillairet, Existence of weak solutions up to collision for viscous fluid-solid systems with slip, *Comm. Pure Appl. Math.* **67** (2014) 2022–2075.

13. F. Hecht, New development in freefem++, *J. Numer. Math.* **20** (2012) 251–265.

14. M. Hillairet, Lack of collision between solid bodies in a 2d incompressible viscous flow, *Communications in Partial Differential Equations* **32** (2007) 1345–1371.

15. M. Hillairet, A. Lozinski and M. Szopos, On discretization in time in simulations of particulate flows, 2010.

16. M. Hillairet and T. Takahashi, Collision in 3d fluid structure interactions problems, *SIAM Journal on Mathematical Analysis*.

17. M. Hillairet and T. Takahashi, Blow up and grazing collision in viscous fluid solid interaction systems, *Annales de l'Institut Henri Poincaré C, Analyse non linéaire* **27** (2010) 291–313.

18. M. Hillairet and T. Takahashi, Existence of contacts for the motion of a rigid body into a viscous incompressible fluid with the tresca boundary conditions, *Tunisian Journal of Mathematics* **3** (2021) 447–468.

19. D. Kamensky, F. Xu, C.-H. Lee, J. Yan, Y. Bazilevs and M.-C. Hsu, A contact formulation based on a volumetric potential: Application to isogeometric simulations of atrioventricular valves, *Computer Methods in Applied Mechanics and Engineering* **330** (2018) 522–546.

20. N. Khaledian, P.-F. Villard and M.-O. Berger, Capturing contact in mitral valve dynamic closure with fluid-structure interaction simulation, *International Journal of Computer Assisted Radiology and Surgery* **17** (2022) 1391–1398.

21. V. Martin, J. Jaffré and J. E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, *SIAM Journal on Scientific Computing* **26** (2005) 1667–1691.

22. C. Wang, Strong solutions for the fluid-solid systems in a 2-d domain, *Asymptotic Analysis* **89** (2014) 263–306.