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

1.   Introduction

The numerical simulation of systems involving fluid-structure-contact interaction is of fundamental importance in many engineering and biomedical applications. For instance, contact modeling is a crucial ingredient in the simulation of the dynamics of native or artificial cardiac valves (see, e.g., ^[22,23]).

Modeling contact between solids within a fluid-structure interaction (FSI) framework raises many modeling, mathematical and numerical issues. First, in the case of a ball immersed in a viscous incompressible fluid with no-slip boundary conditions and falling over a fixed horizontal plane, the resulting FSI models are unable to predict contact, both in 2D (see, e.g., ^[17]) and in 3D (see, e.g., ^[19]), 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 seems possible only in very specific configurations like, for example, grazing collision in 3D (see ^[21]). On the other hand, many studies show that 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., ^[20,29] in 2D and ^[12,13,15] in 3D) or by considering rough solids (see, e.g., ^[11] in 2D and ^[12] 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 release 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,7]). 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 ^[7]. 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 a 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 ^[12], 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 ^[12,17], 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 ^[7]. An appropriate velocity scaling of the problem (see, e.g., ^[12]) 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 modeling of the contact wall, 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., ^[12,17,20]), 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) \overset{ \smash{ \mathrm{def} }}{ = } \Omega\setminus S(t)$ and its associated non-cylindrical trajectory

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:

where $({\mathit{\boldsymbol{e}}}_1, {\mathit{\boldsymbol{e}}}_2)$ denotes the canonical basis of $\mathbb{R}^2$, $B({\mathit{\boldsymbol{x}}}, R)$ the disk of radius $R$ centered in ${\mathit{\boldsymbol{x}}}\in \mathbb{R}^2$, and ${d} : \mathbb{R}^+ \rightarrow \mathbb{R}$ stands for the 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 $\mathit{\boldsymbol{u}}:\mathcal T \rightarrow \mathbb{R}^2$ and pressure $p:\mathcal T \rightarrow \mathbb{R}$ fields, which are governed by the following Stokes system

where the fluid Cauchy stress tensor $\boldsymbol{\sigma}$ reads

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

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) and (2.2) have to be coupled with standard geometric, kinematic and dynamic interface conditions on $\partial S(t)$, namely:

where ${\mathit{\boldsymbol{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., ^[17]) in contrast to what is physically observed. In order to circumvent this issue, we consider the alternative modeling approach reported in ^[7], which consists in taking into account the roughness of the contact wall through a thin-walled porous model of fluid seepage.

The matrix of the porous layer is assumed to be rigid and thin-walled, with thickness $\varepsilon_{\rm p} > 0$ around its mid-surface $\Sigma$. We denote by $\boldsymbol{\tau}$ the tangential vector on the surface $\Sigma$; here $\boldsymbol{\tau}$ is simply equal to ${\mathit{\boldsymbol{e}}}_1$. The normal and tangential conductivities of the porous medium are denoted by $\kappa$ and $\eta$, respectively. The averaged porous fluid pressure $\hat{p}: \Sigma\times \mathbb{R}^+\rightarrow \mathbb{R}$ can be described by the following surface Darcy model:

where $u_\Sigma$ denotes the seepage velocity. This model is derived in ^[7,25] by averaging the original Darcy system across the thickness and by applying a trapezoidal rule to approximate the integrals involving the interstitial pressure and the normal velocity across the thickness. Note that (2.4)₁ enforces mass conservation over the porous layer $\Sigma$. In particular, from (2.4), it follows that the net seepage velocity vanishes, viz., ${\int}_\Sigma u_\Sigma = 0$.

Remark 1. As $\Sigma$ is a 1D line, the symbols $\operatorname{div}_{ \boldsymbol \tau}$ and $\boldsymbol \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 equation (2.1) on $\Sigma$, via appropriate kinematic and dynamic relations, namely,

The first relation simply enforces that the seepage velocity in the porous layer has to be the normal component of the bulk fluid velocity. The second corresponds to the balance at the interface of normal stress between the bulk and porous fluids. In particular, the quantity $\hat{p}+\frac{\varepsilon_{\rm p}}{4 \kappa} \mathit{\boldsymbol{u}}\cdot {\mathit{\boldsymbol{n}}}$ represents the reconstruction of the interfacial porous pressure from the averaged porous pressure $\hat{p}$. Finally, free slip tangential stress is enforced with the last relation.

Remark 2. It should be noted that, for the sake of simplicity, a pure slip condition has been considered in (2.5)₃. In practice, it is however more appropriate to consider the so-called Beavers-Joseph-Saffman (BJS) law:

where the coefficient $\alpha\geq 0$ depends on the structure of the porous layer (see, e.g., ^[3,26,27]). As indicated in Remark 11, the analysis below can be adapted straightforwardly to the case in which (2.5)₃ is replaced by the BJS condition (2.6) and, in particular, it does not modify the collision result of Section 3.3.1.

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

with ${d}(0) = {d}_0$ and $\dot {d} (0) = \dot {d}_0$ as initial conditions.

Velocity scaling. System (2.7) 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 ^[12], which leverages the linearity of problem (2.7) with respect to $\dot d$ and basically consists in performing the following change of variables for $t > 0$:

By inserting the above relations into (2.7) 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

the triplet $\mathit{\boldsymbol{u}}_ {d} : \Omega_ {d} \rightarrow \mathbb{R}^2$, $p_ {d} : \Omega_ {d} \rightarrow \mathbb{R}$, $\hat{p}_ {d} : \Sigma \rightarrow \mathbb{R}$ is solution of the following pure (steady) fluid coupled problem:

Note that the velocity scaling (2.8) yields a constant unit velocity as kinematic constraint in (2.9d), 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 (\mathit{\boldsymbol{u}}_{d}, p_ {d}, \hat{p}_{d})$ is well defined, so that the change of variables (2.8) is feasible.

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

with the notation

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

We will show in Section 3.3 that $\mathcal{F}_{{d}}\geq 0$, so that from (2.11) we can infer that the fluid-solid coupling (2.7) 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.10), which only depends on the instantaneous gap ${d}$, trough the solution of the pure fluid problem (2.9). The analysis of the system (2.9) 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., ^{[12,13,17,19,29]}) 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

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

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 a weak solution to the coupled system (2.9). 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

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.9), 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.9) 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 a weak solution to system (2.9) 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^2_0(\omega)$ the space of $L^2(\omega)$ functions with zero mean value in $\omega$ and by $H^1_{0, \Upsilon}(\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 $\Vert \cdot \Vert_{0, \omega}$ and the $H^1(\omega)$ norm by $\Vert \cdot \Vert_{1, \omega}$. We recall that the $H^1$ semi-norm $\vert \cdot \vert_{1, \omega}$ is a norm in $H^1_{0, \Upsilon}(\omega)$, since the Poincaré inequality holds true for functions in $H^1_{0, \Upsilon}(\omega)$. On $\Upsilon$, we also consider the $H^s$ norms, denoted by $\|\cdot\|_{s, \Upsilon}$, for every $0 < s\leqslant 1$, and the $H^1$ semi-norm, denoted by $\vert \cdot \vert_{1, \Upsilon}$.

Since only the gradient of $p_ {d}$ and of $\hat{p}_{d}$ is involved in (2.9a) and (2.9g), respectively, both $p_ {d}$ and $\hat{p}_{d}$ are in principle defined up to a constant. Nonetheless, owing to (2.9e), 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

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,

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

for all $(\mathit{\boldsymbol{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 1. For any ${d}\in(0, {\tilde l}-2)$, problem (3.2) admits a unique solution $(\mathit{\boldsymbol{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 $(\mathit{\boldsymbol{u}}_{d}, p_ {d}, \hat{p}_{d})$ satisfies the following a priori estimates

Proof. The sketch of the proof is as follows. We first introduce a lifting of the non-homogeneous Dirichlet boundary condition (2.9d) and we derive the associated homogeneous weak formulation. Then, we establish an inf-sup condition in order to prove existence and uniqueness of a solution to problem (3.2). The end of the proof deals with the derivation of estimates (3.3)–(3.5). For the sake of clarity, we drop the ${d}$ subscript of all variables in the remaining of the proof, since ${d} > 0$ is fixed. The steps of the proof are detailed hereafter.

Homogeneous problem. We first consider a lifting of the non-homogeneous boundary condition (2.9d). Since ${\int}_{\partial S_ {d}}{\mathit{\boldsymbol{e}}}_2\cdot {\mathit{\boldsymbol{n}}} = 0$, it follows from Theorem 5.1 and Remark 5.3 of ^[14] that there exists a unique $({\mathit{\boldsymbol{w}}}, q) \in H^1(\Omega_ {d})\times L^2_0(\Omega_ {d})$, solution of the Stokes systems (2.9a) and (2.9b) with the following Dirichlet boundary conditions:

Moreover, there exists a constant $C_1 > 0$ which only depends on $\Omega_ {d}$ such that

Therefore, we define for every ${\mathit{\boldsymbol{u}}}\in \mathcal{U}_ {d}$ the new velocity $\bar{\mathit{\boldsymbol{u}}} \in \mathcal{U}_ {d}^0$ given by

Rewriting (3.2) with $\mathit{\boldsymbol{u}} = \bar{\mathit{\boldsymbol{u}}} + {\mathit{\boldsymbol{w}}}$, we obtain that the triplet $\bar{\mathit{\boldsymbol{u}}} : \Omega_ {d} \rightarrow \mathbb{R}^2$, $p : \Omega_ {d} \rightarrow \mathbb{R}$, $\hat{p} : \Sigma \rightarrow \mathbb{R}$ is solution of the following homogeneous problem: Find $(\bar{\mathit{\boldsymbol{u}}}, \hat{p}, p)\in \mathcal{U}_ {d}^0 \times \mathcal{D} \times L^2(\Omega_ {d})$ such that

for all $(\mathit{\boldsymbol{v}}, \hat{q}, q)\in \mathcal{U}_ {d}^0 \times \mathcal{D} \times L^2(\Omega_ {d})$. Problems (3.9) and (3.2) are obviously equivalent. We will then prove existence and uniqueness of a solution $(\bar{\mathit{\boldsymbol{u}}}, \hat{p}, p)$ to problem (3.9) and deduce that problem (3.2) is also well-posed.

For the sake of simplicity, we rewrite the weak formulation (3.9) with usual notations for saddle point problem so the reader can easily find his way through the standard arguments to prove existence and uniqueness. We consider the Hilbert spaces

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

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

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., [10, Theorem 3.77]), it follows that the 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

for all $q\in L^2(\Omega_ {d})\setminus \{0\}$. From Theorem 3 of ^[9] and by adapting Bogovskii's lemma stated in ^[5] 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 ${\mathit{\boldsymbol{w}}}_q\in \mathcal{U}_ {d}^0$ such that $\operatorname{{div}}{\mathit{\boldsymbol{w}}}_q = q$ and $\Vert {\mathit{\boldsymbol{w}}}_q\Vert_{1, \Omega_ {d}} \leqslant C \Vert q\Vert_{0, \Omega_ {d}}$. As a result

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

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., [14, Corollary 4.1]). 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 $\mathit{\boldsymbol{u}}$ and $\hat{p}$, we first consider $(\bar{\mathit{\boldsymbol{u}}}, \hat{p})\in X$ as a test function in the weak formulation (3.9), which yields

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

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

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

Thanks to (3.3), we can derive a sharper estimate on $\hat{p}_{d}$. Indeed using the variational formulation (3.2) with $\mathit{\boldsymbol{v}} = {\mathit{\boldsymbol{0}}}$ and $\hat{q} = \hat{p}_{d}$, we obtain

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

with $C > 0$. From (3.3), we get

from which we finally deduce estimate (3.4) with $C_2 = {2}C\sqrt{C_1}$.

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

for all ${\mathit{\boldsymbol{v}}}\in \mathcal{U}_ {d}^0$. Therefore, owing to (3.6) and (3.8), we have

for all ${\mathit{\boldsymbol{v}}}\in \mathcal{U}_ {d}^0$. By combining this expression with the inf-sup condition (3.12), we get the following pressure estimate

so that, using Cauchy-Schwarz inequality, the continuity of the trace operator and that $\Vert (\mathit{\boldsymbol{v}}, \hat{q}) \Vert_{ \mathcal{U}_ {d}^0 \times \mathcal{D}} \geq \Vert \mathit{\boldsymbol{v}} \Vert_{ \mathcal{U}_ {d}^0 }$, we finally have

with $\tilde C > 0$ a positive constant. Estimate (3.5) then follows from (3.3) and (3.4), which concludes the proof of Theorem 1. □

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 $\mathit{\boldsymbol{u}}_{d}$ is bounded in $H^1(\Omega_ {d})$, independently of $\kappa$ and $\eta$. On the other hand, we loose 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^2_0\oplus\mathbb{R}$, so that

Proposition 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 to system (2.9), satisfies the a priori estimate

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

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

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

for all

According to Bogovskii's Lemma ^[5], there exists $C_4 > 0$ such that for all $p^*\in L^2_0(\Omega_ {d})$ there exists $\tilde{{\mathit{\boldsymbol{v}}}}\in H^1_0(\Omega_ {d})$ such that $\operatorname{div}\tilde{{\mathit{\boldsymbol{v}}}} = p^*$ and

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

This completes the proof. □

Remark 3. 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}$

for all $\mathit{\boldsymbol{v}} \in \mathcal{U}_ {d}^0$, where $H^{\frac{1}{2}}_{00}(\Sigma)$ denotes the Lions-Magenes space and $\big(H^{\frac12}_{00}(\Sigma)\big)^\prime$ its dual (see Remark 4 below). Thanks to the $H^1$ regularity of $\hat{p}_{d}$ and the $H^{\frac{1}{2}}$ regularity of the trace of $\mathit{\boldsymbol{u}}_{d}\cdot {\mathit{\boldsymbol{n}}}$ and using the Hahn-Banach theorem and density arguments, we infer that

We can get rid of the Darcy pressure term by integrating this relation over $\Sigma$ and by using that $\hat{p}_{d}\in L^2_0(\Omega_ {d})$ and ${\int}_\Sigma \mathit{\boldsymbol{u}}_{d}\cdot{\mathit{\boldsymbol{n}}} = 0$, which yields ${\int}_{\Sigma} \boldsymbol{\sigma}(\mathit{\boldsymbol{u}}_{d}, p_ {d}^*){\mathit{\boldsymbol{n}}} \cdot{\mathit{\boldsymbol{n}}} -2Lc_ {d} = - {\int}_\Sigma \hat{p}_{d}-\frac{\varepsilon_{\rm p}}{4 \kappa} {\int}_\Sigma \mathit{\boldsymbol{u}}_{d}\cdot {\mathit{\boldsymbol{n}}} = 0,$ so that the pressure constant is given by

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

Remark 4. The space $H^{\frac{1}{2}}_{00}(\Sigma)$ can be characterized as the trace on $\Sigma$ of functions belonging to $H^1_{0, \Gamma}(\Omega)$ (see ^[24,28]). It is hence a sub-space of $H^{\frac{1}{2}}(\Sigma)$ but with further regularity. In particular, for $\xi \in H^{\frac{1}{2}}_{00}(\Sigma)$ we have that $x\in(- {l}, {l})\mapsto\frac{\xi(x)}{\sqrt{x+ {l}}}$ belongs to $L^2(\Sigma)$.

3.2.   Asymptotic analysis

In this section, we use the a priori estimates of Theorem 1 to study the asymptotic behavior of system (2.9) with respect to $\kappa$ or $\eta$. These results are stated in Theorem 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.9). We look for the fluid velocity $\mathit{\boldsymbol{u}}_ {d}:\Omega_ {d} \rightarrow \mathbb{R}^2$ and the fluid pressure $p_ {d}:\Omega_ {d} \rightarrow \mathbb{R}$ solution of (2.9a)–(2.9d), that we recall here,

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

● Navier boundary conditions:

● Neumann boundary conditions:

● Robin boundary conditions:

● Darcy boundary conditions without correction term

In what follows, we shall denote the weak solutions of systems (2.9a)–(2.9d) with boundary conditions of type Navier (3.19), Neumann (3.20), Robin (3.21) or Darcy without correction (3.22) by

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.22), the well-posedness of the solution can be proved by using the same arguments as for system (2.9) in Section 3.1.

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

and the divergence-free subspace $V_ {d} \overset{ \smash{ \mathrm{def} }}{ = } \left\{ \mathit{\boldsymbol{u}}\in H^1(\Omega_ {d})^2: \ \operatorname{{div}} \mathit{\boldsymbol{u}} = 0\right\}$.

3.2.2.   Main result

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

Theorem 2 (Asymptotic analysis). Let ${d} > 0$ be given and let $(\mathit{\boldsymbol{u}}_{d}, \hat{p}_{d}, p_ {d})\in \mathcal{U}_ {d}\times \mathcal{D}\times L^2(\Omega_ {d})$ be the solution to (2.9) which depends on the conductivity parameters $\kappa$ and $\eta$. The following propositions hold true:

$(i)$ When either $\kappa$ or $\eta$ goes to zero, $(\mathit{\boldsymbol{u}}_{d}, p_ {d}-\frac{1}{|\Omega_ {d}|} {\int}_{\Omega_ {d}} p_ {d})$ weakly converges towards $(\mathit{\boldsymbol{u}}_ {d}^{Na}, p_ {d}^{Na})$ in $H^1(\Omega_ {d})\times L^2(\Omega_g)$ solution to (2.9a)–(2.9d) with Navier boundary conditions (3.19) on $\Sigma$;

$(ii)$ When both $\kappa$ and $\eta$ go to infinity, $(\mathit{\boldsymbol{u}}_{d}, p_ {d})$ weakly converges towards $(\mathit{\boldsymbol{u}}_{d}^{Ne}, p_ {d}^{Ne})$ in $H^1(\Omega_ {d})^2\times L^2(\Omega_g)$ solution to (2.9a)–(2.9d) with Neumann boundary conditions (3.20) on $\Sigma$;

$(iii)$ When $\kappa$ fixed and $\eta\to \infty$, $(\mathit{\boldsymbol{u}}_{d}, p_ {d})$ weakly converges towards $(\mathit{\boldsymbol{u}}_{d}^{R}, p_ {d}^{R})$ in $H^1(\Omega_ {d})^2\times L^2(\Omega_g)$ solution to (2.9a)–(2.9d) with Robin boundary conditions (3.21) on $\Sigma$;

$(iv)$ When $\kappa\to \infty$ and $\eta$ fixed, $(\mathit{\boldsymbol{u}}_{d}, p_ {d}, \hat{p}_{d})$ weakly converges towards $(\mathit{\boldsymbol{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 to (2.9a)–(2.9d) coupled with Darcy without correction system (3.22) on $\Sigma$.

Remark 5. Theorem 2 mainly focuses on the asymptotic behavior of the pure fluid part of (2.9). Some asymptotic results for $\hat p$ in the limit cases mentioned above are discussed in Remark 6. In particular, it is worth mentioning that, in the different cases studied in Theorem 2 (except in case (iv)), the fluid unknowns $(\mathit{\boldsymbol{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 $\mathit{\boldsymbol{u}}_{ \kappa, \eta}$, $p_{ \kappa, \eta}$ and $\hat{p}_{ \kappa, \eta}$ that make the dependence on the conductivity parameters explicit. Thanks to a priori estimate (3.3) on $\mathit{\boldsymbol{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 $\mathit{\boldsymbol{u}}_{ \kappa, \eta}$ the subsequence and by $\mathit{\boldsymbol{u}}_l$ its limit:

By continuity of the trace and of the divergence operators, we have $\mathit{\boldsymbol{u}}_l\in V_ {d} \cap \mathcal{U}_ {d}$. We now consider the different cases.

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

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

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

which implies that $\mathit{\boldsymbol{u}}_l\cdot {\mathit{\boldsymbol{n}}}$ is constant on $\Sigma$. Morever, since $\mathit{\boldsymbol{u}}_l\cdot{\mathit{\boldsymbol{n}}}\in H^{\frac{1}{2}}_{00}(\Sigma)$ (see Remark 4), we recover $\mathit{\boldsymbol{u}}_l\cdot {\mathit{\boldsymbol{n}}} = 0$ also in this case. We therefore have that $\mathit{\boldsymbol{u}}_l\in V_ {d}\cap \mathcal{U}_ {d}^{Na}$ either when $\kappa$ or $\eta$ goes to zero.

Testing (3.2) with $\mathit{\boldsymbol{v}} \in \mathcal{U}_ {d}^{0, Na}\cap V_ {d}$ and both $\hat{q} = 0$ and $q = 0$, yields

so that, passing to the limit,

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

For the pressure, we consider its zero mean part

which, according to Proposition 1, is bounded independently of $\kappa$ and $\eta$. We then have a subsequence $(\mathit{\boldsymbol{u}}_{ \kappa, \eta}, p^*_{ \kappa, \eta})$ in $\mathcal{U}_ {d} \times L^2_0(\Omega_g)$ which weakly converges towards $(\mathit{\boldsymbol{u}}^{Na}, p^*_l)$. Testing (3.2) with $\mathit{\boldsymbol{v}} \in \mathcal{U}_ {d}^{0, Na}$ and both $\hat{q} = 0$ and $q = 0$, yields

which, passing to the limit, gives

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

This concludes the proof of (i).

Proof of (ii) and (iii). When neither $\eta$ nor $\kappa$ vanish, we have control on $\mathit{\boldsymbol{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

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

which is exactly the weak formulation associated to Robin boundary conditions (3.21). For $\kappa$ fixed, we therefore recover, when $\eta \to \infty$, the solution to the problem with Robin boundary condition, which is unique, i.e.,

When both $\eta\to\infty$ and $\kappa\to\infty$, we have that the term $\frac{\varepsilon_{\rm p}}{4 \kappa}(\mathit{\boldsymbol{u}}_{ \kappa, \eta} \cdot {\mathit{\boldsymbol{n}}}, \mathit{\boldsymbol{v}} \cdot {\mathit{\boldsymbol{n}}})_{0, \Sigma}$ disappears, so that we retrieve the Neumann boundary conditions at the limit and have

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{\varepsilon_{\rm p}}{4 \kappa}(\mathit{\boldsymbol{u}}_{ \kappa, \eta} \cdot {\mathit{\boldsymbol{n}}}, \mathit{\boldsymbol{v}} \cdot {\mathit{\boldsymbol{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 1. This concludes the proof. □

Remark 6. The following asymptotic behaviors for the Darcy pressure can also be proved:

Indeed, from (3.3), we have $\vert \hat{p}_{d}\vert_{1, \Sigma}^2\leqslant \frac{C_1}{\varepsilon_{\rm p} \eta}$, so that $\hat{p}_{d}$ strongly convergences 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

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

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

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.23) and (3.24), so that

By density arguments, we conclude that

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

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}$ converge at the same speed, $\hat{p}_{d}$ is bounded in $H^1$ and therefore admits a subsequence that weakly converges 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 ${\mathit{\boldsymbol{u}}}\cdot{\mathit{\boldsymbol{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.7), enables contact even with Dirichlet boundary conditions on the disk. The proof of this result is based on an estimate (Theorem 3 below) of the asymptotic behavior of the drag force (2.10) 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 ^[12,17]. For Navier or Dirichlet boundary conditions, similar arguments have previously shown that the drag force becomes singular as ${d} \to 0$ (see, e.g., ^[17,29] in 2D and ^[13,19] in 3D). In some cases, such as Dirichlet boundary conditions on both the disk and the wall or the combination of Navier and Dirichlet boundary conditions, the singularity in the drag force scales as $d^{-\frac32}$ 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.9) is bounded as $d\to0$.

Theorem 3 (Drag force estimate). The drag force $\mathcal F_ {d}$ given by (2.10) is nonnegative and bounded from above independently of ${d}$. More precisely, we have the following estimate

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

Corollary 1 (Collision result). Assume that Eq (2.11), 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 <- C {d}_0$ we have $t_c<+\infty$ with ${d}(t_c) = 0$ and $\dot {d}(t_c)<0$.

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

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

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

Assume that, at time $t = 0$, the disk is located at position ${\mathit{\boldsymbol{x}}}_0 = {d}_0 {\mathit{\boldsymbol{e}}}_2$ and has the initial velocity $\mathit{\boldsymbol{\dot x}}_G = \dot {d}_0 \mathit{\boldsymbol{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\leqslant t\leqslant t_1$, so that from (2.11) and (3.27) we get

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

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

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

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

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

Remark 7. Note that in the current model nothing prevents the disk of penetrating the bottom wall. Indeed, if the initial velocity $\vert \dot {d}_0\vert$ 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.7). 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

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. 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.26). Since the proof builds on ^[17], 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

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

Definition 1 (Operator $\operatorname{A}$) For any $f\in L^2_0(\Sigma)$ the problem

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

In particular, for all $\mathit{\boldsymbol{v}}\in \mathcal{A}_ {d}$, $\left. \mathit{\boldsymbol{v}}\cdot{\mathit{\boldsymbol{n}}}\right|_{\Sigma}\in L^2_0(\Sigma)$ as $\mathit{\boldsymbol{v}}$ is divergence free. We have that $\operatorname{A} (\mathit{\boldsymbol{v}} \cdot {\mathit{\boldsymbol{n}}})$ is the solution of the following strong form of problem (3.32) with $f = \mathit{\boldsymbol{v}}\cdot{\mathit{\boldsymbol{n}}}$:

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

Proposition 2 (Energy minimization problem). We introduce the energy functional $\mathcal{E}_ {d}: \mathcal{A}_ {d}\rightarrow \mathbb{R}^+$ defined by

for all $\mathit{\boldsymbol{u}} \in \mathcal{A}_ {d}$. The drag force (2.10) satisfies

Moreover, the energy minimization problem (3.35) admits a unique solution in $\mathcal{A}_{{d}}$, $\mathit{\boldsymbol{u}}_{d}$, given by (3.2).

Proof. We first prove that

Owing to the weak formulation (3.2), we have

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

Since $\operatorname{{\mathit{\boldsymbol{d}}}iv} (\boldsymbol{\sigma}(\mathit{\boldsymbol{u}}_{d}, p_ {d}))\in L^2(\Omega_ {d})$, the trace of $\boldsymbol{\sigma}(\mathit{\boldsymbol{u}}_{d}, p_ {d}) {\mathit{\boldsymbol{n}}}$ can be defined by duality as follows:

for all $\xi\in H^{\frac12}(\partial \Omega_ {d})$, where the symbol ${\mathcal{L}}$ denotes a standard lifting operator from $H^{\frac12}(\Sigma)^2$ to $H^1(\Omega_ {d})^2$. Taking ${\mathcal{L}} \boldsymbol \xi = \mathit{\boldsymbol{u}}_{d}$ and using that $\mathit{\boldsymbol{u}}_{d}\in \mathcal{U}_ {d}\cap V_d$ and the free sliding condition in (2.5), we get

where the drag force appears in the term on $\partial S_ {d}$. It is standard to check that $\boldsymbol{\sigma}(\mathit{\boldsymbol{u}}_ {d}, p_ {d}) {\mathit{\boldsymbol{n}}}\in L^2(\partial S_ {d})$. Indeed, the $\mathcal{C}^\infty$ regularity of the disk and the Dirichlet boundary condition on $\partial S_ {d}$ guarantee $H^2$ regularity for the fluid velocity $\mathit{\boldsymbol{u}}_{d}$ and $H^1$ regularity for the pressure $p_ {d}$ near the ball (see, e.g., ^[6]), so that we have

As regards the first term in (3.37), we use the relation (3.17), the free sliding condition (2.5)₃ and the fact that $\mathit{\boldsymbol{u}}_{d}\vert_\Gamma = \boldsymbol 0$ to get

Since by construction $\hat{p}_{d} = \operatorname A(\mathit{\boldsymbol{u}}_{d} \cdot {\mathit{\boldsymbol{n}}})$, taking $\hat{q} = A(\mathit{\boldsymbol{u}}_{d} \cdot {\mathit{\boldsymbol{n}}})$ in (3.32) yields

so that

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

which yields (3.36).

In order to show that $\mathit{\boldsymbol{u}}_{d}$ is a minimizer of the functional $\mathcal{E}_ {d}$, we introduce the tangent Hilbert space to $\mathcal{A}_ {d}$, namely:

Straightforward computations yield

On the other hand, from the definition of $\operatorname*{A}$ (3.32), the following symmetry relation holds

Therefore, since $\hat{p}_{d} = \operatorname{A}(\mathit{\boldsymbol{u}}_{d}\cdot{\mathit{\boldsymbol{n}}})$ and using the definition of the bi-linear form $a$ in (3.10), we have that

and, by inserting this identity into (3.41), we get

Yet, testing (3.2) with $\mathit{\boldsymbol{v}} \in \mathcal{A}_ {d}^0 \subset \mathcal{U}_ {d}^0$, $q = 0$ and $\hat q = \operatorname{A}(\mathit{\boldsymbol{v}}\cdot{\mathit{\boldsymbol{n}}})$, gives

so that, from (3.42), we finally have

which, since $\mathcal{E}_ {d}(\mathit{\boldsymbol{v}})\geqslant 0$, yields (3.35).

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

and evaluate the second Fréchet differential of $J_ {d}$ to show that

for all non zero $\mathit{\boldsymbol{v}}\in \mathcal{A}_ {d}^0$ and ${\mathit{\boldsymbol{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.35) 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.

Lemma 1 (Stream function). The admissible space $\mathcal{A}_ {d}$, defined in (3.31), is also given by

Moreover, we have

for all $\mathit{\boldsymbol{u}} \in \mathcal{A}_ {d}.$

Proof. We start by proving (3.43). Given $\mathit{\boldsymbol{u}}\in \mathcal{A}_ {d}$, $\mathit{\boldsymbol{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 $\mathit{\boldsymbol{u}} = -\partial_2\varphi {\mathit{\boldsymbol{e}}}_1 + \partial_1 \varphi{\mathit{\boldsymbol{e}}}_2$ (see, e.g., [14, Theorem 3.1]). In (3.43), 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 $\mathit{\boldsymbol{u}}\vert_{\partial S_ {d}} = {\mathit{\boldsymbol{e}}}_2, \quad \mathit{\boldsymbol{u}}\vert_{\Gamma} = \boldsymbol 0,$ enforce

Conversely, given $\varphi\in H^2(\Omega_ {d})$, if $\varphi$ satisfies (3.45) then obviously $\mathit{\boldsymbol{u}} = -\partial_2 \varphi {\mathit{\boldsymbol{e}}}_1 +\partial_1 \varphi {\mathit{\boldsymbol{e}}}_2$ belongs to $\mathcal{A}_ {d}$.

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

Integrating the first relation yields

On the other hand, using (3.46)₂, we have

As a result, since from (3.43) we have $\varphi(- {l}, 0) = 0$, we get

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.35).

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

for all ${d}>0$ and $0<\delta<1$, and where $\gamma : x \in [-\delta, \delta] \longmapsto 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

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

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., ^[12,17,29]) 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 space between the disk and the bottom wall. The terms involving $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.4.

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.,

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

Lemma 2 (Symmetry of the solution). Solution of (2.9) satisfies the following symmetries:

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

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

Lemma 3. Let $\tilde{\mathcal{A}}_{{d}}$ be the relaxed space given by (3.50). There holds

Proof. Let $\mathit{\boldsymbol{u}} = -\partial_2 \varphi {\mathit{\boldsymbol{e}}}_1 +\partial_1 \varphi {\mathit{\boldsymbol{e}}}_2\in \mathcal{A}_ {d}$. The boundary conditions $\mathit{\boldsymbol{u}}\vert_{{\partial S}} = {\mathit{\boldsymbol{e}}}_2$ implies

from which we deduce

so that

As we already fixed the constant by imposing $\varphi({- {l}, 0}) = 0$ in (3.43), the constant $c$ cannot be arbitrarily chosen. We have $c = \varphi(0, {d})$. Owing to Lemma 2, by a symmetry argument on $\mathit{\boldsymbol{u}}$, we have $u_1(0, y) = 0$ in $\left(0, {d}\right)$, so that

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

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

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}$

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.51). This completes the proof. □

The relaxed minimization problem reads therefore as follows: Find

such that

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

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

Proposition 3 (Explicit expression of $\tilde{\varphi}_{{{d}}}$). The relaxed minimization problem (3.52) admits a unique solution $\tilde {\mathit{\boldsymbol{u}}}_{{{d}}} = -\partial_2 \tilde{\varphi}_d {\mathit{\boldsymbol{e}}}_1 + \partial_1 {\tilde{\varphi}_d} {\mathit{\boldsymbol{e}}}_2 \in \tilde{\mathcal{A}}_{{d}}$, given by

with

Proof. As in the proof of Proposition 2, we can rewrite (3.52) as a minimization problem in the tangent Hilbert space $\tilde{\mathcal{A}}_{{d}}^0$:

So, let $\mathit{\boldsymbol{u}}_0 \in \tilde{\mathcal{A}}_{{d}}$ be given, instead of (3.52), we consider the auxiliary equivalent problem

Let us denote $\mathfrak{D}^2\tilde J_ {d}(\mathit{\boldsymbol{u}})$ the second Fréchet differential of $\tilde J_ {d}$ in $\mathit{\boldsymbol{u}}$. For all $\mathit{\boldsymbol{u}}\in \tilde{\mathcal{A}}_{{d}}^0$, we have

Let $\mathit{\boldsymbol{v}}\in \tilde{\mathcal{A}}_{{d}}^0$ and $\psi$ its associated streamfunction. The relation ${\mathfrak{D}}^2\tilde J_ {d}(\mathit{\boldsymbol{u}})(\mathit{\boldsymbol{v}}, \mathit{\boldsymbol{v}}) = 0$ implies that

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

The minimum $\tilde {\mathit{\boldsymbol{u}}}_{{{d}}} \in \tilde{\mathcal{A}}_{{d}}$ of (3.52) can be characterized by the Euler equations associated to the auxiliary minimization problem (3.56), which writes ${\mathfrak{D}}J_ {d}(\tilde {\mathit{\boldsymbol{u}}}_{{{d}}}- \mathit{\boldsymbol{u}}_0) \mathit{\boldsymbol{v}} = 0$ for all $\mathit{\boldsymbol{v}} \in \tilde{\mathcal{A}}_{{d}}^0$, or equivalently

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

After two integrations by parts, we obtain

From (3.51) and (3.57) we obtain, using standard density and trace arguments, the following PDE system for $\tilde{\varphi}_d$

whose solution is given by (3.54). Indeed, the streamfunction ${\tilde{\varphi}_d}$ defined by (3.54) is a polynomial of 3$^{\rm rd}$ degree in $y$, so that it satisfies (3.58)₁. For the second relation in (3.58), as

one can check that $\partial_2^3 \tilde{\varphi}_d(x, 0) + \frac{1}{\varepsilon_p \eta }\tilde{\varphi}_d(x, 0) = \frac{6x}{6+2\alpha_d(x)}\left(-\frac{ \alpha_d(x)}{(d+\gamma(x))^{3}} +\frac{1}{\varepsilon_p \eta}\right) = 0.$ The rest of the computations are straightforward. The solution (3.54) therefore characterizes the unique minimizer $\tilde{ \mathit{\boldsymbol{u}}}_{{d}}$ of (3.52). This completes the proof. □

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

Theorem 4. The relaxed drag force (3.53) is given by

In particular, we have that

Proof. In order to prove (3.59), we combine the definition (3.53) of the relaxed drag force $\tilde{\mathcal{F}}_ {d}$ with the explicit expression of the minimizer given in terms of the stream function $\tilde{\varphi}_d$ in (3.54). Straightforward computations then yield

Now, by inserting into (3.59) the expression for $\alpha_ {d}$ given in (3.55), we get

Since $\lim _{{d} \rightarrow 0} f_{{d}}(x) = {1/\big({3\varepsilon_{\rm p} \eta +\gamma(x)^{3}}\big)}$, the dominated convergence theorem gives

which is clearly bounded when $\varepsilon_{\rm p}\eta$ is fixed. Indeed, as $\gamma(x)\geqslant 0$ for any $x\in [-\delta, \delta]$, we have

which yields (3.60) and completes the proof. □

We conclude this section with a series of remarks.

Remark 8 (Darcy versus Navier or Dirichlet). A salient difference 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., $\mathit{\boldsymbol{u}}\cdot {\mathit{\boldsymbol{n}}}$ can be non-zero on $\Sigma$, so that $\partial_1 \varphi$ can play an important role in the estimation of the drag force. On the contrary, the work reported in ^[12] only focuses on the horizontal velocity of the fluid. Note that we neglect the impact of $\kappa$ in the relaxed energy (3.49) in order to facilitate the solution of (3.52) (Proposition 3). However, the dependance on $\kappa$ is retrieved 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^2_1 \varphi$ cannot be simplified in the energy functional (3.44), so that we do not have $\tilde{\mathcal{F}}_ {d} \leqslant \mathcal{F}_ {d}$ as in the Navier or Dirichlet case (see ^[12]).

Remark 9 (Asymptotic of $\tilde{\mathcal{F}}_ {d}$ for $\varepsilon_{\rm p} \eta = \mathcal O({d}^p)$, $p\geqslant 3$). If $\varepsilon_{\rm p} \eta = \mathcal O({d}^p)$ with $p\geqslant 3$, we obtain an asymptotic behavior for the relaxed drag force (3.59) similar to the case with Dirichlet boundary conditions on both the disk and the wall

Indeed, by using $\varepsilon_{\rm p} \eta = \mathcal O({d}^p)$, $p\geqslant 3$, in (3.55) we infer that $\alpha_ {d}\geqslant 1$. Since $f : x \mapsto \frac{x}{x+3}$ is continuous in $[1, +\infty)$ and $\lim_{x\to \infty} f(x) = 1$, we have

from which we deduce

We therefore recover Dirichlet asymptotic behavior as ${\int}_{-\delta}^\delta \frac{x^{2}}{({d}+\gamma(x))^{3}} \mathrm d x \sim {d}^{-\frac32}$. 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 $\mathit{\boldsymbol{u}} = \frac{{\mathit{\boldsymbol{v}}}}{\sqrt{{d}}}$ to exhibit ${d}^{-\frac32}$ out of the integral (see ^[12] for the no-slip case).

Note that the obtained asymptotic (3.61) of the relaxed drag force is consistent with the asymptotic analysis carried out at the Section 3.2. Indeed, as stated in Theorem 2, when $\eta$ tends to zero the solution of problem (2.9) converges weakly towards solution of a problem with Navier boundary conditions on the wall and Dirichlet boundary conditions on the disk. So, for $\varepsilon_{\rm p} \eta = \mathcal O({d}^p)$, $p\geqslant 3$ we expect to recover (3.61) 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., ^[13]). Yet, the proved asymptotic behavior (3.61) only concerns the relaxed drag force $\tilde{\mathcal{F}}_ {d}$. Nothing can be said on the asymptotic of the whole drag force due to the lack of lower bound on $\mathcal{F}_ {d}$ (see Remark 8).

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

by combining the Taylor expansion $\gamma(x) = \frac{x^2}{2} + O_{x\to 0}(x^4)$ with the change of variable $t = x/d^{\frac{p}{6}}$. This asymptotic behavior could be useful in a more general setting, in which the compressibility of the porous matrix can induce a change in the conductivity parameters. For instance, when the ball gets closer to the wall, it compresses the porous matrix and thus might reduce its conductivity.

3.3.4.   Proof of Theorem 3

The main idea to complete the proof consists in extending $\tilde{\varphi}_d$ to the whole domain $\Omega_ {d}$ in order to obtain an admissible velocity $\check{ \mathit{\boldsymbol{u}}}_ {d}\in \mathcal{A}_ {d}$. To obtain a bound on $\mathcal{F}_ {d}$, we have to estimate the whole energy $\mathcal{E}_ {d}(\hat{ \mathit{\boldsymbol{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 ^[17] and use similar notations for the auxiliary functions.

Proof. We consider the following extension of $\tilde{\varphi}_d$ (see ^[17]):

where the auxiliary function $\chi_{2\delta}$ is defined as

with

and the remaining auxiliary function $\hat{\varphi}_ {d}$ is given by

with ${\mathit{\boldsymbol{x}}}_ {d} \overset{ \smash{ \mathrm{def} }}{ = } (0, 1+ {d})$ denoting the center of the disk. The auxiliary function $\chi_{\delta}$ has a smooth $\mathcal{C}^2$ decreasing transition between $1+\frac{\delta}{2}$ and $1+\delta$. Figure 4 illustrates how $\check{\varphi}_ {d}$ is contructed by combining $\hat\varphi_ {d}$ and $\tilde{\varphi}_d$. The different zones have been depicted in Figure 5 for illustration purposes.

Let us introduce the function

We first check that $\check{ \mathit{\boldsymbol{u}}}_ {d}$ is admissible, i.e., that $\check{ \mathit{\boldsymbol{u}}}_ {d}\in \mathcal{A}_ {d}$, with $\mathcal{A}_ {d}$ the admissible space defined in (3.31). As $\chi_\delta$ belongs to $\mathcal{C^\infty}((1, 1+\delta))$ by construction, the function $\hat{\varphi}_ {d}$ also belongs to $\mathcal{C}^\infty(\Omega_ {d})$. The extension $\check{\varphi}_ {d}$ has therefore the same $H^2$-regularity than $\tilde{\varphi}_d$, so that $\check{ \mathit{\boldsymbol{u}}}_ {d}$ belongs to $H^1(\Omega_ {d})$. Then, one can easily check that $\check{{\mathit{\boldsymbol{u}}}}_ {d} = {\mathit{\boldsymbol{e}}}_2$ on $\partial S_ {d}$.

Finally, as the support of $\check\varphi_ {d}$ is strictly included in $\Omega_ {d}$, we do have $\check{ \mathit{\boldsymbol{u}}}_ {d} = {\mathit{\boldsymbol{0}}}$ on the boundary $\Omega_ {d}\setminus(\Sigma\cup\partial S_d)$. According to Proposition 2, since $\check{ \mathit{\boldsymbol{u}}}_ {d}\in \mathcal{A}_ {d}$, we have the following bound on $\mathcal{F}_ {d}$

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

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

Under the disk. We consider the following energy:

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

with $\alpha_ {d}$ defined in (3.55) and the auxiliary functions

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

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

We finally obtain the estimate of the first term

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

Finally, we obtain the following estimate

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

Away from the disk $\Omega_ {d}\setminus\Omega_{{d}, 2\delta}$. Outside of $\Omega_{{d}, 2\delta}$, $\check{\varphi}_ {d} = \hat{\varphi}_ {d}$ and $\hat{\varphi}_ {d}$ depends on the gap distance ${d}$ only through a translation (see ^[17]HY__HY, Proof of Lemma 8]). Thus, energy in $\Omega \diagdown \Omega_{{d}, 2\delta}$ of $\check{ \mathit{\boldsymbol{u}}}_ {d} = \nabla^\perp \check{\varphi}_ {d}$ does not depend on ${d}$ and is therefore bounded. So, there exists a constant $C_9(\delta) > 0$ such that

for ${d} < {\tilde l}$.

Transitional zone $\Omega_{{d}, 2\delta}\setminus\Omega_{{d}, \delta}$. Since $\vert\chi_\delta\vert\leqslant 1$, using Young's inequality, we have

So, the estimate of the transitional term can simply be obtained by combing the bounds (3.66) and (3.67).

Finally, (3.26) then follows by combining the previous bounds with (3.60), which concludes the proof of Theorem 3. □

Remark 11. The analysis above can be adapted straightforwardly to the case in which the pure slip condition (2.5)₃ is replaced by the BJS law (2.6). Indeed, in this case it suffices to add the following tangential friction term

to the relaxed energy (3.49). By building on similar arguments than those used for the derivation of (3.59), one can show that the associated drag force can be bounded by $C(\varepsilon_{\rm p}\eta)^{-1}$, with $C > 0$ a constant independent of the friction coefficient $\alpha$. Similarly, the extension contributions can also be bounded irrespectively of $\alpha$. In particular, this means that the drag force does not blows when $\alpha\to\infty$, as ${\mathit{\boldsymbol{u}}}\cdot\boldsymbol{\tau} \rightarrow 0$ and the fluid can escape through the porous layer along the normal direction.

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.9), for different values of $d$ (viz., the distance between the disk and the 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, and the weak convergence results stated in Theorem 2 when either $\kappa$ or $\eta$ vanishes. Then, we numerically investigate the ability of the porous layer to enable collision as stated Corollary 1. Finally, we provide numerical insight on the validity of the fluid-reduced Darcy model (2.9) 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.9) 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.9a)–(2.9d) supplemented with Navier boundary conditions (3.19). The values for ${d}$ are chosen to be equidistant in log scale

This experiment setup is very similar to the one considered in ^[18]. Therein the authors retrieve the $d^{-\frac32}$ 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.9). 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}^f_ {d}$ denote a triangulation of the fluid domain $\Omega_ {d}$ and let $\mathcal{T}^p_ {d}$ be the 1D triangulation of the porous mid-surface $\Sigma$, given simply as the edges of $\mathcal{T}^f_ {d}$ 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,

and the space of continuous piecewise linear functions for the fluid pressure

We also introduce the homogeneous counterpart of $\mathcal{U}_{{d}, h}^{Na}$

For any ${d} > 0$, the corresponding finite element approximation of (2.9a)–(2.9d) with the Navier boundary conditions writes

Here, the small perturbation of the incompressibility constraint $\varepsilon_s\left(p_{h}, q_{h}\right)_{\Omega_{{d}}}$, with $\varepsilon_s = 10^{-10}$, serves to enforce zero mean on the discrete pressure $p_h$ (this can be inferred by taking $(\mathit{\boldsymbol{v}}_h, q_h) = (\boldsymbol 0, 1)$ in (4.2) and by integrating by parts in the divergence term). This is a standard approach for fixing the constant of the pressure for Stokes problems with only Dirichlet boundary conditions (see, e.g., ^[4,16]).

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$,

with its homogeneous counterpart

For the Darcy pressure, we choose the space of continuous piecewise polynomial functions of degree 2 on $\Sigma$, namely,

The finite element approximation of (3.2) finally writes

The last term enforces the mean value of the Darcy pressure $\hat{p}_{h}$ to be 0 by penalization, 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

The meshes are not structured uniformly 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 ^[16]) integrated mesh generator then we optimize them with the Mmg remeshing software (see ^[8]). The resulting meshes have between $10^4$ and $10^5$ nodes (see Figure 6).

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 2 and the estimate of the drag force stated in Theorem 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 surface porous model (2.4)-(2.5) is derived in ^[7] by averaging the original Darcy system across the thickness. In order to investigate the validity of this approximation, we simulate the fall of the disk over the 2D porous layer $\Omega_{\rm p} \overset{ \smash{ \mathrm{def} }}{ = } (- {l}, {l}) \times (0, -\varepsilon_{\rm p})$. The porous fluid velocity in $\Omega_{\rm p}$ is given by the Darcy relation ${\mathit{\boldsymbol{u}}}_{\rm p} \overset{ \smash{ \mathrm{def} }}{ = } -K\boldsymbol{\nabla} p_{\rm p}$ (see, e.g., ^[2,25]), where the porous fluid pressure $p_{\rm p}:\Omega_{\rm p}\rightarrow \mathbb{R}$ is described by the divergence free constraint:

with the conductivity matrix $K \overset{ \smash{ \mathrm{def} }}{ = } \mbox{diag}(\kappa, \eta)$. The coupling between (4.4) and the Stokes system (2.9a)–(2.9d) is then operated by the following kinematic and dynamic interface conditions on $\Sigma$:

Figures 7 and 8 report the velocity and pressure fields obtained in both cases and indicate that the behavior of the solution of the 2D Darcy problem (4.4) is very similar to the solution of the Darcy surface model (2.4) considered in this paper. In both cases, the physical parameters are $\kappa = \eta = 1$, $\varepsilon_{\rm p} = 0.01$ and ${d} = 0.001$. In Figures 7(b) and 8(b), the porous layer thickness has been amplified by a factor 25 to better visualize the Darcy pressure.

Drag force asymptotics. In Figure 9, we report the drag force $\mathcal F_ {d}$ for ${d}\in[10^{-3}, 1]$ arising from system (2.9) where the fluid is coupled to the porous layer model. For the sake of simplicity, we call it the Darcy drag force. We also plot the drag force asymptotic arising from system (2.9a)–(2.9d) with Navier boundary conditions (3.19), that we similarly call Navier drag force asymptotic.

We test various set of values for the conductivity parameters $\kappa$ and $\eta$, notably to investigate the behavior when the conductivity parameters gets smaller. Thus, we can compare with the asymptotic results of Theorem 2 and the estimate of the Darcy drag force (3.26) stated in Theorem 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, whereas for Navier the drag force explodes. The very good recovery of the theoretical -1.5 slope for Navier drag force is clearly visible.

As expected from the estimate on the drag force (3.26), the plateau value of the Darcy asymptotic increases when $\kappa$ and/or $\eta$ goes to 0. We can also note that the conductivity parameters $\kappa$ and $\eta$ play a role in the switch from the Navier 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 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 10 and 11 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$).

In the case of the coupling with a porous layer, Figure 11, 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 10 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{\cdot}10^{-3}$ from the wall and we look at the evolution of the incoming fluid flow ${\mathit{\boldsymbol{u}}}_h \cdot {\mathit{\boldsymbol{n}}}$ on $\Sigma$ and of the Darcy pressure $\hat{p}_h$ when $\kappa$ or $\eta$ vanishes, see Figures 12 and 13. We zoom around the contact point at $x = 40$.

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 11. 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 ${\mathit{\boldsymbol{u}}}_h\cdot {\mathit{\boldsymbol{n}}}$ to 0, which is consistent with weak convergence results of Theorem 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 2, with a huge peak of the Darcy pressure up to $8{\cdot}10^3$ at $\eta = 10^{-6}$, see Figure 13. This singularity makes it difficult to run the simulation for smaller values of $\eta$ because of numerical instabilities. It also gives numerical evidence on why in Figure 9 the Darcy drag force assymptotic seems to converge much faster towards the Navier one for $\eta\to 0$ than $\kappa\to 0$. 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 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 12 and 13: in the first asymptotic we observe that the (reentrant) fluid velocity on both sides of the disk tends to flatten as $\kappa$ tends to 0, while for the second asymptotic, two narrow peaks start to appear before flattening out.

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

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.5). In the following, we fix the conductivity parameters $\kappa$ and $\eta$ to be equal to 1. In Figure 9, 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 $\vert \dot d_0\vert$ going from 1 to 100, see Figure 14. As expected from Corollary 1, we have collision for $\vert \dot d_0\vert\leqslant Cd_0$.

We now compare the effect of Darcy and Navier on the ball dynamics in Figure 15, which provides the trajectories obtained in both cases for different initial velocities. For $\dot d_0 = 1$, the trajectories are indistinguishable. The porous layer has no significant effect on the dynamics whenever the disk stays far away from the wall. 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 disk maintains a distance ${d} = 0.1$ from the wall. The porous layer completely modifies the disk dynamics at contact.

5.   Conclusions

In this paper, we have analyzed the contact capabilities of a fluid-structure interaction model with seepage reported in ^[7]. 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.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

MF and MV were partially supported by the French National Research Agency (ANR), through the SIMR project (ANR-19-CE45-0020). The authors thank the anonymous reviewers for their constructive comments that helped to improve the manuscript.

Conflict of interest

The authors declare no conflicts of interest.