David Gérard-Varet

Anne-Laure Dalibard

Filippo Santambrogio

Vincent Millot

Élie Bretin

Pablo Alvarez- Caudevilla

Yannick Privat

Florian Omnès

Je tiens tout d'abord à remercier David Gérard-Varet, qui a accepté de rédiger un rapport préliminaire sur ce mémoire et m'a adressé de nombreux conseils, encouragements et commentaires bienveillants sur mon travail. Je remercie également Benoît Merlet et Michael Vogelius qui m'ont fait l'honneur d'écrire un rapport sur mon manuscrit. Un grand merci à Simon Masnou d'avoir accepté le rôle de président du jury, ainsi qu'à Éric Bonnetier, Dorin Bucur et Anne-Laure Dalibard qui ont répondu avec enthousiasme à mon invitation à participer à la soutenance. Je tiens également à remercier Xavier Blanc pour l'aide qu'il m'a apportée concernant les démarches administratives liées à la préparation de mon habilitation.

Les travaux qui sont rassemblés dans ce mémoire n'auraient pas pu exister sans l'apport de tous mes collaborateurs. Je tiens en particulier à remercier Antoine Lemenant, qui m'a initié très tôt aux méthodes de champ de phase lors de mon arrivée au LJLL, et Francisco Suárez-Grau, dont les multiples séjours à Paris ont nourri des discussions fructueuses autour des effets de parois sur les fluides.

Préambule

Ce mémoire présente certains résultats que j'ai obtenus depuis la fin de ma thèse, et qui concernent principalement l'approximation variationnelle de problèmes de minimisation de longueur issus du calcul des variations, et l'étude de modèles d'écoulements visqueux au voisinage de parois provenant de la mécanique des fluides ou de la biologie. Ces deux thèmes ont en commun l'utilisation de l'analyse asymptotique, qui se traduit formellement par l'introduction d'un petit paramètre ε et d'une suite d'objets qui lui sont associés, et dont on cherche à comprendre et à justifier le comportement à la limite quand ε tend vers zéro. Cette démarche peut avoir différents objectifs, parmi lesquels :

• construire une suite de problèmes bien posés et résolubles en pratique, dont les solutions constituent une approximation régulière de la solution du problème initial. L'approche variationnelle par Γ-convergence en est un exemple ;

• simplifier la géométrie d'un modèle présentant des échelles spatiales différentes, tout en exprimant l'effet des échelles que l'on a filtrées, à l'aide d'un modèle effectif. Les méthodes d'homogénéisation telles que la convergence double-échelle entrent dans cette catégorie.

Dans ce mémoire, nous abordons les méthodes variationnelles sous l'angle de l'approximation par champ de phase du problème de Steiner, consistant à déterminer l'ensemble connexe 1-dimensionnel de longueur minimale connectant une famille finies de points dans le plan. Ce problème peut être conçu comme un modèle pour développer des méthodes systématiques de résolution de problèmes de minimisation de longueur sous contrainte de connexité.

Nous détaillons au Chapitre 2 des résultats théoriques d'approximation variationnelle basés sur un résultat de Γ-convergence du à Modica et Mortola, et sur l'utilisation de distances géodésiques pour pénaliser la disconnexité des compétiteurs. Nous développons également la mise en oeuvre numérique de ces méthodes d'approximation.

Les méthodes d'homogénéisation sont au coeur des démonstrations et éléments de preuves développés aux Chapitres 3, 4 et 5.

Le Chapitre 3 traite de modèles d'interaction entre une paroi solide et un fluide visqueux, pour des écoulements à très faible échelle, visant à expliquer et à quantifier le glissement observé expérimentalement pour ce type de configuration. Le premier modèle utilise pour cela la géométrie et les propriétés chimiques particulières de parois dites superhydrophobes ; le second s'appuie sur l'hypothèse que la viscosité du fluide n'est pas réellement constante dans tout le domaine de l'écoulement, mais qu'il existe une fine couche près de la paroi dans laquelle elle est drastiquement réduite.

Le Chapitre 4 est centré sur un modèle d'écoulement turbulent au voisinage d'une paroi, appartenant à une famille de modèles très utilisés en pratique : les modèles dits LES, pour Large Eddy Simulation. Ces derniers sont basés sur le filtrage des petites échelles spatiales et sur la modélisation des interactions entre ces dernières et les grandes échelles, c'est-à-dire celles qui sont effectivement simulées. Nous proposons dans ce chapitre une étude asymptotique et numérique du modèle de Smagorinky permettant la prise en compte de motifs rugueux sur la paroi. Notre approche repose sur une méthode d'homogénéisation dérivée de la convergence double échelle, appelée unfolding method.

Le Chapitre 5 présente un modèle de lubrification par un type particulier de fluide complexe, appelé fluide micropolaire, pour lequel des conditions aux limites non standard ont été introduites assez récemment dans la littérature mécanicienne, puis formalisées mathématiquement. Nous nous intéressons dans ce chapitre à l'influence conjointe d'une géométrie rugueuse des parois et de ces conditions aux limites nouvelles.

Nous concluons ce mémoire par la présentation d'une méthode numérique d'optimisation de forme pour un modèle d'aquaporine (Chapitre 6). Il s'agit d'une sorte de canal nanoscopique issu du monde du vivant et qui joue un rôle prépondérant dans les échanges ioniques entre cellules. Nous mettons l'accent sur la mise en oeuvre d'une méthode de descente de gradient permettant de réduire la dissipation, et qui repose sur un procédé d'extension et régularisation de la déformation normale construite par un calcul de dérivée de forme.

Les chapitres sont généralement organisés de la manière suivante. Après une introduction présentant le contexte scientifique, nous détaillons le modèle et nos principaux résultats avant d'aborder les aspects des démonstrations qui nous paraissent les plus marquants. Nous proposons ensuite une étude numérique du modèle et concluons par une brève synthèse et quelques perspectives.

Publication list Preprints

[B1] M. Bonnivard, J. Olivier. From adherence to slip in nanofluidics: a mathematical justification based on a drop of viscosity. The scalar case.

[B2] M. Bonnivard, J. Olivier. From adherence to slip in nanofluidics: a mathematical justification based on a drop of viscosity. The 3d case.

[B3] M. Bonnivard, I. Pažanin, F. J. Suárez-Grau. A generalized Reynolds equation for micropolar flows past a ribbed surface with nonzero boundary conditions.

Published journal papers

[B4] M. Bonnivard, É. Bretin, A. Lemenant. Numerical approximation of the Steiner problem in dimension 2 and 3. Math. Comput., 89-321 (2020), 1-43.

[B5] P. Álvarez-Caudevilla, M. Bonnivard, A. Lemenant. Asymptotic limit of linear parabolic equations with spatio-temporal degenerated potentials. ESAIM: Control, Optimisation and Calculus of Variations, 26-50 (2020).

[B6] M. Bonnivard, V. Millot, A. Lemenant. On a phase field approximation of the planar Steiner problem: Existence, regularity, and asymptotic of minimizers. Interfaces Free Bound., 20-1 (2018), 69-106.

Chapter 1

Introduction

This thesis presents results that I have obtained since the end of my PhD, concerning two main topics: topic 1: the approximation of length minimization problems coming from calculus of variations;

topic 2: the study of fluid/structure interaction models from applied fluid mechanics.

Both topics have in common the frequent use of asymptotic analysis, either as a regularization and relaxation tool for the related minimization problems, or as a homogenization tool for extracting the main order behaviour of a model, with respect to a small parameter or geometric scale typically denoted ε. Even though the models or mathematical techniques may vary from one chapter to another, we have tried to present the results in a sort of unified manner, insisting on modeling issues as well as rigorous mathematical statements and efficient numerical simulations.

Let us now give an overview of the content of the next chapters.

What is contained in this thesis

Topic 1 is developed in Chapter 2, where we have synthesized results obtained in collaboration with É. Bretin, A. Lemenant, V. Millot and F. Santambrogio in [B11, B6, B4]. Our main contributions concern the variational approximation of length minimization problems with connectedness constraint, based on a phase field approximation and the penalization of unconnected competitors using geodesic distances.

To simplify the presentation, this chapter focuses on the Steiner problem. Starting from the original functional studied in [B11], we give the main ingredients of the proofs of convergence, that rely on Γ-convergence tools and geometric measure theory. Then, we present a relaxation of the functional introduced in [B6], that leads to a qualitative study of the numerical approximation of minimizers in [B4].

Topic 2 comprises Chapters 3 to 6. These chapters gather convergence results and numerical simulations of various models of viscous flows in the vicinity of a wall.

Chapter 3 deals with the justification by homogenization techniques of the slippage measured in micro-and nanofluidics, using two different models proposed by physicists:

• a model of rough hydrophobic surface;

• a depletion layer model. These models are applied to a similar setting: a laminar viscous flow in an infinite channel, in 3d.

The hydrophobicity model was studied in collaboration with D. Gérard-Varet and A.-L. Dalibard in [B12]. Our main contributions are the obtention of effective slip lengths and their validation using numerical simulations and comparisons with the existing literature.

The depletion layer model was addressed in a recent work with J. Olivier [B1, B2]. This model assumes that the viscosity is dramatically reduced in a small vicinity of solid walls, called the "depletion layer". Our main result is the identification of an effective slip, reached when passing to the limit in the size of this layer, using techniques from thin domain flow analysis.

Chapter 4 describes the main achievements of [B9] obtained with F. J. Suárez-Grau on the impact of surface roughness on Large Eddy Simulation of turbulent flow. The main results are the homogenization of a timedependent and non-linear model of turbulence due to Smagorinsky, using an adaptation of the two-scale convergence called the unfolding method, and the numerical simulation of the effective system in the particular case of riblets, which is of great interest in applications.

Chapter 5 focuses on results recently obtained in collaboration with I. Pažanin and F. J. Suárez-Grau [B3], on the lubrication with a micropolar fluid over a ribbed surface. Our main contributions are the obtention of a generalized Reynolds equation, which is a valuable tool in the fluid mechanics literature to compare the performances of lubricants, and a numerical study of the properties of the system depending on its most relevant mechanical and geometric parameters.

Finally, Chapter 6 presents results obtained in collaboration with F. Omnès and Y. Privat in [B8], on the shape optimization of a model of aquaporin. Aquaporins are proteins that play a crucial role in water transport between cells. In this chapter, we present a model initially proposed by soft matter physicists, and develop a numerical algorithm to optimize the dissipated energy by deforming the entrance and outrance regions of the channel. Our main contribution is an extension and regularization procedure on the normal deformation of these regions, allowing for a smooth deformation of the mesh during the gradient descent.

What is not contained in this thesis

Besides [B15, B14, B13] which were part of my PhD thesis, this habilitation thesis does not contain some results I obtained since then. Let me now quickly describe them.

In [B5], we studied in collaboration with P. Álvarez-Caudevila and A. Lemenant the asymptotic limit of parabolic problems in cylindrical domains, with degenerated potentials. We provide a strong convergence result for the solution by Γ-convergence techniques, and establish an exponential decay estimate. We decided not to include this article since the featured model seemed poorly related to the rest of the thesis.

[B7] addresses a model of micropolar lubrication, which is rather similar to the model described in Chapter 5, and is studied in the same spirit. Hence, we chose to omit this contribution to avoid redundancy.

Finally, [B10] is about the influence of wavy riblets on viscous flows, as opposed to classical straight riblets. Since the techniques of proofs are similar to the ones developed in Chapter 4, we found it necessary to choose between [B10] and [B9] to avoid redundancy.

General notation and conventions

General notation We denote by (e i ) 1≤i≤n the canonical basis of R n , with n = 2 or n = 3. H d stands for Hausdorff measure of dimension d in R n . When periodicity is involved, we write T = R/Z and T 2 = R 2 /Z 2 .

If Ω ⊂ R n is a Lipschitz domain, we denote by ν the outward-pointing normal to ∂Ω. The tangential part of a vector field ζ on ∂Ω is then defined by

[ζ] tan = ζ -(ζ • ν)ν .
We reserve bold letters for vector fields: for instance, if n = 3, we note u the vector field u = (u 1 , u 2 , u 3 ). In Chapters 3, 4 and 5, the variable x = (x 1 , x 2 , x 3 ) ∈ R 3 will be decomposed as x = (x , x 3 ) where x = (x 1 , x 2 ) is the horizontal variable. Analogously, any vector field u will be decomposed as u = (u , u 3 ) with u 3 = u • e 3 and u = u -u 3 e 3 , where (e 1 , e 2 , e 3 ) is the canonical basis of R 3 .

We denote O ε a sequence converging to zero as ε goes to zero.

Functional spaces Let U be an open subset of R n and k, m ∈ N * . We denote by

• C 1 c (U ) the space of C 1 functions with compact support in U ;

• H k (U, R m ) the Sobolev space W k,2 (U, R m ). For m = 1 we simply write H k (U ) ;

• H 1 0 (U ) the completion of C 1 c (U ) for the H 1 norm;

• H -1 (U ) the dual space of H 1 0 (U )

Differential operators

We write ∂ i u the i-th partial derivative of a function u. If functions depend on different variables, for instance x = (x 1 , x 2 , x 3 ) or y = (y 1 , y 2 , y 3 ), we sometimes use the notation ∂ x i u or ∂ y i v to highlight this dependence. If u = (u i ) 1≤i≤n is a vector field, we define:

• ∇u: the Jacobian matrix of u, defined by (∇u) i,j = ∂ j u i ,

• D(u): the symmetric part of ∇u, given by D(u) = ∇u+(∇u) T 2 ,

• ∇ × u: the curl of vector field u, i.e. (for n = 3)

∇ × u =    ∂ 2 u 3 -∂ 3 u 2 ∂ 3 u 1 -∂ 1 u 3 ∂ 1 u 2 -∂ 2 u 1    .
For every θ = (θ 1 , θ 2 ), we define

[θ ] ⊥ = (-θ 2 , θ 1 ), rot

x 3 θ = ∂ x 3 [θ ] ⊥ , Rot x θ = ∂ x 1 θ 2 -∂ x 2 θ 1 .
Chapter 2

Phase field approximation and numerical simulations of the Steiner problem

This chapter presents a variational approach to the approximation of length minimization problems among compact connected sets of Hausdorff dimension 1, and studies possible strategies to implement the method efficiently. It gathers results obtained in collaboration with É. Bretin, A. Lemenant, V. Millot and F. Santambrogio in [B11, B6, B4]. Although several applications are considered in these papers, we focus here on the Steiner tree problem, which is easy to formulate, well-understood from a theoretical point of view, but remains challenging from a numerical perspective.

A phase field approximation enforcing connectedness

In the pioneering work by Modica and Mortola [START_REF] Modica | Il limite nella Γ-convergenza di una famiglia di funzionali ellittici[END_REF] (see also [START_REF] Alberti | Variational models for phase transitions, an approach via Γ-convergence[END_REF]), it is proved that the functional

M ε (ϕ) = ε Ω |∇ϕ| 2 dx + 1 ε Ω ϕ 2 (1 -ϕ) 2 dx
converges in a certain sense to the perimeter functional

P (ϕ) = 1 3 Per(A, Ω) if ϕ = 1 A ∈ BV(Ω)
, and + ∞ otherwise.

9

The original motivation for Modica and Mortola was a mathematical justification of convergence for a phase transition model by Cahn and Hilliard. Later, this phase-field approximation based on the minimization of elliptic functional M ε was used by many authors in the last decades, with quite satisfactory results. Let us mention Ambrosio and Tortorelli [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF] for the Mumford-Shah functional from image segmentation [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associated variational problems[END_REF], Oudet [START_REF] Oudet | Approximation of partitions of least perimeter by γ-convergence: Around Kelvin's conjecture[END_REF] and Bourdin, Bucur and Oudet [START_REF] Bourdin | Optimal partitions for eigenvalues[END_REF] for partition problems, Oudet and Santambrogio [START_REF] Santambrogio | A Modica-Mortola approximation for branched transport[END_REF][START_REF] Oudet | A Modica-Mortola approximation for branched transport and applications[END_REF] for branched transport, among many others.

The Steiner tree problem is one of the simplest problems of this kind. Given a bounded and convex open set Ω ⊂ R 2 , with smooth boundary, and a collection of points a 0 , a 1 , . . . , a N ∈ Ω, it consists in finding a compact connected set K ⊂ Ω containing all the a i 's and having minimal length. In other words, it amounts to solving the minimization problem min{H 1 (K); K ⊂ Ω, connected, and containing all points a i } ,

where H 1 (K) stands for the one-dimensional Hausdorff measure of K.

It is known that minimizers for (2.1) exist, need not be unique, and are trees composed by a finite number of segments joining with only triple junctions at 120 • . However, computing a minimizer is very hard. In fact, it is known as a (N P )-hard problem in combinatorial theory (see e.g. [START_REF] Karp | Reductibility among Combinatorial Problems[END_REF]). Finding efficient algorithms to compute an approximate solution is still an active research field in graph theory (see for instance [START_REF] Forte | Iterated local search algorithms for the Euclidean Steiner tree problem in n dimensions[END_REF] and references therein).

One major difficulty raised by the numerical approximation of minimizers of the Steiner tree problem, is enforcing the connectedness of the competitors. The strategy developed in [B11] consists in adding in the Modica-Mortola functional, a penalization term relying on the weighted geodesic distance, defined as follows.

Let Ω ⊂ R 2 . Given a nonnegative Borel measurable function ϕ : Ω → [0, ∞), we define the geodesic distance between two points a, b ∈ Ω relative to the metric ϕ to be

D(ϕ; a, b) := inf Γ:a b Γ ϕ dH 1 ∈ [0, +∞] , (2.2)
where Γ : a b means that Γ is a rectifiable curve in Ω of finite length connecting a and b (i.e. a Lipschitz image of [0, 1] contained in Ω running from a to b).

From a numerical point of view, the geodesic distance can be computing using the so-called "fast marching method" [93]. This algorithm was improved in [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF] to compute at the same time the geodesic distance and its gradient with respect to the metric, which is especially useful when one wants to optimize the metric itself.

We propose to approximate problem (2.1) by minimizing the functional

1 4ε Ω (1 -ϕ) 2 dx + ε Ω |∇ϕ| 2 dx + 1 λ ε N i=1 D(ϕ; a 0 , a i ) + ε q+1 Ω |∇ϕ| q (2.3) among functions ϕ ∈ W 1,q (Ω) satisfying λ 2 ε ≤ ϕ ≤ 1 and ϕ = 1 on ∂Ω.
Here λ ε is a parameter satisfying λ ε → 0 as ε → 0, and q > 2 is a fixed exponent.

The key point which enforces connectedness is that whenever

N i=1 D(ϕ; a 0 , a i ) = 0 ,
the set {ϕ = 0} must be arcwise connected and contain the {a i }. Heuristically, last condition results from bounding the third term 1 λε N i=1 D(ϕ; a 0 , a i ). The first two terms are a simple variant of the standard Modica-Mortola functional, already used in [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF]. The very last term is only needed to guarantee the continuity of ϕ and easily prove the existence of a minimizer, but the power of ε is chosen so that it disappears at the limit.

Main result

In this section, we state our first result of approximation of the Steiner problem, in a slightly more general framework than (2.1). To this aim, we introduce a probability measure µ on R 2 , with compact support, and define K µ := {K ⊂ R 2 ; compact, connected, and s.t. spt (µ) ⊂ K} .

We then investigate the problem

inf{H 1 (K) ; K ∈ K µ }.
(2.4)

Notice that here µ is only important through its support. If the infimum in the above problem is finite, then the problem admits a solution as a direct consequence of Blaschke and Golab's Theorem; notice that in general the minimal set is not unique. In the case of problem (2.1), where we consider a finite set of points {x i } =: D, one can choose any measure µ such that spt (µ) = D, for instance µ = 1 D i δ x i . To approximate (2.4) we introduce an open set Ω containing the convex hull of D. Observing that a minimizer for problem (2.4) will always stay inside Ω, this problem is equivalent to min{H 1 (K) ; K ∈ K µ and K ⊂ Ω} .

(2.5)

Now, to approximate (2.5) we first choose an arbitrary point x 0 ∈ D. Then recalling the definition of D(ϕ; •, •) in (2.2) we introduce the family of functionals defined on L 2 (Ω) by

F ε (ϕ) : = 1 4ε Ω (1 -ϕ) 2 dx + ε Ω |∇ϕ| 2 dx + 1 λ ε Ω D(ϕ; x, x 0 ) dµ(x) + ε q+1 Ω |∇ϕ| q dx (2.6) if ϕ ∈ W 1,q (Ω) satisfies λ 2 ε ≤ ϕ ≤ 1, ϕ = 1 on ∂Ω, and F ε (ϕ) = +∞ otherwise.
Our approximation result is as follows.

Theorem 2.1.1. For all ε > 0 let ϕ ε be a minimizer of F ε , where λ ε → 0.

Consider the sequence of functions d ϕε (x) = D(ϕ ε ; x, x 0 ) and suppose that they uniformly converge to a certain function d. Then the set K := {d = 0} is compact and connected and is a solution to Problem (2.4).

Remark 2.1.1. Notice that the assumption of d ϕε converging to a function d is not restrictive since they are all 1-Lipschitz functions, and thus they always converge uniformly, up to a subsequence.

Elements of proof

Regularity properties of compact connected sets of finite length Proposition 2.1.2. Let K ⊂ R 2 be a compact connected set such that L := H 1 (K) < +∞. Then :

(i) There exists a Lipschitz surjective mapping f : [0, 1] → K. In particular K is arcwise connected.

(ii) For H 1 -a.e. x ∈ K there exists a tangent line T x ⊂ R 2 , in the sense that

lim r→0 d H K ∩ B(x, r), T x ∩ B(x, r) r = 0 ,
where d H stands for the Hausdorff distance.

PHASE-FIELD APPROXIMATION ENFORCING CONNECTEDNESS 13

(iii) For all ε ∈ (0, 1/2) there exists r 0 > 0 such that

π(K ∩ B(x, r)) ⊇ T x ∩ B(x, (1 -ε)r) ∀r ≤ r 0 , (2.7)
where π : R 2 → T x denotes the orthogonal projection onto the line T x identified with R with origin at x.

A lower bound on a new geometric integral quantity

For any closed set A ⊂ R 2 , λ > 0 and direction ν ∈ S 1 , we denote by A λ,ν the λ-enlargement of A in the direction ν defined by

A λ,ν := {x + tν ; |t| ≤ λ and x ∈ A}.
Next, we denote by L 2 the two dimensional Lebesgue measure on R 2 and consider for ν ∈ S 1 the function ν → L 2 (A λ,ν ). It is Borel-measurable and the following quantity is well defined:

I λ (A) := 1 2λπ S 1 L 2 (A λ,ν ) dH 1 (ν) . (2.8)
Lemma 2.1.3. Let (A n ) n∈N be a sequence of compact connected subsets of R 2 converging for the Hausdorff distance to a compact and connected set A.

Then for any 0 ≤ L < H 1 (A) and for every λ > 0 small enough (depending on A), one can find n 0 ∈ N such that

I λ (A n ) ≥ 1 5π L , ∀n ≥ n 0 .

An elementary inequality on the total variation

When I = (a, b) ⊂ R is an interval, we denote by Var(f, I) the total (pointwise) variation of the one-variable function f :

I → R defined by Var(f, I) := sup{ i |f (t i ) -f (t i+1 )|; a = t 0 < t 1 < • • • < t N = b}. If J ⊂ R is open, we define Var(f, J) = I Var(f, I)
, where the sum is taken over all the connected components of J.

Lemma 2.1.4. Let J ⊂ R be an open and bounded set and f : J → R a continuous function on J. Take a finite number of intervals

I ± i ⊂ J, 1 ≤ i ≤ N , satisfying the following properties: i) Ī- i = [a - i , c i ] and Ī+ i = [c i , a + i ] with a - i < c i < a + i , ii) a + i < a - i+1 for all 1 ≤ i ≤ N -1.
Then denoting by m ± i the average of f on I ± i , there holds

Var(f, J) ≥ 1≤i≤N m + i + m - i -2f (c i ).

A useful recovery sequence

Here we recall a standard construction for the Γ-limsup inequality, that is obtained by studying the optimal profile in one dimension for minimizers of an energy of Modica-Mortola type.

Lemma 2.1.5. Let Ω ⊂ R 2 be open and K ⊂ Ω be a compact and connected set. For any r > 0 we will set K r := {x ∈ R 2 ; dist(x, K) ≤ r}. Let k ε → 0 be given. Then for all ε > 0 there exists a function

ϕ ε ∈ W 1,2 loc (R 2 ) ∩ C 0 (R 2 ) equal to k ε on K ε 2 , equal to 1 on R 2 \ K ε 2 +2ε| ln(ε)| , and such that lim sup ε→0 R 2 ε|∇ϕ ε | 2 + (1 -ϕ ε ) 2 4ε ≤ H 1 (K) (2.9)
(the integral being indeed performed on a neighborhood of Ω, since ϕ ε is the constant 1 outside).

Proof. Since K is connected, it automatically satisfies the lower Ahlforsregularity inequality inf x∈K,r<r 0

H 1 (K ∩ B(x, r)) 2r > 0 (2.10)
which is enough to guarantee that the Minkowsky content and the Hausdorff measure coincide (see for e.g. [8, Theorem 2.104, page 110]), namely,

lim r→0 L 2 (K r ) 2r = H 1 (K). (2.11)
Then, for some suitable infinitesimals a ε , b ε that will be fixed below, we take a function ϕ ε similar to the one considered in the proof of Theorem 3.1. in [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF], namely

ϕ ε =      k ε on K bε 1 on Ω \ K bε+aε k ε + λ ε 1 -exp bε-dist(x,K) 2ε on K bε+aε \ K bε where λ ε = 1-kε 1-exp(-aε/(2ε))
, which ensures that ϕ ε is continuous. Fixing b ε = ε 2 and a ε = -2ε ln(ε), so that a ε /ε → +∞ and λ ε → 1, one obtains (2.9). Remark 2.1.2. Notice that from the explicit formulas for ϕ ε , one can also estimate higher summability for |∇ϕ ε |, with bounds depending on ε. Without any intention to be sharp, it is easy to check |∇ϕ ε | ≤ C/ε, and hence

ε q+1 |∇ϕ ε | q ≤ Cε → 0. Semicontinuity of ϕ → D(ϕ; x, y)
The existence of minimizers for the functional F ε defined by (2.6) is based on the following lemma.

Lemma 2.1.6. Let ϕ n be a sequence of continuous functions on a compact domain Ω, such that ϕ n (x) ≥ c > 0 for every x and every n, and ϕ n → ϕ uniformly. Then D(ϕ; x, y) ≤ lim inf n D(ϕ n ; x, y) for every pair (x, y) ∈ Ω.

The main liminf inequality

First, we define

S ε (ϕ, µ, y) := 1 4ε Ω (1 -ϕ) 2 dx + ε Ω |∇ϕ| 2 dx + 1 λ ε Ω D(ϕ; x, y) dµ(x) ,
for any triple (ϕ, µ, y) where ϕ is a continuous function on Ω, µ is a finite positive measure on Ω, and y ∈ Ω. Here the parameter λ ε is considered as a fixed parameter depending on ε.

Lemma 2.1.7. Let Ω ⊂ R 2 be a bounded open set, µ ε be a family of measures on Ω and ϕ ε ∈ H 1 (Ω) ∩ C(Ω) be a family of functions satisfying ϕ ε = 1 on ∂Ω and 0 ≤ ϕ ε ≤ 1 in Ω, λ ε → 0 and x ε a sequence of points of Ω such that

(i) x ε → x 0 for some x 0 ∈ Ω, (ii) µ ε * µ for some measure µ on Ω, (iii) D(ϕ ε ; x, x ε ) converges uniformly to some 1-Lipschitz function d(x) on Ω, (iv) sup ε>0 S ε (ϕ ε , µ ε , x ε ) ≤ C < +∞. Then (a) the compact set K := {x ∈ Ω ; d(x) = 0} is connected, (b) x 0 ∈ K, (c) spt (µ) ⊂ K, (d) H 1 (K) ≤ lim inf ε→0 S ε (ϕ ε , µ ε , x ε ).
Proof. Up to extracting a subsequence in ε, we may assume that the liminf in (d) is a true limit. Our first goal is to realize the set {d = 0} as the Hausdorff converging limit of level sets of D(ϕ; x, x 0 ). For every ε and δ > 0, let us define

K ε,δ := {x ∈ Ω ; D(ϕ ε ; x, x ε ) ≤ δ}.
Since these sets are all compact sets contained in Ω, up to a subsequence we can assume that K ε,δ converges to some K δ for the Hausdorff distance when ε → 0. Next we define

K := δ>0 K δ .
Using that K is connected, one can prove that K = K. The fact that spt (µ) ⊆ K is an easy consequence of Ω D(ϕ ε ; x, x ε ) dµ ε → 0 , which gives, thanks to the uniform convergence D(ϕ ε ; x, x ε ) → d(x) and the weak convergence of the measures, Ω d(x) dµ(x) = 0.

At this stage, (a), (b), (c) are proved and it remains to prove (d), which is achieved in two main steps. The first step consists in finding a weaker version of (d), namely

H 1 (K) < +∞.
(2.12)

This will automatically imply that K is rectifiable because it is also compact and connected. Once we know that K is rectifiable, we will use the existence of a tangent line almost everywhere to deduce the more precise claim (d). This will be the purpose of the second step.

Step 1. Rectifiability of K. We first prove (2.12). Our strategy is to find a certain region where ϕ ε admits a transition between the values 0 and 1, in order to take into account its total variation. At the same time we need this region to be connected and converge to the set K, in order to bound its length by use of a version of Golab theorem that relies rely on the quantity I λ defined in (2.8), instead of the usual H 1 length.

Let us now be more precise. Fix δ 0 , τ 0 > 0, and let {z 1 , z 2 , . . . , z N } ⊆ K be a τ 0 -network in K, which means

K ⊆ 1≤i≤N B(z i , τ 0 ).
Due to the uniform convergence d ϕε (•, x ε ) → d and to the fact that d(z i ) = 0 for 1 ≤ i ≤ N , there exists ε 1 > 0, depending on δ 0 and τ 0 , such that the following holds: for any ε < ε 1 , there exists a C 1 regular curve Γ ε i (of finite length) connecting z i to x ε and satisfying

Γ ε i ϕ ε (s) dH 1 (s) < δ 0 , ∀1 ≤ i ≤ N.
(2.13)

Now we consider

Γ ε := 1≤i≤N Γ ε i .
Our goal is to estimate the quantity I λ (Γ ε ) (defined by (2.8)).

In view of applying Lemma 2.1.3, we denote by Γ 0 the Hausdorff limit of Γ ε (which exists up to extraction). Let λ > 0 be a small enough parameter, and ε 0 > 0 be given by Lemma 2.1.3 applied with L < H 1 (Γ 0 ), in such a way that

L ≤ 5πI λ (Γ ε ) ∀ε < ε 0 . (2.14)
Let now ν ∈ S 1 be an arbitrary direction: for any t ∈ R we denote by L t the affine line Rν + tν ⊥ . Since Γ ε is a finite union of curves of finite length, we know that

H 0 (L t ∩ Γ ε ) < +∞ , for a.e. t ∈ R.
For any such t, we note {x j } j∈J the finite set L t ∩ Γ ε , where J = {1, . . . , j 0 }. We identify x j with its coordinate on the line L t and we assume that they are labelled in increasing order, i.e. x j < x j+1 . Next, we decompose the relative interior Int((Γ ε ) λ,ν ∩ L t ) as follows :

Int((Γ ε ) λ,ν ∩ L t ) = j∈J I - j ∪ I + j ,
where

I + j := [x j , x j + λ) if |x j+1 -x j | ≥ 2λ [x j , x j +x j+1 2
] otherwise

I - j := (x j -λ, x j ] if |x j -x j-1 | ≥ 2λ ( x j +x j-1 2 , x j ] otherwise.
We also define x 0 = x 1 -diam(Ω) and x j 0 +1 = x j 0 + diam(Ω) so that I - 1 and I + j 0 are well defined.

Let P (s) := s -s 2 /2 be the primitive of s → (1 -s) satisfying P (0) = 0 and P (1) = 1/2. Extending ϕ ε by 1 to the whole R 2 and arguing as Modica and Mortola [START_REF] Modica | Il limite nella Γ-convergenza di una famiglia di funzionali ellittici[END_REF], using the inequality 1 4ε a 2 + εb 2 ≥ ab, we infer that

C ≥ S ε (ϕ ε , µ ε , x ε ) ≥ 1 4ε R 2 (1 -ϕ ε ) 2 dx + ε R 2 |∇ϕ ε | 2 dx ≥ R 2 (1 -ϕ ε )|∇ϕ ε | dx ≥ R 2 |∇(P (ϕ ε ))| dx.
Hence, we can go on with

C ≥ R 2 |∇ ν (P (ϕ ε ))| dx = R Var(f t , L t ∩ Ω ) dt, (2.15)
where f t := P (ϕ ε )| Lt . On the other hand, applying Lemma 2.1.4 we can write

Var(f t , L t ∩ Ω) ≥ Var(f t , Int((Γ ε ) λ,ν ∩ L t )) ≥ j∈J m + j + m - j -2f t (x j ),
where m ± j denotes the average of f t on I ± j . Since H 1 (I ± j ) ≤ λ for all j, and f t ≥ 0, we deduce that

Var(f t , L t ∩ Ω) ≥ 1 λ (Γε) λ,ν ∩Lt f t (s) ds -2 j∈J f t (x j ).
Integrating over t ∈ R, applying Fubini's Theorem and using (2.15) it follows that 1 λ (Γε) λ,ν

P (ϕ ε (x)) dx -2 R j∈J f t (x j ) dt ≤ C. (2.16)
Now we estimate the second term in the left hand side of (2.16). The coarea formula (see for instance [8, Section 2.12]) applied on the 1-rectifiable set Γ ε provides

R j∈J f t (x j ) dt = Γε P (ϕ ε (x))C ν (x) dH 1 (x),
where C ν (x) denotes the one dimensional co-area factor associated with the orthogonal projection on ν ⊥ . Since the latter mapping is 1-Lipschitz, it is easy to verify that

|C ν (x)| ≤ 1 yielding R j∈J f t (x j ) dt ≤ 1≤i≤N Γ ε i P (ϕ ε (x)) dH 1 (x) ≤ 1≤i≤N Γ ε i ϕ ε (x) dH 1 (x)
(2.17)

≤ N δ 0 . (2.18)
Here (2.17) is a consequence of P (ϕ ε ) = ϕ ε -ϕ 2 ε /2 ≤ ϕ ε , and (2.18) of (2.13).

Returning now to (2.16), we have proved that

1 λ (Γε) λ,ν P (ϕ ε (x)) dx ≤ C + 2N δ 0 . (2.19)
where C is a uniform constant, but N depends on τ 0 . From (2.19) we get

1 2λ (L 2 ((Γ ε ) λ,ν ) = 1 λ (Γε) λ,ν P (1) dx ≤ C + 2N δ 0 + 1 λ R 2 |P (1) -P (ϕ ε )| dx,
and finally, averaging over ν ∈ S 1 and using (2.14),

L ≤ 10πC + 40πN δ 0 + 20π λ R 2 |P (1) -P (ϕ ε )| dx . (2.20)
Observing that P (ϕ ε ) → P (1) strongly in L 1 (R 2 ), we can let ε → 0 in the previous inequality, and using that L being an arbitrary number smaller than H 1 (Γ 0 ), we obtain

H 1 (Γ 0 ) ≤ 10πC + 40πN δ 0 .
(2.12) can be deduced from the previous inequality by letting δ 0 and τ 0 successively go to zero, and taking the corresponding Hausdorff limit.

Step 2. More precise upper bound on H 1 (K). Now that we know that K is of finite length, and hence rectifiable since it is connected, we will improve the previous upper bound and prove (d). We will use a similar argument as the one used in step 1, with the major difference that instead of being averaged in all directions, the variation of ϕ ε will now be estimated locally along vertical directions orthogonal to the tangent line.

For this purpose we consider the following family of measures supported on Ω:

m ε = 1 4ε (1 -ϕ ε ) 2 + ε|∇ϕ ε | 2 L 2 | Ω .
Thanks to the bound in (iv), we can assume that m ε is weakly- * convergent to a measure m supported on Ω. Applying Lemma 2.1.2 to the set K, we know that H 1 -a.e. point x ∈ K admits a tangent line T x . We assume without loss of generality that x = 0 and T x = Re 1 , and denote by π the orthogonal projection on T x . Let 0 < λ < 1 be fixed (very close to 1). Proposition 2.1.2 guarantees that for some r 0 > 0 and for all r ≤ r 0 ,

π(K ∩ B(x, r)) ⊇ [-λr, λr].
Then we consider the rectangle

C λ (r) := [-λr, λr] × [-hr, hr],
with h = √ 1 -λ 2 , so that C λ (r) ⊂ B(x, r). We want to estimate m ε (C λ (r)), for r small.

Let β > 0 be a very small parameter compared to h. Up to taking a smaller r 0 we may also suppose that for all r ≤ r 0 , K ∩ C λ (r) ⊂ W (β, r), where W (β, r) is a small strip near the tangent of relative width β, namely

W (β, r) := {y ∈ C λ (r); d(y, T x ) ≤ rβ}.
(2.21)

Let us define

δ ε := ||D(ϕ ε ; •, x ε ) -d|| L ∞
, which is a sequence converging to 0 as ε → 0, and consider the sets K ε,δε , which converge in the Hausdorff topology to K. Using the definition of the tangent line, for all t ∈ [-λr, λr] we can find a point z t that belongs to π -1 (t) ∩ K (and hence also to π -1 (t) ∩ K ε,δε ). Also, for every z t there exists a curve Γ ε (t) connecting z t to x ε and such that

Γε(t) ϕ ε (s) ds ≤ δ ε .
Due to the fact that x ε lies outside C λ (r), one can also prove that the curve Γ ε (t) must exit C λ (r). More precisely, denoting by C + and C -the two connected components of ∂C λ (r) ∩ W (β, r), C ± := {y ∈ ∂C λ (r); d(x, (±r, 0)) ≤ βr},

we must have Γ ε (t) ∩ C + = ∅ or Γ ε (t) ∩ C -= ∅.
Let us define

t R := inf{t ∈ [-r, r]; Γ ε (t) ∩ C + = ∅} t L := sup{t ∈ [-r, r]; Γ ε (t) ∩ C -= ∅}.
Notice that we necessarily have t R ≤ t L , otherwise for all the points t ∈ ]t L , t R [, the curve Γ ε (t) would meet neither C + nor C -, which is impossible. Then, take t L < t L and t R < t R such that |t R -t L | < ε and define

Γ ε := Γ ε (t L ) ∪ Γ ε (t R ).
The set Γ ε is not necessarily connected, but we have .22) For every t ∈ I λ,r,ε we denote by g t a point in Γ ε ∩ π -1 (t).

π(Γ ε ) ⊃ I λ,r,ε := [-λr, λr] \ ]t L , t R [. ( 2
Let us now estimate

m ε (C λ (r)) = 1 4ε C λ (r) (1 -ϕ ε ) 2 dx + ε C λ (r) |∇ϕ ε | 2 dx ≥ C λ (r) (1 -ϕ ε )|∇ϕ ε | dx ≥ C λ (r) |∇(P (ϕ ε ))| dx ≥ C λ (
m ε (C λ (r)) ≥ 1 r(h + β) I λ,r,ε ×[-hr,hr] P (ϕ ε (x)) dx -2 I λ,r,ε f t (g t ) dt ≥ 1 r(h + β) C λ (r) P (ϕ ε (x)) dx - εh 2(h + β) -2 I λ,r,ε f t (g t ) dt .
(2.24)

To estimate the last term in the right hand side of (2.24) we use the same argument as for (2.16) relying on the co-area formula, and get

I λ,r,ε f t (g t ) dt ≤ 2δ ε . (2.25)
Coming back to (2.24), we obtain for all r ≤ r 0

m ε (C λ (r)) ≥ 1 r(h + β) C λ (r) P (ϕ ε (x)) dx - ε 2 -4δ ε . (2.26)
Passing to the limsup in ε → 0, using δ ε → 0, together with the facts that ϕ ε → 1 strongly in L 1 , that P (1) = 1/2, that C λ (r) is closed and that m ε converges weakly- * to m we get

m(C λ (r)) ≥ lim sup ε→0 m ε (C(x, r)) ≥ 1 2r(h + β) L 2 (C λ (r)).
Recalling that C λ (r) ⊆ B(x, r) we get m(B(x, r)) ≥ Applying [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]Theorem 2.56], we find that m ≥ H 1 | K and conclude that

H 1 (K) ≤ lim inf ε→0 1 4ε Ω (1 -ϕ ε ) 2 dx + ε Ω |∇ϕ ε | 2 dx .
Proof of Theorem 2.1.1

The existence of a minimizer for F ε is easy, since the term ε q+1 Ω |∇ϕ| q guarantees W 1,q bounds on any minimizing sequence, and hence one can deduce uniform convergence. The lower bound λ 2 ε ≤ ϕ allows us to use Lemma 2.1.6 for the term involving d ϕ , and the other terms of F ε are trivially semicontinuous.

Defining

τ n := sup k≥n |D(ϕ ε k ; x, x 0 ) -d(x)| ∞ , it is easy to see that K n := {x ∈ Ω; D(ϕ εn ; x, x 0 ) ≤ τ n }
converges, for the Hausdorff distance, to the set K = {d(x) = 0}. Now let K 0 ∈ K µ be a minimizer for the Steiner problem (2.5). In particular, K 0 ⊂ Ω because it is contained in the convex hull of the support of µ.

Then, let a ε and b ε be the same parameters as in Lemma 2.1.5, and let ψ ε be the family of function given by Lemma 2.1.5, with k ε = λ 2 ε . For ε small enough we have ψ ε = 1 on ∂Ω. Since the ϕ ε are minimizers of F ε , we have F εn (ϕ εn ) ≤ F εn (ψ εn ). Notice that D(ψ ε ; x, x 0 ) ≤ k ε H 1 (K 0 ) for all x, thus lim ε→0 1 λε Ω D(ψ ε ; x, x 0 ) dµ(x) = 0 and that, thanks to Remark 2.1.2, we also have lim ε→0 ε q+1 Ω |∇ψ ε | q = 0. Hence, we infer that

lim inf n F εn (ϕ εn ) ≤ lim sup n F εn (ψ εn ) ≤ H 1 (K 0 ) . (2.27)
In particular S εn (ϕ εn , µ εn , x 0 ) ≤ F εn (ϕ εn ) ≤ C so that Lemma 2.1.7 applies, thus we obtain that K ∈ K µ and that

H 1 (K) ≤ lim inf n→+∞ S εn (ϕ εn , µ εn , x 0 ) ≤ lim inf n→+∞ F εn (ϕ εn ) .
Gathering with (2.27) we deduce that H 1 (K) ≤ H 1 (K 0 ) which proves that K is a minimizer.

Minimization over phase field/curve pairs (u, γ)

In the aim of approximating numerically the minimization problem min

u∈1+H 1 0 (Ω) F ε (u) , (2.28)
where F ε is the functional defined in (2.6), a natural idea is to compute a time-discretization of the L 2 -gradient flow of F ε (u), i.e,

u t = -∇F ε (u) .
However, this strategy raises some difficulties in the computation of the gradient of the geodesic terms with respect to u. Even if numerical methods have been developed to differentiate geodesics with respect to the metric [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF], in practice, these computations may exhibit some anisotropy issues, or more generally, a dependence on the spatial discretization of the domain. Besides, the cost of the numerical differentiation of the geodesic distance can become prohibitive if the number of Steiner points N increases.

In order to avoid the differentiation of the geodesic distance with respect to the metric, we have developed in [B6] a different approach, consisting in dissociating the minimization problem (2.28) by introducing an extra variable γ := (γ i ) 1≤i≤N , where each γ i is a Lipschitz curve joining a base point a 0 to the point a i in Ω.

Dissociation and regularization of functional F ε

The first step is to regularize the geodesic distance term of functional F ε , defined by (2.6). Fix a positive finite measure µ, compactly supported in Ω, and a base point a 0 ∈ Ω. For a given set of parameters ε, λ ε , δ ε ∈ (0, 1), we consider the functional

F µ ε : H 1 (Ω) ∩ L ∞ (Ω) → [0, +∞) defined by F µ ε (u) := ε Ω |∇u| 2 dx + 1 4ε Ω (u -1) 2 dx + 1 λ ε Ω D δ ε + u 2 ; a 0 , x dµ ,
where, in the D-term, δ ε + u 2 denotes the precise representative of the Sobolev function

δ ε + u 2 ∈ W 1,1 (Ω) ∩ L ∞ (Ω).
In this way, the value of D δ ε + u 2 ; a 0 , x only depends on a 0 , x, and the equivalence class of δ ε + u 2 . Moreover, the function x → D δ ε +u 2 ; a 0 , x turns out to be (δ ε + u 2 L ∞ (Ω) )-Lipschitz continuous, so that F µ ε is well defined. In the definition of the new functional F µ ε , the parameter δ ε can be seen as an elliptic regularisation term. In turn, the use of term u 2 enables us to derive a linear elliptic equation for u (see (2.37)). As will be developed below, a large part of the arguments of the proofs rests on this equation and rather classical linear estimates.

We are interested in the minimization problem min

u∈1+H 1 0 (Ω)∩L ∞ (Ω) F µ ε (u) . (2.29)
Our first main result deals with existence and regularity of solutions.

Theorem 2.2.1. Problem (2.29) admits at least one solution. In addition, any solution u ε belongs to W 1,p (Ω) for every p < ∞ (in particular, u ε ∈ C 0,α (Ω) for every α ∈ (0, 1)), and 0 ≤ u ε ≤ 1.

We now describe the asymptotic behavior of minimizers of F µ ε as ε → 0. For this issue, we shall assume (for simplicity) that the two parameters λ ε and δ ε satisfy the following relation:

λ ε -→ ε→0 0 and δ ε = λ β ε for some β ∈ (1, 2
} k∈N ⊂ 1 + H 1 0 (Ω) be such that F µ ε k (u k ) = min 1+H 1 0 (Ω) F µ ε k for each k ∈ N .
There exist a (not relabeled) subsequence and a compact connected set K * ⊂ Ω such that {u k ≤ t} → K * in the Hausdorff sense for every t ∈ (0, 1). In addition, K * solves the Steiner problem (2.31) relative to {a 0 } ∪ spt µ , and the following holds:

(i) F µ ε k (u k ) → H 1 (K * ); (ii) D δ ε k + u 2 k ; a 0 , x → dist(x, K * ) uniformly on Ω 0 ; (iii) u k → 1 in C 2 loc (Ω \ K * ).
Throughout this section, we assume that the measure µ has finite support, i.e.,

µ = N i=1 β i δ a i (2.32)
for some distinct points a 1 , . . . , a N ∈ Ω and coefficients β i > 0. We fix a a 0 ∈ Ω (possibly equal to one of the a i 's), and to the resulting collection of points, we associate the following space of Lipschitz curves

P(a 0 , µ) := - → γ = (γ i ) N i=1 : γ i ∈ P(a 0 , a i ) ,
where we have set

P(a, b) := γ ∈ Lip([0, 1]; Ω) : γ(0) = a and γ(1) = b .
We endow P(a 0 , µ) with the topology of uniform convergence. For -→ γ ∈ P(a 0 , µ), we write

Γ(γ i ) := γ i ([0, 1]) and Γ( - → γ ) := N i=1 γ i ([0, 1]) .
Dissociation of the geodesic term in F µ ε

In order to get more flexibility in the minimization of F µ ε , we come back to the definition (2.2) of the geodesics and relax problem (2.29) in the space H 1 (Ω) × P(a 0 , µ).

For a given -→ γ ∈ P(a 0 , µ), we consider the functional

E µ ε (•, - → γ ) : H 1 (Ω) → [0, +∞] defined by E µ ε (u, - → γ ) := ε Ω |∇u| 2 dx + 1 4ε Ω (u -1) 2 dx + 1 λ ε N i=1 β i Γ(γ i ) (δ ε + u 2 ) dH 1 , (2.33
) where each term Γ(γ i ) (δ ε + u 2 ) dH 1 is understood as the integration of the precise representative of δ ε + u 2 with respect to the measure H 1 Γ(γ i ).

By the very definition of F µ ε , the functional E µ ε relates to F µ ε through the formula

F µ ε (u) = inf -→ γ ∈P(a 0 ,µ) E µ ε (u, - → γ ) ∀u ∈ H 1 (Ω) ∩ L ∞ (Ω) . (2.34)
Last identity is the key ingredient to investigate existence and regularity of minimizers of F µ ε .

Existence and regularity of minimizing pairs

The study of the regularity of minimizers of E µ ε requires the introduction of a mild constraint, formulated in terms of Ahlfors constant. We denote by Al(K) the Ahlfors constant of a rectifiable set K, defined by

Al(K) := sup H 1 (K ∩ B(x, r)) r : r > 0 , x ∈ K .
Then, for any constant Λ ≥ 2, we set

P Λ (a 0 , µ) := - → γ ∈ P(a 0 , µ) : Al Γ(γ i ) ≤ Λ for each i . ( 2
E µ ε (u, - → γ ε ) . Corollary 2.2.4. Assume that µ is of the form (2.32). The functional F µ ε admits at least one minimizer u ε in 1 + H 1 0 (Ω) ∩ L ∞ (Ω).
In addition, any such minimizer belongs to W 1,p (Ω) for every p < ∞ (in particular, u ε ∈ C 0,α (Ω) for every α ∈ (0, 1)). Moreover, there exists -→ γ ε ∈ P(a 0 , µ)

such that (u ε , - → γ ε ) is a minimizing pair of E µ ε in (1 + H 1 0 (Ω)) × P(a 0 , µ).

Elements of proof

The minimization problem with prescribed curves

In this paragraph, we investigate the minimization problem min

u∈1+H 1 0 (Ω) E µ ε (u, - → γ ) (2.36)
for a prescribed set of curves -→ γ belonging to P Λ (a 0 , µ), defined by (2.35).

The introduction of this constraint is motivated by the following lemma, consequence of a classical result due to Meyers and Ziemer [START_REF] Meyers | Integral inequalities of Poincaré and Wirtinger type for BV functions[END_REF].

Lemma 2.2.5. If - → γ ∈ P Λ (a 0 , µ), then the functional B µ [ - → γ ] : (u, v) ∈ H 1 (Ω) × H 1 (Ω) → N i=1 β i Γ(γ i ) uv dH 1
defines a symmetric, nonnegative, and continuous bilinear form on

H 1 (Ω) satisfying B µ [ - → γ ] ≤ C Ω µ Λ ,
for some constant C Ω depending only on Ω.

Given -→ γ ∈ P Λ (a 0 , µ), we now rewrite for u ∈ H 1 (Ω),

E µ ε (u, - → γ ) = ε Ω |∇u| 2 dx + 1 4ε Ω (u -1) 2 dx + 1 λ ε B µ [ - → γ ](u, u) + δ ε λ ε N i=1 β i H 1 (Γ(γ i ))
.

By the previous lemma,

E µ ε (u, - → γ ) < ∞ for every u ∈ H 1 (Ω), and E µ ε (•, - → γ ) is
lower semicontinuous with respect to weak convergence in H 1 (Ω). Owing to the strict convexity of the functional E µ ε (•, -→ γ ), we conclude to the following Theorem 2.2.6. Given -→ γ ∈ P Λ (a 0 , µ), problem (2.36) admits a unique solution u-→ γ .

For -→ γ ∈ P Λ (a 0 , µ), we shall refer to u-→ γ as the potential of -→ γ . It satisfies the Euler-Lagrange equation

     -ε 2 ∆u-→ γ = 1 4 (1 -u-→ γ ) - ε λ ε B µ [ - → γ ](u-→ γ , •) in H -1 (Ω) , u-→ γ = 1 on ∂Ω . (2.37)
The next lemmas gather some regularity estimates on u-→ γ with explicit dependence on the parameters. Their proof is based on the comparison of function v = 1 -u-→ γ with explicit supersolutions of the equation -4ε 2 ∆v + v = 0 with homogeneous Dirichlet boundary conditions, in a ball centered at x 0 and contained in Ω, in the spirit of [START_REF] Bethuel | Asymptotics for the minimization of a ginzburg-landau functional[END_REF].

Lemma 2.2.7. Let - → γ ∈ P Λ (a 0 , µ). The potential u-→ γ satisfies 0 ≤ u-→ γ ≤ 1 a.e. in Ω, and u-→ γ ∈ C ∞ Ω \ Γ( - → γ ) . Lemma 2.2.8. Let - → γ ∈ P Λ (a 0 , µ). At every x 0 ∈ Ω \ Γ( - → γ ) satisfying dist(x 0 , Γ( - → γ )) ≥ ε, we have 0 ≤ 1 -u-→ γ (x 0 ) ≤ exp - dist(x 0 , Γ( - → γ )) 3ε . Proposition 2.2.9. If - → γ ∈ P Λ (a 0 , µ), then u-→ γ ∈ C 0,α (Ω) for every 0 < α < 1, and
u-→ γ C 0,α (Ω) ≤ C α (1 + Λ µ λ -1 ε ) ε α ,
for some constant C α depending only on α, Ω, and on the distance between points {a i } and ∂Ω.

Proof of Theorem 2.2.3

A fundamental step in the proof is a replacement procedure revealing the upper Alhfors regularity of optimal curves. For the next statements, we thus introduce the upper Alhfors threshold

Λ ε := 2 + 3 δ ε . (2.38) Lemma 2.2.10. Let u ∈ 1 + H 1 0 (Ω) ∩ L ∞ (Ω) be such that u L ∞ (Ω) ≤ 1, and let - → γ = (γ 1 , . . . , γ N ) ∈ P(a 0 , µ). If for some i 0 ∈ {1, . . . , N }, x ∈ Γ(γ i 0 ),
and r > 0,

H 1 Γ(γ i 0 ) ∩ B(x, r) ≥ Λ ε r , (2.39)
where Λ ε is defined in (2.38), then there exists -→ γ ∈ P(a 0 , µ) of the form

- → γ = (γ 1 , . . . , γ i 0 -1 , γ i 0 , γ i 0 +1 , . . . , γ N )
such that

E µ ε (u, - → γ ) ≤ E µ ε (u, - → γ ) - β i 0 r λ ε .
Next lemma provides the existence of a minimizer -→ γ in P Λε (a 0 , µ) associated to some fixed smooth function u.

Lemma 2.2.11.

Let u ∈ 1 + H 1 0 (Ω) ∩ C 1 (Ω) be such that 0 ≤ u ≤ 1. There exists - → γ = (γ 1 , . . . , γ N ) ∈ P Λε (a 0 , µ) satisfying E µ ε (u, - → γ ) ≤ E µ ε (u, - → γ ) ∀ - → γ ∈ P(a 0 , µ) , (2.40)
and such that each

γ i : [0, 1] → Ω 0 is injective if a i = a 0 , and constant if a i = a 0 .
Next result will allow us to replace an arbitrary pair (u, -→ γ ) by a regular one, with controlled energy.

Lemma 2.2.12. For every σ > 0, u ∈ 1 + H 1 0 (Ω), and -→ γ ∈ P(a 0 , µ), there

exist u σ ∈ 1 + H 1 0 (Ω) ∩ C 1 (Ω) and - → γ σ ∈ P Λε (a 0 , µ) such that 0 ≤ u σ ≤ 1 and E µ ε (u σ , - → γ σ ) ≤ E µ ε (u, - → γ ) + σ .
Now we can move on to the proof of Theorem 2.2.3.

Proof of Theorem 2.2.3

Step 1 (existence). Let {(u n , -→ γ n )} n∈N be a minimizing sequence for E ε over

(1 + H 1 0 (Ω)) × P(a 0 , µ), i.e. lim n→∞ E µ ε (u n , - → γ n ) = inf (1+H 1 0 (Ω))×P(a 0 ,µ) E µ ε .
By Lemma 2.2.12, there is no loss of generality assuming that (

u n , - → γ n ) ∈ C 1 (Ω) × P Λε (a 0 , µ) and 0 ≤ u n ≤ 1.
In addition, by Lemma 2.2.11 we can even assume that, setting -→ γ n = (γ n 1 , . . . , γ n N ), all γ n i 's are injective curves for a i = a 0 , and constant for a i = a 0 . Then we consider the sequence

{(u-→ γ n , - → γ n )} n∈N , where u-→ γ n is the potential of - → γ n , i.e.
, the minimizer of

E µ ε (•, - → γ n ) over 1 + H 1 0 (Ω). Obviously, {(u-→ γ n , - → γ n )} n∈N is still a minimizing sequence by minimality of u-→ γ n . By Proposition 2.2.9, u-→ γ n C 0,α (Ω) ≤ C α (ε) ∀α ∈ (0, 1) ,
for some constant C α (ε) independent of n. By Arzelà-Ascoli Theorem, we can extract a subsequence such that u-→ γ n → u ε uniformly in Ω and weakly in H 1 (Ω) for some function

u ε ∈ 1 + H 1 0 (Ω) ∩ C 0,α (Ω)
, for every α ∈ (0, 1). On the other hand, the energy being invariant under reparametrization, we can assume that each γ n i is a constant speed parametrization of its image Γ(γ n i ). In particular, each

γ n i is a H 1 (Γ(γ n i ))-Lipschitz curve. Since H 1 (Γ(γ n i )) ≤ λ ε δ ε E µ ε (u-→ γ n , - → γ n ) ≤ C(ε) ,
we infer that each sequence {γ n i } n∈N is equi-Lipschitz. Therefore, we can extract a further subsequence such that, for each i ∈ {1, . . . , N }, γ n i → γ ε i uniformly on [0, 1] and weakly* in W 1,∞ (0, 1) for some γ ε i ∈ P(a 0 , a i ). Then we set -→ γ ε := (γ ε 1 , . . . , γ ε N ) ∈ P(a 0 , µ). By the uniform convergence of u-→ γ n towards u ε , the injectivity of each γ n i (for a i = a 0 ) and the weak* convergence

γ n i * γ ε i in W 1,∞ (0, 1), we deduce that lim inf n→∞ Γ(γ n i ) (δ ε + u 2 -→ γ n ) dH 1 ≥ Γ(γ ε i ) (δ ε + u 2 ε ) dH 1 ∀i ∈ {1, . . . , N } .
Since the diffuse part of the energy is clearly lower semicontinuous with respect to weak H 1 -convergence, we conclude that

E µ ε (u ε , - → γ ε ) ≤ lim n→∞ E µ ε (u n , - → γ n ) , and thus (u ε , - → γ ε ) is a minimizer of E µ ε .
Step 2 (regularity). Now we consider an arbitrary minimizer (u

ε , - → γ ε ) of E µ ε in (1+H 1 0 (Ω))×P(a 0 , µ).
Arguing as in the proof of Lemma 2.2.7, we obtain 0

≤ u ε ≤ 1 by minimality of u ε for E µ ε (•, - → γ ε ). In turn, the minimality of - → γ ε for E µ ε (u ε , •) implies that - → γ ε ∈ P Λε (a 0 , µ) by Lemma 2.2.10. Now Theorem 2.2.6 shows that u ε is the potential of - → γ ε . Proof of Theorem 2.2.2
The proof of this theorem is based upon Lemmas 2.1.5 and 2.1.7 from [B11].

Besides, the term u 2 is the key new ingredient to derive a linear elliptic equation for u ε , and obtain rather classical estimates in this framework. Finally, relation (2.30) between λ ε and δ ε is imposed for the following reasons: on one hand the condition δ ε = o(λ ε ) is necessary to derive the Steiner problem in the limit; on the other hand the condition λ 2 ε = o(δ ε ) allows us to apply [B11] in a straightforward way, even if it is probably unnecessary.

Numerical approximation of the Steiner problem

Following the approach described in Section 2.2, we start from the minimization problem (2.29). By (2.34), this problem can be formulated as

inf E ε (u, γ), u ∈ 1 + H 1 0 (Ω), γ ∈ P(a 0 , µ) ,
which suggests to compute a minimizing sequence for F µ ε by considering the following L 2 -gradient flow of E ε :

u t = -∂ u E ε (u, γ) , γ = Argmin γ {E ε (u, γ)} .
Then, an approximation (u n , γ n ) of the flow (u, γ) at time t = nδ t can be obtained by the following time-decoupled scheme:

Step 1 -Computation of γ n as

γ n = Argmin γ∈P(a 0 ,µ) E ε (u n , γ).
Step 2 -Computation of the function u n+1 defined by u n+1 = v(δ t ) where v is the solution to the following PDE:

v t = -∂ u E ε (v, γ n ) , v(•, 0) = u n (•) .
One major advantage of this approach is that v can be computed without requiring the differentiation of the geodesic distance with respect to the metric, which is known to be a very expensive procedure.

However, two major difficulties remain in the computation of step 2:

• the lack of regularity of the solution u, which is merely Hölder-continuous;

• the differentiation of the geodesic term R ε , defined by

R ε (u) = N i=1 Γ(γ i ) (δ ε + u 2 ) dH 1 , (2.41)
whose differential with respect to u is in general a measure.

In this section, we explain how we can slightly modify the previous phase field model to improve the regularity of the solution u, and facilitate the derivation of a simple unconditionally-stable scheme.

Modified phase field models

To enhance the regularity of u, the idea consists in considering a secondorder Cahn-Hilliard type energy to approximate H 1 (K), of the form

P ε (u) = Ω ε 3 |∆u| 2 + 1 ε V (u) dx ,
where V is the single-well potential defined by V (s) = 1 4 (1 -s) 2 . These higher order energies are convenient to regularize the profiles of the minimizers near the singular sets.

Another remarkable property of the second-order approximations of the length is that they can be used in dimension 3 as well. In that case, for instance, one can consider the modified energy

P ε (u) = Ω ε 2 |∆u| 2 + 1 ε 2 V (u) dx ,
where the usual gradient term is replaced by a Laplacian which, as opposed to the gradient, is able to detect 1-dimensional objects in R 3 .

Recall that the numerical computation of the solution v to the equation

v t = -∂ u E ε (v, γ) = 2ε∆v - 1 ε V (v) - 2 λ ε N i=1 H 1 | Γ(γ i ) v,
raises some numerical stability issues, due to the presence of singular measures

H 1 | Γ(γ i ) .
A natural way of regularizing these measures is to consider the convolution with a kernel ρ ε α of size ε α , where α > 0 is a parameter, i.e. ρ ε α = 1 ε 2α ρ(•/ε α ). This leads to the regularized geodesic term

R ε (u, γ) = 1 λ ε Ω N i=1 (ρ ε α * H 1 | Γ(γ i ) ) (δ ε + u 2 ) dx . (2.42)
In the rest of the section, we will denote by ω ε [γ] the function

ω ε [γ] = N i=1 (ρ ε α * H 1 | Γ(γ i ) ). (2.43)
Moreover, in practice and in our numerical experiments, we can observe (see for instance Figure 2.3) that the width of the sublevel sets of u can vary along the support of γ. This results from the use of a penalization term ω ε [γ], whose weight is not uniform along the support of γ. It makes the method inefficient for a large number of attachment points, since the geodesics may not be precisely detected by the phase field u.

Hence, we propose to reduce these adverse effects by considering the weight

ω ε max [γ] = max N i=1 ρ ε α * H 1 | Γ(γ i ) ,
and the associated geodesic term

R ε,max (u, γ) = 1 λ ε Ω ω ε max [γ](δ ε + u 2 ) dx .
(2.44)

Numerical experiments

The motivation of this section is to explain how to compute the solution of the PDE

u t = -∂ u E ε (u, γ) γ = Argmin γ {E ε (u, γ)}.
where E ε (u, γ) is equal to

• Case 1 : Cahn-Hilliard of order 1 and classical smooth geodesic term

E ε (u, γ) = P ε (u) + Rε (u, γ),
• Case 2 : Cahn-Hilliard of order 2 and classical smooth geodesic term

E ε (u, γ) = Pε (u) + Rε (u, γ),
• Case 3 : Cahn-Hilliard of order 2 and concentrated smooth geodesic term E ε (u, γ) = Pε (u) + Rε,max (u, γ).

In each situation, we consider the solution at all time t ∈ [0, T ] in a computation box Q with periodic boundary conditions, associated with the initial condition u(x, 0) = 1.

As explained previously, our numerical approach is based on a splitting method, which constructs the approximative sequence (u n , γ n ) of (u, γ) at time nδ t recursively, as follows:

Step 1 : Computation of γ n as

γ n = min γ∈P(a 0 ,µ) {E ε (u n , γ)} .
Step 2 : Computation of function u n+1 defined by u n+1 = v(δ t ) where v is solution of the following PDE:

v t = 2ε∆v -1 2ε (v -1) -2 λε ω ε [γ n ]v v(•, 0) = u n (•) Recall that ω ε [γ n ] is defined by (2.43).

Numerical scheme

In this paragraph, we give more details on the computation of the geodesic γ n and how to obtain an approximation of u n+1 = v(δ t ).

Computation of γ n Recall that, in cases 1 and 2, the geodesics γ n = γ n i are defined by

γ n = min γ∈P(a 0 ,µ) {E ε (u n , γ)} = min γ∈P(a 0 ,µ) Rε (u n , γ) ,
where Rε is defined by (2.42). Hence, we can compute the geodesics γ n i in two steps.

• We first compute by a Fast Marching Method [93] the weighted distance function x → d x 0 ,ω (x), corresponding to the distance function from a point x 0 to x, associated to the weight w = (δ ε + u 2 ) * ρ ε α .

• The second step consists in computing each geodesic

γ i : [0, 1] → Ω satisfying γ i (0) = x 0 , γ i (1) = x i , associated with the distance d x 0 ,ω .

Computation of u n+1

We adopt a Fourier semi-implicit approach [START_REF] Chen | Applications of semi-implicit Fourier-spectral method to phase field equations[END_REF] to compute u n+1 , by setting

u n+1 -u n δ t = 2ε∆u n+1 - 1 2ε (u n+1 -1)-αu n+1 - 2 λ ε ω ε [γ n ] -α u n . (2.45)
Here α can be viewed as a stabilization parameter for the scheme. Indeed, it is well known [START_REF] Shen | Second-order convex splitting schemes for gradient flows with Ehrlich-Schwoebel type energy: Application to thin film epitaxy[END_REF] that if the explicit part is the L 2 -gradient flow of a concave functional, then the global numerical scheme is stable without any condition on the time step δ t .

The advantage of this semi-implicit approach is that the operator associated to the linear equation is now homogeneous, and remains identical during the iterations. Moreover, using additional periodic boundary conditions, Equation (2.45) can be solved in the Fourier space, which ensures a high accuracy of the spatial approximation of u n+1 .

Spatial discretization and Fourier space

We recall that the P Fourier approximation of a 2D function u in a box

Q = [0, L 1 ] × [0, L 2 ] is given by u P (x) = P/2 p 1 ,p 2 =-P/2+1 c p e 2iπξp•x
where p = (p 1 , p 2 ) and ξ p = (p 1 /L 1 , p 2 /L 2 ). Here c p represents the P 2 first discrete Fourier coefficients of u.

First experiments and comparison of the two different Cahn-Hilliard functionals

The first numerical experiments concern simple configurations of 3 and 4 points. The results are plotted in Figure 2.1 and 2.2. As expected, the solution associated to the operator -ε 3 ∆ 2 is smoother than the one obtained with the classical Cahn-Hilliard functional.

Comparison of the geodesic terms R ε and R ε,max

Since the weight ω ε [γ] at point x depends on the number of geodesics that cross this point, expression (2.42) shows that R ε can be sensitive to the number of endpoints a i , which raises some numerical difficulties to fix the value of γ ε in practice. To reduce these drawbacks, we propose to use the new weight ω ε max [γ], replacing the sum of the contributions of all the geodesics by the maximum contribution, and the associated penalized term R ε,max defined by (2.44).

We can clearly observe (see Figure 2.3) that when R ε is used, the width of the diffuse interface depends on the number of points x i , which is not the case with the new model R ε,max .

Numerical experiments with a large number of points

In view of the previous numerical results, we consider the phase field model

E ε (u, γ) = P ε (u) + R ε,max (u, γ).
We present two numerical experiments in Figure 2.4, with respectively 50 and 100 points randomly distributed in the computation box Q.

We observe a good behavior of the solution which becomes more precise as ε goes to 0. However, we can not guarantee that the iterations converge to a global minimizer of the Steiner problem; in particular, we can observe the presence of a cycle on the final solution obtained using 100 points, which means that this solution is clearly not optimal. This comes from the fact that, at ε fixed, the minimization procedure that we implement relies on the alternate minimization of E ε over γ and u, and thus, may converge to a local minimizer of this functional.

Numerical experiments in dimension 3

We finally propose to test our numerical method in the 3-dimensional case. More precisely, we consider the phase field energy where P ε now satisfies

E µ ε (u, γ) = P ε (u) + R ε,max (u, γ).
P ε (u) = Q ε 2 2 (∆u) 2 + 1 ε 2 V (u) dx.
We plot in figure 2.5 the results of two first experiments that consider 3 and 10 points respectively. Each red ball represents a point a i and the green surface corresponds to the level set (1 + min(u))/2 of u, defined as the boundary of the set {x ∈ Q; u(x) ≤ (1 + min(u))/2}. The approximations of u are computed with the following set of parameters: P = 2 7 , δ t = ε 2 and λ ε = ε 2 . More precisely, we compute u n until we reach a stationary solution, and then divide ε by 2 and compute a new stationary solution.

We conclude this section by a last example where the points a i are located on the corners of a cube. This experiment is very sensitive since for reasons of symmetry, there exist many global minimizers, and consequently, a lot of local minimizers as well. Using the following set of parameter P = 2 6 , δ t = ε 2 and λ ε = 2 , we represent in figure 2.6 the evolution of u during the iterations and until convergence to a stationary solution. The last picture concerns the stationary solution obtained with ε = 2.5/P . We observe that the algorithm seems to converge to an approximation of a global minimizer, which highlights the great potential of the method also in dimension 3.

Conclusion and perspectives

In this chapter, we have presented a variational approximation of Modica-Mortola type for the Steiner tree problem, seen as a model problem for length minimization problems involving a connectivity constraint on onedimensional competitors. We have obtained rigorous convergence results while developing a stable numerical scheme to make the method efficient and reliable in practice.

Hopefully, the results obtained on the Steiner problem should enable us to handle more complex models, where one-dimensional objects are coupled with PDE's. In that sense, we have started to work with R. Ducasse, A. Lemenant and A Zilio on a reaction-diffusion model based on paper by Berestycki, Roquejoffre and Rossi [START_REF] Berestycki | The influence of a line with fast diffusion on Fisher-KPP propagation[END_REF]. In this model, a volumic diffusion is coupled with a diffusion on a segment K. The motivation is to understand the impact of the "road" K on the propagation of a species. Our objective is to generalize the results in [START_REF] Berestycki | The influence of a line with fast diffusion on Fisher-KPP propagation[END_REF] to the case where K is a network, and possibly, to investigate the existence and uniqueness of optimal networks for certain relevant criteria, and approximate them numerically.

Introduction

With the recent development of micro-and nanofluidics, drag reduction for low Reynolds number flows, notably at solid walls, has become a stimulating issue. Therefore, the interaction between a fluid and a solid boundary has been investigated thoroughly, both at the experimental and theoretical levels. At the nanoscale, flows present some striking features, such as the capacity for the material to slip much more easily that one would expect [START_REF] Secchi | Massive radius-dependent flow slippage in carbon nanotubes[END_REF][START_REF] Tunuguntla | Enhanced water permeability and tunable ion selectivity in subnanometer carbon nanotube porins[END_REF]. The origin of this slip is the subject of current debate in the physics community [START_REF] Faucher | Critical knowledge gaps in mass transport through single-digit nanopores: A review and perspective[END_REF][START_REF] Lauga | Microfluidics: The No-Slip Boundary Condition[END_REF].

From a mathematical perspective, one successful strategy initiated in the late 1990's [START_REF] Achdou | Effective boundary conditions for laminar flows over periodic rough boundaries[END_REF][START_REF] Casado-Díaz | Why viscous fluids adhere to rugose walls: A mathematical explanation[END_REF][START_REF] Jäger | On the roughness-induced effective boundary conditions for an incompressible viscous flow[END_REF] to explain the occurrence of slip on solid walls, consists in modeling micro-asperities on the surface and analyzing their effect on the flow by an homogenization process, imposing only a mild non penetration boundary condition on the rugous wall. This so-called "rugosity effect" has been studied quite extensively in the last decades, which has led to a rather complete description of the asymptotic effect of rough patterns on viscous flow [START_REF] Bucur | Boundary behavior of viscous fluids: influence of wall roughness and friction-driven boundary conditions[END_REF][START_REF] Bucur | On the asymptotic limit of the Navier-Stokes system on domains with rough boundaries[END_REF][START_REF] Dalibard | Effective boundary condition at a rough surface starting from a slip condition[END_REF][START_REF] Gérard-Varet | Relevance of the slip condition for fluid flows near an irregular boundary[END_REF][START_REF] Jäger | Couette flows over a rough boundary and drag reduction[END_REF]. The models described in Chapters 4 and 5 fall into this category. In all these examples, the authors introduce sequences of oscillating boundaries Γ ε associated with a small parameter ε controlling the amplitude/period of the roughness, and derive an effective flow by homogenizing the system when ε tends to zero.

In this chapter, we focus on two different models of fluid/solid interaction at the micro-to nanoscale, that consider a fixed, flat boundary Γ instead of rugous walls Γ ε , but resort to:

• either an alternance of no-slip and perfect slip boundary conditions on Γ, modeling rough hydrophobic surfaces (Section 3.2),

• or a reduction of the viscosity coefficient in a small vicinity of the wall (Section 3.3), proposed in [START_REF] Myers | Why are slip lengths so large in carbon nanotubes? Microfluid[END_REF] to model a viscous flow in a nanotube.

Let us point that, although one could expect the averaging procedure usually applied to obtain macroscopic flow equations to fail at the nanoscale, Stokes equations remain surprisingly efficient in nanofluidics [START_REF] Gravelle | Large permeabilities of hourglass nanopores: From hydrodynamics to single file transport[END_REF].

To clarify the presentation, we consider the same physical situation for both models: a three-dimensional Stokes flow between two parallel, infinite plates Γ t (top) and Γ b (bottom). Since we are interested in the fluid solid interaction, we impose periodic lateral boundary conditions in the fluid domain Ω and set

Ω = T 2 × (0, 1) , Γ t = T 2 × {1} , Γ b = T 2 × {0} .
The flow is driven by a given source term f ∈ L 2 (Ω, R 2 ). The solid wall/fluid interaction will be modeled on the lower boundary Γ b , while we impose noslip boundary condition u = 0 on Γ t .

Homogenization of a model of rough hydrophobic surfaces

In the recent years, promising results have been obtained concerning a class of rough hydrophobic surfaces, see for instance [START_REF] Vinogradova | Surface roughness and hydrodynamic boundary conditions[END_REF]. Indeed, by the combination of the chemical and geometrical properties of these surfaces, the hollows of the roughness get filled with gas. Hence, the viscous fluid above does not penetrate: it slips above the hollows, and only sticks at the bumps, reaching the so-called Cassie state.

We focus in this section on a model proposed by Ybert et al. in [START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: Scaling laws for generic geometries[END_REF], where the rough boundary is replaced by a flat plane, divided in small periodic cells (say of side ε 1). Each cell is divided in two zones:

• a no-slip zone, corresponding to a plane projection of the sticky part of the roughness (bumps);

• a slip-zone, corresponding to a plane projection of the slippery part.

Using homogenization techniques, we derive an effective boundary condition as ε goes to zero, depending on the characteristic scale a ε of the no-slip zones. We provide in this way scaling laws for the slip coefficients, for various configurations (patches, riblets). Such laws are in global agreement with the formal computations led in [START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: Scaling laws for generic geometries[END_REF]. All our theoretical results are grounded by numerical computations.

Modelisation of the rough hydrophobic surface

We consider Stokes equation

-∆u + ∇p = f in Ω, div u = 0 in Ω, (3.1)
completed with the no-slip condition

u = 0 on Γ t . (3.2)
On the lower surface Γ b , we assume that u satisfies alternately "perfect slip" and "no slip" boundary conditions, corresponding respectively to the hollows and bumps of the rough hydrophobic surface. More precisely, let ε > 0 and

S ε := [0, ε) 2 ∼ (R/ (εZ)) 2 ,
the elementary square of side ε . For simplicity, we shall assume all along that ε -1 is an integer. Let T ε be a Lipschitz subdomain of S ε , modeling an elementary no-slip zone. Details about T ε will be given right below. From this elementary no-slip zone, we define a global one inside [0, 1) 2 ∼ T 2 :

T ε := k∈[|0,...,ε -1 |] 2 (εk + T ε ) .
Finally, the boundary condition on Γ b is

u 3 = 0 on Γ b , ∂ 3 u = 0 on (T ε ) c × {0}, u = 0 on T ε × {0} . (3.3) It is easily proved that (3.1)-(3.2)-(3.3) has a unique solution (u ε , p ε ) ∈ H 1 (Ω, R 3 ) × L 2 (Ω)/R.
Definition of T ε We will distinguish between two types of no-slip pattern T ε :

• Patches: we assume that

T ε := ( ε 2 , ε 2 ) + a ε T , (3.4)
where ( ε 2 , ε 2 ) is the center of the square S ε , and where the domain T is relatively compact in the square (-1/2, 1/2) 2 , and contains a disk of radius α > 0, centered in the origin (see Figure 3.1). The parameter a ε is a positive number such that a ε < ε. In this case, the no-slip zone is a union of periodically distributed patches.

• Riblets: we assume that

T ε := (εT) × ε 2 + a ε I (3.5)
where

I ⊂ (-1 2 , 1 2
) is an open interval (see Figure 3.2). In this case, the no-slip zone is a union of stripes, invariant in the x 1 -direction. We shall focus on two particular cases:

-f =
= (k 1 , k 2 ) ∈ [|0, . . . , ε -1 |] 2 , the intersection of the no-slip zone T ε with the cell [εk 1 , ε(k 1 +1))×[εk 2 , ε(k 2 +1)) is defined by εk + T ε = εk + a ε T .

Main results

Theorem 3.2.1. (Asymptotic behavior for patches)

Consider the patches configuration (3.4). Let u ε ∈ H 1 (Ω, R 3 ) be the solution of (3.1), (3.2), (3.3). One must distinguish between three cases:

1. Subcritical case: if a ε ε 2 , then u ε ū0 in H 1 (Ω, R 3 ); 2. Supercritical case: if a ε ε 2 , then u ε ū∞ in H 1 (Ω, R 3 );
3. Critical case: there exists a symmetric, positive definite matrix

M 0 such that if a ε /ε 2 → C 0 > 0, then u ε ūC 0 M 0 in H 1 (Ω, R 3 ).
A similar result holds for riblets, that we merely state in the critical case:

Theorem 3.2.2. (Asymptotic behaviour for riblets)

Consider the riblet configuration (3.5). Suppose that lim ε→0 -ε ln(a ε ) = C 0 > 0, and furthermore that f does not depend on x 1 .

Then, u ε ūM rib , where

M rib = π C 0 0 0 2π C 0 . (3.8)
Additionally, when f = e 1 or f = e 2 , the limit system can be simplified:

• if f = e 1 (riblets parallel to the main flow), then ūM rib 2 = 0 and ūM rib

1 satisfies ∂ 3 ūM rib 1 = π C 0 ūM rib 1 on Γ b .
Hence, the slip length is C 0 /π; Hence, the slip length is C 0 /(2π).

• if f = e 2 (

Elements of proof

By classical arguments, the sequence (u ε , p ε ) is bounded in H 1 (Ω, R 3 ) × L 2 (Ω)/R, and consequently, there exists a couple (ū, p) in the same space such that

u ε ū weakly in H 1 (Ω, R 3 ), p ε p weakly in L 2 (Ω)/R.
Using the weak formulation of Eqs. (3.1) and the continuity of the trace operator, one obtains easily that the weak-limit (ū, p) satisfies Eqs. (3.1) and boundary condition (3.2). Similarly, we get ū3 = 0 on Γ b .

Patches configuration

Supercritical case: a ε ε 2 The proof in the supercritical case relies on a quantitative Poincaré inequality, of the form

T 2 ×{0} |u ε | 2 ≤ η(ε) Ω |∇u ε | 2 (3.9)
where the positive function η(ε) goes to zero with ε. Since u ε is bounded in H 1 (Ω, R 3 ), we immediately infer that in the sense of trace, ū vanishes on Γ b , hence ū = ū0 .

Critical and subcritical cases: a ε ε 2 We follow here the strategy of articles [START_REF] Allaire | Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes. I. Abstract framework, a volume distribution of holes[END_REF][START_REF] Allaire | Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes. II. Noncritical sizes of the holes for a volume distribution and a surface distribution of holes[END_REF] by Allaire. These articles deal with the homogenization of the Stokes equations across a network of balls, with a Dirichlet condition at the surface of the obstacles. Notably, in section 4 of [START_REF] Allaire | Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes. I. Abstract framework, a volume distribution of holes[END_REF], the balls are assumed to be distributed along a hypersurface. In the setting considered here, the rough idea is to extend the Stokes solution to the lower half-space by an appropriate symmetry, reducing the problem to the homogenization of Stokes equations across a planar network of patches. Hence, the ideas of [START_REF] Allaire | Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes. I. Abstract framework, a volume distribution of holes[END_REF], based on the construction of correctors and the method of oscillating test functions, essentially apply.

Riblet configurations

In the case of riblets, we recall that T ε is invariant by translation in x 1 . Since f = (f 1 , f 2 , f 3 ) is also independent on the x 1 variable, the solution (u ε , p ε ) of system (3.1)-(3.2)-(3.3) depends only on (x 2 , x 3 ). As a result, the first component of u ε satisfies:

-∆ 2,3 u ε,1 = f 1 in T × (0, 1), u ε,1 = 0 on T × {1}, ∂ 3 u ε,1 = 0 on (T × {0}) \ (ΠT ε × {0}), u ε,1 = 0 on ΠT ε × {0}, (3.10)
where ∇ 2,3 and ∆ 2,3 stand for the gradient (resp. the Laplacian) with respect to the (x 2 , x 3 ) variables and Π is the projection operator defined by Π(x 1 , x 2 ) = x 2 . In the same fashion, (u ε,2 , u ε,3 , p ε ) satisfy the following Stokes problem:

-∆ 2,3 u ε,2 u ε,3 + ∇ 2,3 p ε = f 2 f 3 in T × (0, 1), ∇ 2,3 • u ε,2 u ε,3 = 0 in T × (0, 1), u ε,2 = u ε,3 = 0 on T × {1}, u ε,3 = 0 on T × {0}, ∂ 3 u ε,2 = 0 on (T × {0}) \ (ΠT ε × {0}), u ε,2 = 0 on ΠT ε × {0}. (3.11)
Hence, the original 3d problem reduces to the study of two independent systems, set in the 2d domain T×(0, 1). This change from a 3d to a 2d setting explains the change of scalings between Theorem 3.2.1 and Theorem 3.2.2.

• To handle the Stokes equations (3.11), we proceed similarly as in the case of patches, pondering this time on the 2d results of [START_REF] Allaire | Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes. II. Noncritical sizes of the holes for a volume distribution and a surface distribution of holes[END_REF] about periodic network of disks along a line. We obtain the following limit system:

-∆ 2,3 ū2 ū3 + ∇ 2,3 p = f 2 f 3 in T × (0, 1), ∇ 2,3 • ū2 ū3 = 0 in T × (0, 1), ū2 = ū3 = 0 on T × {1}, ū3 = 0 on T × {0}, ∂ 3 ū2 = 2π C 0 ū2 on T × {0}, (3.12)
where we recall that C 0 := lim ε→0 -ε ln |a ε |.

• As regards the Laplace equation (3.10), the situation is simpler and has been analyzed for a longer time. One may start from the work of Cioranescu and Murat [START_REF] Cioranescu | Un terme étrange venu d'ailleurs[END_REF], instead of [4, section 4]. In our setting, the limit system reads -∆ 2,3 ū1 = f 1 in T × (0, 1), ū1 = 0 on T × {1},

∂ 3 ū1 = π C 0 ū1 on T × {0}. (3.13)
We deduce from systems (3.13) and (3.12) that ū = ūM rib , M rib being given by (3.8). The sub-cases where f = e 1 or f = e 2 follow easily.

Numerical computations

This section is devoted to simulations of system (3.1)-(3.2)-(3.3). For simplicity, we consider a constant source term, say

f = 2e, e ∈ span(e 1 , e 2 ).
The idea is to recover numerically the scalings for the slip length given in Section 3.2.2. However, to observe significant slip implies to consider very small scales: patches of size less than ε 2 , in a grid of side ε. This forbids direct computations. To overcome this difficulty, we rely on a boundary layer approximation of the Stokes flow. Such approximation, often implicitly used in physics papers, has been fully justified in the context of wall laws [START_REF] Jäger | On the roughness-induced effective boundary conditions for an incompressible viscous flow[END_REF][START_REF] Dalibard | Effective boundary condition at a rough surface starting from a slip condition[END_REF][START_REF] Amirat | Effect of rugosity on a flow governed by stationary Navier-Stokes equations[END_REF].

Boundary layer approximation

The boundary layer approximation reads

u ε ≈ u P (x) + εv ε,bl ( x ε ) -εx 3 (V ε,∞ e, 0) (3.14)
where

• u P is the reference Poiseuille flow with no-slip condition at both planes, given by u P (x) = -x 3 (x 3 -1)e ;

• v ε,bl is the solution of the boundary layer problem

-∆v + ∇p = 0, in T 2 × R + , div v = 0, in T 2 × R + , v 3 = 0, y 3 = 0 , v = 0, y ∈ ε -1 T ε × {0}, ∂ y 3 v = -e, y ∈ ε -1 (T ε ) c × {0} . (3.
15) The boundary layer field v ε,bl has constant horizontal average v ε,∞ := T 2 v ε,bl (y) dy with respect to y 3 , and converges to (v ε,∞ , 0) exponentially fast as y 3 goes to infinity. Furthermore, by linearity of (3.15), one may denote v ε,∞ = V ε,∞ e for a two by two, symmetric positive definite matrix V ε,∞ ;

• the last term is a macroscopic Couette flow that corrects the "boundary layer constant" V ε,∞ e at infinity.

Averaging in the small scale, we find

u ε|x 3 =0 ≈ εV ε,∞ e, ∂ 3 u ε|x 3 =0 ≈ ∂ 3 u P |x 3 =0 ≈ e ,
ending up with the approximate boundary condition

u ε = εV ε,∞ ∂ 3 u ε on Γ b . (3.16)
On the basis of the previous reasoning, we implement the following strategy for the numerical computation of the slip length:

• Compute numerically (say with e = e 1 and e = e 2 ) the solution of (3.15), in order to determine the matrix V ε,∞ .

• Check for the asymptotics of εV ε,∞ , for various shapes and sizes of the no-slip zone T ε . This allows us to compare with theoretical results of Theorems 1 and 2. Indeed, sending ε to zero in (3.16) yields

ū = lim ε→0 (εV ε,∞ ) ∂ 3 ū on Γ b , (3.17)
so that the matrix M in the theorems satisfies

M -1 = lim ε→0 (εV ε,∞ ).
Numerical approximation of the matrix V ε,∞ In the numerical simulations, we solve the system (3.15) associated to different shapes of the no-slip zone T ε : circular or rectangular patches, and riblets parallel or orthogonal to the flow.

Notice that for such configurations, the matrix V ε,∞ is diagonal, so condition (3.16) simply reads

u ε,i = ε(V ε,∞ e i • e i ) ∂ 3 u ε,i on Γ b , for i = 1, 2. (3.18)
In the rest of this section, the quantity V ε,∞ e i • e i will be refered to as the average slip length in direction e i . It can be computed numerically by solving the cell problem (3.15) in a truncated domain T 2 × (0, H), using a finite element method. This step is very similar to the one described in Chapter 4, Section 4.5.

Computation of the average slip length

To illustrate our numerical results, we consider the case of patches, in the critical regime. We define the solid fraction φ s ε as the relative area of the no-slip zone T ε in the elementary square of size ε. In the case of patches, definition (3.4) yields φ s ε = aε ε 2 |T |, where |T | stands for the area of the domain T .

We have plotted V ε,∞ e 1 • e 1 against 1/ φ s ε , considering circular and squared patches in Fig. 3.3. We observe that the dependence is affine, and the coefficients obtained by linear regression are very close to the ones computed in [START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: Scaling laws for generic geometries[END_REF]. Consequently, since lim ε→0 φ s ε = 0,

V ε,∞ e 1 • e 1 ∼ α φ s ε as ε → 0. (3.19)
By Theorem 3.2.1, in the critical case a ε /ε 2 → C 0 > 0, there exists a 2×2 matrix M 0 (depending on the pattern T ) such that lim ε→0 εV ε,∞ = 1 C 0 M -1 0 . For circular or squared patterns centered in the unit square, by symmetry arguments, one obtains the existence of λ 0 > 0 such that M 0 = λ 0 0 0 λ 0 . Also,the following relation holds:

lim ε→0 εV ε,∞ e 1 • e 1 = 1 C 0 λ 0 .
Besides, using the definition of φ s ε in the case of patches, the asymptotic relation (3.19)

yields lim ε→0 εV ε,∞ e 1 • e 1 = α C 0 √ |T |
, which is consistent with Theorem 3.2.1.

Asymptotic analysis of a depletion layer model

In this section, we consider a completely different interpretation of the apparent slip length measured in nanoscopic devices, proposed in [START_REF] Myers | Why are slip lengths so large in carbon nanotubes? Microfluid[END_REF], where the author postulates that the source of this slip arises from a "depletion layer with reduced viscosity near the wall". This hypothesis is supported by experimental evidence [START_REF] Poynor | How water meets a hydrophobic surface[END_REF] and Molecular Dynamics simulations [START_REF] Joseph | Why are carbon nanotubes fast transporters of water?[END_REF] bringing out that the viscosity drops near the wall of the nanotube.

In [START_REF] Myers | Why are slip lengths so large in carbon nanotubes? Microfluid[END_REF], the flow is modeled as a Stokes flow in an infinite cylindrical pipe, with no-slip boundary condition, but two viscosities: a "bulk" one at the center of the pipe and a "wall" one near the walls which is smaller than the bulk one. By solving the equation in this simple geometry, the author shows that the resulting flow coincides with the one associated with an effective slip length and a constant viscosity equal to the one in the bulk.

From a mathematical perspective, the depletion layer model is strongly related to the reinforcement problems introduced by Sanchez-Palencia in [START_REF] Sanchez-Palencia | Problèmes de perturbations liés aux phénomènes de conduction à travers des couches minces de grande résistivité[END_REF], where an elastic medium is reinforced by the adjunction of a thin layer of very strong material. Brézis, Caffarelli and Friedman solved the interior and boundary reinforcement problems for elliptic equations, in the case of Dirichlet boundary conditions on a C 2 boundary and using strong solutions in [START_REF] Brézis | Reinforcement problems for elliptic equations and variational inequalities[END_REF]. A few years later, geometric measure theory and Gamma-convergence were successfully applied to boundary reinforcement problems (see for instance Acerbi and Buttazzo [START_REF] Acerbi | Reinforcement problems in the calculus of variations[END_REF], Buttazzo and Kohn [START_REF] Buttazzo | Reinforcement by a thin layer with oscillating thickness[END_REF], Buttazzo, Dal Maso and Mosco [START_REF] Buttazzo | Asymptotic Behaviour for Dirichlet Problems in Domains Bounded by Thin Layers[END_REF]).

As regards viscous flow models, Zhao and Yao recently studied in [START_REF] Zhao | Optimal boundary conditions for the Navier-Stokes fluid in a bounded domain with a thin layer[END_REF] the effects of reinforcement on strong solutions of Navier-Stokes equations, in a regular 2d domain. However, to our best knowledge, no results have yet been obtained in the case of a three-dimensional viscous flow with minimal regularity assumption on the boundary of the depletion layer, in the general framework of weak solutions. Let us present our contributions to this topic.

We have proposed in [B1, B2] a new approach based on a rescaling of the solution in the depletion layer, in the spirit of the unfolding method described in Chapters 4 and 5. Starting from natural energy bounds, we gain compactness on the rescaled velocity and pressure by a careful adaptation of arguments from [START_REF] Casado-Díaz | Asymptotic behavior of the Navier-Stokes system in a thin domain with Navier condition on a slightly rough boundary[END_REF], related to viscous flows in thin domains, and derive the Reynolds equation for the rescaled pressure. Then, we use a sequence of well-adapted test functions to determine the boundary condition satisfied by the limit velocity field, which depends, as expected, on the ratio between the value of the viscosity and the thickness of the depletion layer.

Drop of viscosity model

We assume that the depletion layer B ε is located in a small neighborhood of Γ b , of thickness of order ε, and set

B ε = {(x , x 3 ) ∈ Ω, 0 < x 3 < ε d(x )} where d ∈ W 1,∞ (T 2 ) is a positive function. We note Ω ε the bulk region Ω ε = Ω\B ε and Γ ε the interface defined by Γ ε = {(x , x 3 ) ∈ Ω, x 3 = ε d(x )}.
To the macroscopic variable x = (x , x 3 ) ∈ B ε , we associate the microscopic variable y = (y , y 3 ) defined by

y = x , y 3 = x 3 ε . (3.20)
The depletion layer is then described in microscopic variable by B ε = {(y , εy 3 ), (y , y 3 ) ∈ ω}, where ω is defined by

ω = (y , y 3 ) ∈ T 2 × (0, 1), 0 < y 3 < d(y ) .
Equations and boundary conditions Let u ε be the velocity of the fluid and p ε the pressure, and let f be a right-hand side defined on Ω. The viscosity µ ε is defined by

µ ε (x , x 3 ) = 1 Ωε (x , x 3 ) + ε α 1 Bε (x , x 3 ) (3.21)
where 1 Ωε and 1 Bε stand for the indicator functions of Ω ε and B ε respectively. We will work with Stokes equation

-div(2µ ε D(u ε )) + ∇p ε = f in Ω (3.22) div u ε = 0 in Ω (3.23)
completed with no-slip condition

u ε = 0 on Γ t , (3.24)
and periodic conditions on Γ . Since we aim at identifying the influence of the depletion layer B ε on the effective slippage, we start from the least favorable situation:

u ε = 0 on Γ b . (3.25)

Functional spaces and useful inequalities

Solutions of problem (3.22)-(3.25) will be defined in the space

V D = v ∈ H 1 0 (Ω, R 3 ), div v = 0 a.e. in Ω .
We will also use the space V N where the prescribed boundary condition on Γ b is non penetration:

V N = v ∈ H 1 (Ω, R 3 ), div v = 0 a.e. in Ω, v = 0 on Γ t , v 3 = 0 on Γ b .
In the above definitions, to lighten the notation, we do not use any special symbol to emphasize the (0, 1) 2 -periodicity in the horizontal variable.

Let us state some well-known inequalities that will be often used in the sequel, where C stands for a nonnegative constant independent on ε.

Poincaré inequalities Since the family Ω ε is uniformly bounded in the vertical direction, there holds for any

v ∈ H 1 (Ω ε ) such that v = 0 on Γ t , Ωε |v| 2 ≤ C Ωε |∇v| 2 .
(3.26)

In the domains B ε of thickness ε, there holds for any

v ∈ H 1 (B ε ) such that v = 0 on Γ b , Bε |v| 2 ≤ Cε 2 Bε |∇v| 2 . (3.27)

Korn inequalities

In Ω ε , a uniform Korn inequality holds since the domain is uniformly Lipschitz (see [START_REF] Bucur | Influence of wall roughness on the slip behaviour of viscous fluids[END_REF]):

for any v ∈ H 1 (Ω ε , R 3 ) s.t. v = 0 on Γ t , Ωε |∇v| 2 ≤ C K Ωε |D(v)| 2 . (3.28)
As a consequence, the following uniform Korn inequality is also valid in B ε :

for any v ∈ H 1 (B ε , R 3 ) s.t. v = 0 on Γ b , Bε |∇v| 2 ≤ C K Bε |D(v)| 2 .
(3.29)

Main results

Existence and uniqueness of u ε From standard theory on Stokes equation, we can rely on the following theorem:

Proposition 3.3.1. Let f be in L 2 (Ω, R 3
). There exists a unique function u ε in V D solution of the following problem

       -div (2µ ε D(u ε )) + ∇p ε = f in Ω div u ε = 0 in Ω u ε = 0 on Γ t ∪ Γ b (3.30) which means that for all φ ∈ V D , Ωε 2 D(u ε ) : D(φ) + ε α Bε 2 D(u ε ) : D(φ) = Ω f • φ . ( 3
u ε 2 L 2 (Ω,R 3 ) = Ωε |u ε | 2 + Bε |u ε | 2 ≤ C Ωε |D(u ε )| 2 + ε 2 Bε |D(u ε )| 2 ≤ C Ωε |D(u ε )| 2 + ε α Bε |D(u ε )| 2 ≤ C u ε L 2 (Ω,R 3 ) f L 2 (Ω,R 3 ) .
Hence, u ε L 2 (Ω,R 3 ) ≤ C and (3.32) is satisfied.

Our main result is the following.

Theorem 3.3.3. Let f ∈ L 2 (Ω, R 3 ), α ∈ (0, 2
) and (u ε ) ε be the family of solutions obtained from Proposition 3.3.1. Then u ε converges strongly in L 2 (Ω) toward the function ū solution of the system

       -div (2D(ū)) + ∇p = f in Ω div ū = 0 in Ω ū = 0 on Γ t (3.33)
completed with a boundary condition on Γ b that depends on α:

if 0 < α < 1 (subcritical case), ū = 0 on Γ b , (3.34)
if α = 1 (critical case),

u 3 = 0, 2[D(ū)e 3 ] tan + 1 d(x ) ū - 1 2 R(ū ) = 0 on Γ b , (3.35)
where R(ū ) is defined by the Reynolds equation (3.50) (see Section 3.3.3), if 1 < α < 2 (supercritical case),

u 3 = 0, 2[D(ū)e 3 ] tan = 0 on Γ b . (3.36)

Compactness implied by the energy bound

In this section, we consider a general sequence u ε ∈ V D satisfying the energy bound (3.32), and gather compactness results that are essential in the proof of Theorem 3.3.3.

Refining u ε in B ε

Combining (3.26), (3.28) and the energy bound, there holds

∀ε > 0 Ωε |u ε | 2 + |∇u ε | 2 ≤ C . (3.37)
However, due to the presence of the term ε α Bε |D(u ε )| 2 in the energy bound, u ε has no reason to be bounded in H 1 (Ω, R 3 ). In order to capture the behaviour of u ε at main order in ε, we refine it and construct a bounded sequence ūε in H 1 (Ω, R 3 ), such that for any ε > 0, ūε coincides with u ε in Ω ε . This step relies on extension operators in Sobolev spaces over Lipschitz domains (see for instance [START_REF] Chenais | On the existence of a solution in a domain identification problem[END_REF]).

Convergence of u

ε in L 2 (Ω, R 3 )
Existence of a limit ū ∈ V N By construction, the sequence (ū ε ) ε is bounded in H 1 (Ω, R 3 ). Hence, there exists a function ū in H 1 (Ω, R 3 ) such that, up to extraction, ūε converges to ū weakly in H 1 (Ω, R 3 ) and strongly in L 2 (Ω, R 3 ). By continuity of trace operators, it is easy to see that ū is at least in V N .

Strong convergence in L 2 Using that u ε is uniformly bounded in L 2 (Ω, R 3 ) and applying a density argument, one can deduce that for any α ∈ (0, 2], u ε converges weakly to ū in L 2 (Ω, R 3 ). Then, the strong convergence results from the inequality

Bε |u ε | 2 ≤ C ε 1-α 2 + √ ε 2-α + ε which implies that for α ∈ (0, 2), lim ε→0 Ω |u ε | 2 = Ω |ū| 2 .

Construction of the rescaled function v ε

In case α ≥ 1, a finer analysis of u ε is required in B ε . We define a rescaled function v ε = (v ε , v ε,3 ), depending on the micro-variable y = (y , y 3 ) ∈ ω by

v ε (y) = ε (α+3)/2 u ε (y , εy 3 ), v ε,3 (y) = -ε (α+1)/2 u ε,3 (y , εy 3 ) . (3.38)
By change of variable (3.20) and the energy bound, v ε satisfies

ω |∇ y v ε,3 | 2 + 1 ε 2 ω |∇ y v ε | 2 + |∂ y 3 v ε,3 | 2 + 1 ε 4 ω |∂ y 3 v ε | 2 ≤ C . (3.39)
Next proposition describes the asymptotic behaviour of v ε . It is obtained by studying the sequence v ε defined in H 1 (0, 1; L 2 ((0, 1) 2 )) by

v ε (z) = v ε (z , d(z )z 3 ) for z ∈ (0, 1) 3 , (3.40)
and arguing similarly as is the proof of [38, Lemma 5.1].

Proposition 3.3.4. Let (v ε ) be a sequence of functions in H 1 (ω, R 3 ) satisfying (3.39) and v ε|y 3 =0 = 0. Then, there exists w ∈ L 2 (ω, R 2 ) with ∂ 3 w ∈ L 2 (ω, R 2 ), and g ∈ H 2 (0, 1; H -1 ((0, 1) 2 )) such that, up to extraction:

(i) v ε 0 weakly in H 1 (ω, R 3 ); (ii) 1 ε v ε,3 0 and 1 ε ∂ 3 v ε,3 0 weakly in L 2 (ω); (iii) 1 ε 2 v ε w and 1 ε 2 ∂ 3 v ε ∂ 3 w weakly in L 2 (ω); (iv) 1 ε 2 v ε • (-∇ y d(z ), 1)
g weakly in H 2 (0, 1; H -1 ((0, 1) 2 )).

Moreover, w satisfies the incompressibility condition ∀η ∈ H 1 ((0, 1) The distinction between cases α < 1 and α ≥ 1 comes from the trace estimate

Γε |u ε | 2 ≤ C ε 1-α , (3.42
Γ b |ū ε | 2 ≤ C ε Bε |∇ū ε | 2 + Γε |ū ε | 2 .
Since ūε is bounded in H 1 (Ω, R 3 ) and ūε and u ε have the same trace on Γ ε , (3.42) implies that

Γ b |ū ε | 2 ≤ C ε + ε 1-α . Hence, if 0 < α < 1, lim ε→0 Γ b |ū ε | 2 = 0, so ū ∈ V D . Case α ≥ 1: ū3 = 0 on Γ b
By property (iv) in Proposition 3.3.4, taking the trace on z 3 = 1 and coming back to variable y, we get that the sequence of functions

y ∈ (0, 1) 2 → - 1 ε 2 v ε (y , d(y )) • ∇ y d(y ) + 1 ε 2 v ε,3 (y , d(y ))
is weakly convergent in H -1 ((0, 1) 2 ). Observing that, by a trace argument, the sequence y ∈ (0, 1)

2 → 1 ε 2 v ε (y , d(y )) is bounded in L 2 ((0, 1) 2 , R 2 ), this implies that 1 ε 2 v ε,3 (y , d(y )
) is also bounded in H -1 ((0, 1) 2 ). By definition (3.38), this means that the sequence

y ∈ (0, 1) 2 → ε (α-3)/2 u ε,3 (y , εd(y ))
is bounded in H -1 ((0, 1) 2 ). As a result, for any α ≤ 2, u ε,3 (•, εd(•)) converges to 0 in H -1 ((0, 1) 2 ), and since this sequence is in L 2 ((0, 1) 2 ), it implies that for any η ∈ C 1 c ((0, 1) 2 ), lim ε→0 (0,1) 2 u ε,3 (y , εd(y ))η(y ) dy = 0 .

Now we are in position to prove that ū3 = 0 on Γ b . Let η = η(y ) ∈ L 2 ((0, 1) 2 ). Integrating on vertical lines yields Observing that the right-hand side of the previous equality converges to 0, and using (3.43), we conclude that ū3 = 0 on Γ b .

Extra information on w in the critical case α = 1

In the critical case, the asymptotic behaviours of the v ε and ū ε are coupled, in the sense that w (y , d(y )) = ū (y , 0) for a.e. y ∈ (0, 1) 2 .

(3.44)

Reynolds equation in critical case α = 1

Existence and rescaling of the pressure By classical results (see for instance [59, paragraph 4.1]), one can associate to the velocity field u ε solution to (3.30), a pressure p ε ∈ L 2 (Ω) (uniquely defined up to a constant) such that

∀ϕ ∈ H 1 0 (Ω, R 3 ) Ω 2µ ε D(u ε ) : D(φ) - Ω p ε div φ = Ω f • φ . (3.45)
Defining q ε : ω → R by q ε (y) = ε (α+1)/2 p ε (y , ε y 3 ) for a.e. y ∈ ω ,

there holds ω |q ε (y

)| 2 dy = ε α Bε |p ε | 2 ≤ C.
Therefore, there exists q ∈ L 2 (ω) such that, up to extraction, q ε q weakly in L 2 (ω) .

(3.47)

Derivation of Reynolds equation on q

Let φ ∈ H 1 0 (ω, R 3 ) and define

φ ε ∈ H 1 0 (B ε , R 3 ) by φ ε (x) = φ (x , x 3 ε ) and φ ε,3 (x) = -ε φ 3 (x , x 3 ε ) .
Setting φ ε = 0 in Ω ε and testing against φ ε in (3.45), we obtain

2ε Bε D(u ε ) : D(φ ε ) - Bε p ε div φ ε = Bε f • φ ε .
Using the convergences established in Sections 3.3.2 and 3.3.3, we can pass to the limit in the previous relation and establish that q = q(y ) and satisfies

-∂ 2 3,3 w + ∇ y q = 0 in H -1 (ω, R 2 ) .
Integrating this equation in the y 3 variable and using that w = 0 on y 3 = 0 and condition (3.44), we get the following expression of w :

w (y) = - y 3 2 (d(y ) -y 3 )∇ y q + y 3 d(y ) ū (y , 0) . (3.48)
Coming back to (3.41), we get Reynolds equation on q:

∀η ∈ H 1 (T 2 ) (0,1) 2 d(y ) 3 6 ∇ y q • ∇ y η dy = (0,1) 2 d(y ) ū (y , 0) • ∇ y η dy .
(3.49) We set R(ū ) = ∇ y q and more generally, R(v ) = ∇ y r if r is the solution to the Reynolds equation associated to v :

∀η ∈ H 1 (T 2 ) (0,1) 2 d(y ) 3 6 R(v ) • ∇ y η dy = (0,1) 2 d(y ) v (y , 0) • ∇ y η dy .
(3.50)

Identification of the limit problem

Using the elements gathered in the previous sections, we are able to characterize ū as the solution to the limit problem (3.33) completed with the boundary condition (3.34), (3.35) or (3.36).

Subcritical case α < 1

Since we have already established that ū = 0 on Γ b , it remains to prove that for any

φ ∈ V D , Ω 2D(ū) : D(φ) = Ω f • φ . (3.51)
By a density argument, one may assume that φ is supported in Ω ε , so by definition of ūε , the variational formulation (3.31) reduces to

Ω 2D(ū ε ) : D(φ) = Ω f • φ .
(3.51) follows by weak convergence of ūε to ū in H 1 0 (Ω, R 3 ).

Critical and supercritical case

α ≥ 1 Let φ ∈ V N and φ ∈ H 1 (ω, R 2 ) such that φ (y , 0) = -φ (y , 0) φ (y , d(y )) = 0
for y ∈ (0, 1) 2 .

(3.52)

By conditions (3.52), the function φ ε defined by

   φ ε (x) = φ (x) + 1 Bε (x) φ (x , x 3 ε ) φ ε,3 (x) = φ 3 (x) for x ∈ Ω , (3.53)
is in H 1 0 (Ω, R 3 ). Hence, we can use it as test function in (3.45) and obtain

2 Ωε D(u ε ) : D(φ) + 2ε α Bε D(u ε ) : D(φ) + 2ε α Bε ∇ x (u ε ) : ∇ y ( φ ) - 1 2 (∂ x 3 u ε + ∇ x u ε,3 ) • 1 ε ∂ y 3 φ - Bε p ε div y φ = Ω f • φ + Bε f • φ . (3.54)
In the previous relation, functions u ε , φ, f and their derivatives are computed at x, function φ and its derivatives with respect to y are computed at (x , x 3 ε ), and all integrals are performed with respect to x.

Passing to the limit in (3.54)

• Using that u ε coincides with ūε in Ω ε and the energy bound (3.32), lim • We use the change of variables (3.20) and the rescaled functions v ε and q ε to turn integrals in x ∈ B ε into integrals in y ∈ ω.

-By (3.47), there holds lim ε→0 Bε

p ε div y φ = lim ε→0 ω q ε div y φ = ω q div y φ . -We have 2ε α Bε ∇ x u ε : ∇ y φ - 1 2 (∂ x 3 u ε + ∇ x u ε,3 ) • ( 1 ε ∂ y 3 φ ) = 2ε α-1 ω ∇ y v ε : ∇ y φ + 1 2ε 2 ∂ y 3 v ε • ∂ y 3 φ + 1 2 ∇ y v ε,3 • ∂ y 3 φ ,
and by Proposition 3.3.4,

lim ε→0 2 Bε ∇ x (u ε ) : ∇ y ( φ ) - 1 2 (∂ x 3 u ε + ∇ x u ε,3 ) • ( 1 ε ∂ y 3 φ ) = ω ∂ y 3 w • ∂ y 3 φ . (3.55)
Gathering the previous results, we can pass to the limit in (3.54) and obtain:

• in the supercritical case α > 1, 2 Ω D(ū) : D(φ) = Ω f • φ .
This proves that ū is the variational solution to the system (3.33), completed with the boundary condition (3.36) on Γ b ;

• in the critical case α = 1, 2 Ω D(ū) : D(φ) + ω ∂ y 3 w • ∂ y 3 φ - ω q div y φ = Ω f • φ . (3.
-∂ y 3 w (y , 0) • φ (y , 0) dy = - ω ∇ y q • φ + (0,1) 2 ∂ y 3 w (y , 0) • φ (y , 0) dy = ω q div y φ + (0,1) 2 1 d(y ) ū (y , 0) - d(y ) 2 ∇ y q(y ) • φ (y , 0) dy .
Coming back to relation (3.56), we conclude that ū is the variational solution to the system (3.33) completed with boundary condition (3.35).

Conclusion and perspectives

In this chapter, we have performed a rigorous asymptotic analysis of two models introduced in fluid mechanics reviews, accounting on different mechanisms to explain the occurence of slip in micro-and nanofluidics and quantify the associated slip lengths:

(a) the Cassie state resulting from superhydrophobicity properties of the surface (Section 3.2) ;

(b) a potential drop of viscosity near the solid wall (Section 3.3).

One main goal of this analysis is to link the mathematical results with the experimental data available in fluid mechanics literature, or with the numerical simulations provided by the teams of researchers at the origin of these models. This approach has guided our numerical study of the Cassie state model (Section 3.2.4), where for instance we have retrieved the same behaviour of the average slip length as a function of the solid fraction as in [START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: Scaling laws for generic geometries[END_REF] (Fig. 3.

3).

As regards the drop of viscosity model, an ongoing work with J. Olivier suggests that the effective slip length identified in Theorem 3.3.3 is compatible with recent experimental results. However, many experiments are performed on axisymmetric domains, such as carbon nanotubes, which is why we are currently working on an adaptation of the model to this type of geometry. Furthermore, the fact that the method relies only on properties of weak solutions suggests that it could be adapted to more general settings, such as time-dependent Navier-Stokes equations or non-Newtonian flows.

Finally, a mid-term objective could be to conceive and analyse a mixed model between models (a) and (b), taking into account a viscosity drop and superhydrophobicity effects as the same time, or more generally, introducing anisotropy in the width of the depletion layer due to geometric irregularities of the wall.

Chapter 4

Roughness effects on Large Eddy Simulation for turbulent flows

This chapter presents results obtained in collaboration with F. J. Suárez-Grau in [B9].

Introduction

The study of roughness effects in fluid-structure interaction has seen considerable development in the last decades. Its main objective is to explain how fine structures on solid walls may influence viscous flows, and how the effects of surface irregularities should be taken into account to enhance the precision of existing models.

One of the first mathematical results in this direction concerns the justification of the no-slip boundary condition. In their pioneering paper [START_REF] Casado-Díaz | Why viscous fluids adhere to rugose walls: A mathematical explanation[END_REF], Casado-Díaz et al. have explained that boundary condition using an asymptotic analysis based on the presence of small-scale irregularities on the wall, and on a minimal boundary condition -the no-penetration condition -expressing the impossibility for the fluid to get through the solid. In a simplified manner, this phenomenon could be interpreted as follows: the rough geometrical structure of the wall causes the normal ν to the boundary to oscillate with high spatial frequency, turning no-penetration condition u•ν = 0 into no-slip condition u = 0, where u stands for the fluid velocity. This re-sult has then been generalized by many authors, leading to a rather complete description of the possible roughness effects on viscous flows. We refer to the Introduction of Chapter 3 for references on this topic.

Among periodical roughness profiles, the riblets, which are assemblies of the same elongated structure, are of considerable interest. They are widely used in engineering, extensively studied in experiments and mechanical models (see [START_REF] Bechert | The viscous flow on surfaces with longitudinal ribs[END_REF][START_REF] Jiménez | Drag reduction by riblets[END_REF][START_REF] Lee | Flow field analysis of a turbulent boundary layer over a riblet surface[END_REF]) and have been the subject to many mathematical studies (see, for instance, [START_REF] Amirat | Riblets and drag minimization[END_REF][START_REF] Bucur | Influence of wall roughness on the slip behaviour of viscous fluids[END_REF][START_REF] Jäger | Couette flows over a rough boundary and drag reduction[END_REF]). However, in many applications such as aeronautics [START_REF] Viswanath | Aircraft viscous drag reduction using riblets[END_REF] or shark-skin drag reduction [START_REF] Dean | Shark-skin surfaces for fluid-drag reduction in turbulent flow: a review[END_REF], riblet effects occur in turbulent regime.

At present, a rigourous mathematical analysis of fluid-structure interaction at high Reynolds number remains out of reach, due to the well-known, tremendous difficulties posed by the three-dimensional Navier-Stokes equations in that context. From a numerical point of view, their direct numerical simulation at high Reynolds number can be very expensive, since it requires solving a very wide range of time and length scales -whose presence is characteristic of turbulence. As a result, lots of tools and models have been developed to describe and simulate turbulent flows, that are commonly used in mechanics [START_REF] Pope | Turbulent Flows[END_REF].

We focus here on the class of Large Eddy Simulation (LES) models, that were introduced by Smagorinsky in 1963 [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF] and are widely used in computational fluid dynamics. The principle of LES is to apply a spatial convolution operator to the Navier-Stokes equations, in order to decompose the velocity, pressure and body forces into two parts: an averaged part describing the large scale behaviour of the fluid, and a remainder corresponding to the small scales. The goal of LES models is to compute only the large scale structure of the flow. However, due to the presence of the nonlinear convective term in the moment equation, the interaction between small and large scales needs to be modeled in order to close the system. This well-known closure problem has led to a large variety of LES models, whose derivation can be motivated by physical or mathematical considerations.

In this chapter, we address the simulation of turbulent flow of a viscous fluid, in the vicinity of a rough wall, using one of the simplest and most popular LES model: the Smagorinsky model. As will be presented below, this model involves a nonlinear stress tensor containing a term of power law type, and thus relates to models of quasi-Newtonian viscous flows, that were studied in the context of periodic roughness effects in [START_REF] Suárez-Grau | Effective boundary condition for a quasi-Newtonian viscous fluid at a slightly rough boundary starting from a Navier condition[END_REF][START_REF] Gérard | Boundary layer for a non-Newtonian flow over a rough surface[END_REF].

As was proved by Coletti in [START_REF] Coletti | Analytical and numerical analysis of k -e and Large eddy simulation turbulence models[END_REF], Smagorinsky model enjoys nice math-ematical properties, in particular, global existence, uniqueness and stability of weak solutions. In order to develop an asymptotic analysis of the model, with respect to the characteristic length of the rough pattern, we rely on a remark by Galdi and Layton. Starting from the observation that the no-slip boundary condition, which is the commonly accepted boundary condition for viscous flows at high Reynolds number, cannot be preserved by spatial filtering, they have proposed in [START_REF] Galdi | Approximation of the Larger Eddies in Fluid Motions II: a Model for Space-Filtered Flow[END_REF] to equip LES models with Navier boundary conditions. As stated in the beginning of this Introduction, such boundary conditions are well adapted to the study of rugosity effects. As a matter of fact, we are able to generalize the multi-scale analysis developed in [START_REF] Casado-Díaz | Asymptotic behavior of a viscous fluid with slip boundary conditions on a slightly rough wall[END_REF], to a timedependant system of equations, involving a nonlinear stress tensor, in the case of a periodic rough wall. We identify a new critical scaling between the amplitude and the period of the oscillations, and obtain the effective system that expresses the effect of the rough patches on the LES model, by adding a nonlinear interaction term to the initial boundary condition of Navier type. The homogenized boundary condition is determined locally by solving a boundary layer problem that involves only the nonlinear part of the stress tensor -that is, the part modelling the large scales/small scales interactions. Finally, we propose and implement a numerical method to simulate the steady-state of the LES model, using three different riblet shapes that are commonly studied in the literature, and discuss the influence of the ribbed geometry on the simulations.

Smagorinsky model for turbulent flow

We consider the flow of a viscous and homogeneous fluid, driven by an external force, in a three-dimensional domain delimited by two fixed horizontal plates. We assume that the upper plate is perfectly smooth, but that the lower one presents a periodic rough surface. We work in dimensionless coordinates and denote by Re the Reynolds number associated with the flow.

Let us start with a description of the geometry.

Geometric notation

For the sake of simplicity, we assume that the flow is periodic with respect to the horizontal coordinates x = (x 1 , x 2 ). To describe the lower (rough) plate, we fix a positive function Ψ ∈ W 2,∞ (T 2 ) and define

Ψ ε (x 1 , x 2 ) = a ε Ψ x 1 ε , x 2 ε ,
where ε > 0 is a small parameter and a ε satisfies

lim ε→0 a ε = 0, lim ε→0 a ε ε = 0. (4.1)
The rough boundary Γ ε is then defined by

Γ ε := {(x 1 , x 2 , x 3 ) ∈ T 2 × R : x 3 = -Ψ ε (x 1 , x 2 )}. (4.2)
By periodicity of Ψ, Ψ ε is periodic of period ε in the x 1 and x 2 directions, and by construction, the maximum value of Ψ ε is proportional to a ε . Thus, hypotheses (4.1) mean that the amplitude of the roughness is negligible with respect to the period, as the period tends to zero. Following the notation from Chapter 3, we denote by Γ t := T 2 × {1} the upper (smooth) boundary, and introduce the fluid domain Ω ε defined by

Ω ε := {(x 1 , x 2 , x 3 ) ∈ T 2 × R : -Ψ ε (x 1 , x 2 ) < x 3 < 1}.
Since the sequence (Ψ ε ) ε>0 converges uniformly to zero, the sequence (Γ ε ) ε>0 converges in the sense of Hausdorff to the flat wall Γ b defined by Γ b = T 2 × {0}. Accordingly, we introduce the limit domain Ω = T 2 × (0, 1) and define the box Λ = T 2 × (-1, 1), which contains all the rough domains Ω ε .

Derivation of the LES model

Let u ε be the velocity of the fluid and p ε the pressure. We assume that the external force is given by the restriction to Ω ε of a vector field f ∈ L 2 (Λ, R 3 ). The fluid motion is described by Navier-Stokes equations

∂ t u ε + div(u ε ⊗ u ε ) - 1 Re ∆u ε + ∇p ε = f in (]0, T ) × Ω ε , (4.3) div u ε = 0 in (0, T ) × Ω ε , (4.4)
completed with the initial condition

u ε (x, 0) = u 0 ε (x), x ∈ Ω ε (4.5)
and no-slip boundary conditions on the plates:

u ε = 0 on Γ t , (4.6
)

u ε = 0 on Γ ε . (4.7)
As is well known, in dimension 3, existence and uniqueness of the solution of (4.3)-(4.5) has only been proved for small Reynolds number [START_REF] Ladyzhenskaya | On the existence and uniqueness of the solution of the nonstationary problem for a viscous, incompressible fluid[END_REF]. For large values of Re, solutions of Navier-Stokes equations become unstable; this is called turbulence. In such regime, direct numerical simulations require a very fine mesh of the space domain, which dramatically increases the storage and computational costs, especially when the boundary presents fast oscillations.

As an alternative to direct numerical simulation, several models of turbulent flows have been developed [START_REF] Mcdonough | Introductory lectures on turbulence. Physics, mathematics and modeling[END_REF]. We focus her on the space filtering approach to large eddy simulation (LES), introduced in 1970 by Deardorff [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF]. Let us recall the main steps of the derivation of this model.

Space filtering approach

The spatial averaging is defined as a convolution with the Gaussian filter

g δ (x) := γ π 3/2 1 δ 3 exp(-γ|x| 2 /δ 2 ) ,
where δ > 0 represents an averaging radius and γ > 0 is a constant. All dependent variables being extended by zero to T 2 × R, the average of any variable r(t, x) is defined by

r(t, x) = (g δ * r)(t, x) = R 3 g δ (x -y) r(t, y) dy.
Applying this averaging operator to system (4.3)-(4.5) and commuting spatial partial derivatives with convolutions leads to the system

∂ t u ε + div(u ε ⊗ u ε ) - 1 Re ∆u ε + ∇p ε = f in (0, T ) × Ω ε , (4.8) div u ε = 0 in (0, T ) × Ω ε , (4.9)
u ε (x, 0) = u 0 ε (x), x ∈ Ω ε . (4.10)
However, the nonlinear term div(u ε ⊗ u ε ) cannot be expressed as a function of u ε only. This notorious closure problem forbids to solve the averaged problem (4.8)-(4.10) directly. It is thus required to derive an approximate expression of u ε ⊗ u ε , involving only the average velocity u ε .

Smagorinsky approximation

We use Smagorinsky approximation [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF]:

u ε ⊗ u ε ≈ u ε ⊗ u ε + δ 2 2γ ∇u ε (∇u ε ) T -C s |∇u ε |∇u ε ,
where C s > 0 is of order O(δ 2 ).

The LES model

Defining (w ε , q ε ) as an approximation of (u ε , p ε ), we obtain the system

∂ t w ε + (w ε • ∇)w ε -div S(∇w ε ) + ∇q ε = f in (0, T ) × Ω ε , (4.11)
div w ε = 0 in (0, T ) × Ω ε , (4.12)

w ε (x, 0) = w 0 ε (x), x ∈ Ω ε , (4.13)
where the stress tensor S(∇w ε ) is defined by

S(∇w ε ) = 1 Re + C s |∇w ε | ∇w ε - δ 2 2γ ∇w ε (∇w ε ) T . (4.14)
In (4.13), the initial datum w 0 ε is obtained by setting w 0 ε (x) = u 0 ε (x) a.e. in Ω ε .

Treatment of the boundary conditions

To complete system (4.11)-(4.13), we need to impose boundary conditions. Since we want to focus on the interaction between the fluid and the rough boundary Γ ε , for simplicity, we impose no-slip boundary condition on the top boundary:

w ε = 0 on Γ t . (4.15)
Concerning the boundary condition on Γ ε , as mentionned in the Introduction, we impose a condition of Navier type, that reads in the context of the nonlinear stress tensor (4.14)

w ε • ν = 0 on Γ ε , [S(∇w ε )ν] tan + β [w ε ] tan = 0 on Γ ε i = 1, 2. (4.16)
Remark 4.2.1. In the tangential condition, β = β(Re, δ) is a nonnegative friction coefficient, that depends typically on the Reynolds number Re and on the averaging radius δ. In particular, the model should impose that lim δ→0 β = +∞, so that formally, conditions (4.16) turn into no-slip conditions as the averaging radius vanishes, recovering (4.7).

In the present study, we make the additional assumption that β(Re, δ) is independant on ε. In fact, as long as β remains bounded when ε goes to zero, it does not play any significant role in the asymptotic analysis of system (4.11)-(4.16).

Remark 4.2.2. In the computations, we will often divide the tensor S(∇w ε ) in two parts: the Newtonian and non-Newtonian parts, setting

S(∇w ε ) = 1 Re ∇w ε Newtonian part + S(∇w ε ) non-Newtonian part , with ∀ζ ∈ R 3×3 S(ζ) = C s |ζ|ζ - δ 2 2γ ζζ T (4.17)
Weak formulation of system (4.11)-(4.16).

We introduce the functional spaces

H ε = L 2 (Ω ε , R 3 ) , W ε = w ∈ W 1,3 (Ω ε , R 3 ), w • ν = 0 on Γ ε , w = 0 on Γ t , V ε = {w ∈ W ε , div w = 0 in Ω ε } ,
and define

J ε = H 1 (0, T ; H ε ) ∩ L 3 (0, T ; V ε ) . (4.18)
Assume that w 0 ε is given in V ε . We say that w ε is a weak solution to problem (4.11)-(4.16) if w ε ∈ J ε , satisfies initial condition (4.13)and if the integral identity

T 0 Ωε [∂ t w ε + (w ε • ∇)w ε ] • ϕ + S(∇w ε ) : ∇ϕ dx dt + β T 0 Γε w ε • ϕ dH 2 dt = T 0 Ωε f • ϕ dx dt (4.19)
holds for every test function ϕ such that

ϕ ∈ C 1 c ((0, T ) × Ω ε , R 3 ), div ϕ = 0 in Ω ε , ϕ • ν = 0 on Γ ε , ϕ = 0 on Γ t .
Notice that, by Sobolev embedding theorem, the space J ε is continuously embedded in C 0 ([0, T ], H ε ), so that initial condition (4.13) is well-defined.

Main results

The first result ensures the existence and uniqueness of a weak solution w ε ∈ J ε to problem (4. [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF], and states the expected energy estimates. It is adapted from Coletti [START_REF] Coletti | Analytical and numerical analysis of k -e and Large eddy simulation turbulence models[END_REF].

Theorem 4.3.1. Let T > 0 and f ∈ H 1 (0, T ; L 2 (Λ, R 3 )). Assume that C s > δ 2 γ .
Then for every ε > 0 and every initial datum w 0 ε ∈ V ε , there exists a unique weak solution w ε ∈ J ε to problem (4.11)-(4.16). Moreover, there exists a constant C = C(T, Re, C s , δ, γ) such that this solution satisfies the energy inequalities:

sup 0≤t≤T Ωε |w ε | 2 dx ≤ Ωε |w 0 ε | 2 dx + C T 0 Ωε |f | 2 dx dt , (4.20) 1 Re T 0 Ωε |∇w ε | 2 dx dt + 2 T 0 Γε β|w ε | 2 dH 2 dt + 2 C s - δ 2 2γ T 0 Ωε |∇w ε | 3 dx dt ≤ Ωε |w 0 ε | 2 dx + C T 0 Ωε |f | 2 dx dt , (4.21) T 0 Ωε |∂ t w ε | 2 + |∇(∂ t w ε )| 2 dx dt + T 0 Γε |∂ t w ε | 2 dH 2 dt ≤ C Ωε |∂ t w 0 ε | 2 + |∇w 0 ε | 3 dx + T 0 Ωε |f | 2 + |∂ t f | 2 dx dt . (4.22)
In the previous energy estimate, ∂ t w 0 ε is defined by

∂ t w 0 ε = -(w 0 ε • ∇)w 0 ε + div S(∇w 0 ε ) -∇q 0 ε + f |t=0 (4.23)
where q 0 ε is the solution to the problem

-∆q 0 ε = div[(w 0 ε • ∇)w 0 ε -div S(∇w 0 ε )] in Ω ε , (4.24
)

∂q 0 ε ∂ν = (div S(∇w 0 ε )) • ν on Γ t ∪ Γ ε . (4.25)
In the next proposition, we construct a pressure field q ε associated with the divergence free constraint on w ε , relying on a classical argument by Tartar [START_REF] Tartar | An Introduction to Navier-Stokes Equation and Oceanography[END_REF].

Proposition 4.3.2. Let T > 0, f ∈ H 1 (0, T ; L 2 (Λ, R 3 )) and assume that C s > δ 2
γ . Let ε > 0, w 0 ε ∈ V ε and let w ε ∈ J ε be the weak solution to problem (4.11)-(4.16). Then, there exists

q ε ∈ W -1,∞ (0, T ; L 3/2 0 (Ω ε )) such that for every ϕ ∈ C 1 c (0, T ; W ε ), T 0 Ωε ∂ t w ε + (w ε • ∇)w ε • ϕ + S(∇w ε ) : ∇ϕ -q ε div ϕ dx dt +β T 0 Γε w ε • ϕ dH 2 dt = T 0 Ωε f • ϕ dx dt . (4.26)
Our main contribution is the obtention of the asymptotic behaviour of the velocity/pressure couple (w ε , q ε ) satisfying (4.26), as ε tends to zero. The statement of the effective system requires some extra notation to define the auxiliary (cell) problem.

Notation related to the cell problem.

We define:

• Y = (-1 2 , 1 2 ) 2 the unitary cube of R 2 ; • Q = Y × (0, +∞); • Q M = Y × (0, M ), for M > 0.
We use the index # to mean periodicity with respect to Y , for example L p # (Y ) denotes the space of functions u ∈ L p loc (R 2 ) which are Y -periodic, while L p # ( Q) denotes the space of functions u ∈ L p loc (R 2 × (0, +∞)) such that Q | u| p dy < +∞, u(y +k , y 3 ) = u(y), ∀k ∈ Z 2 , a.e. y ∈ R 2 ×(0, +∞).

We denote by V the space of functions v : R 2 × (0, +∞) → R such that

∀M > 0 v ∈ W 1,3 # ( Q M , R) , and ∇ v ∈ L 3 # ( Q, R 3 ) .
Asymptotic behaviour of (w ε , q ε ).

Depending on the scaling of the amplitude a ε with respect to ε, the asymptotic behaviour of system (4.11)-(4.16) is described by three different regimes, that we call critical, subcritical and supercritical. The relevant scaling is controlled by parameter λ, which is a nonnegative number, possibly infinite, defined by λ := lim We state the result in the critical case; the other cases will be discussed in Remark 4.3.2.

Theorem 4.3.3. Let T > 0, f ∈ H 1 (0, T ; L 2 (Λ, R 3 )) and assume that C s > δ 2
γ . We consider a family of initial data w 0 ε ∈ V ε , satisfying the uniform bounds

∀ε > 0 w 0 ε Wε ≤ C , div S(∇w 0 ε ) L 2 (Ωε,R 3 ) ≤ C . (4.28)
Without loss of generality, we may assume the existence of w 0 ∈ W 1,3 (Ω, R 3 ) such that w 0 ε|Ω converges weakly to w 0 in W 1,3 (Ω, R 3 ). For every ε > 0, let w ε ∈ J ε be the weak solution to problem (4.11)-(4.16) and q ε ∈ W -1,∞ (0, T ; L

3/2 0 (Ω ε )) be the associated pressure. Assume a ε ≈ ε 5 3 with a ε /ε 5 3 → λ, 0 < λ < +∞.
Then, there exists (w, q) such that

w ε|Ω w in L 3 (0, T ; W 1,3 (Ω, R 3 )) , q ε|Ω q in W -1,∞ (0, T ; L 3 2 (Ω)) .
(4.29) Moreoever, (w, q) satisfies the effective problem

         ∂ t w + (w • ∇)w -div (S(∇w)) + ∇q = f in (0, T ) × Ω, div w = 0 in (0, T ) × Ω, w(0, x) = w 0 (x), x ∈ Ω, w = 0 on (0, T ) × Γ t , (4.30)
completed with the boundary conditions:

w 3 = 0 on (0, T ) × Γ b , -S(∇w) 1,3 + β w 1 + λ 3 R 1 (w ) = 0 on (0, T ) × Γ b , -S(∇w) 2,3 + β w 2 + λ 3 R 2 (w ) = 0 on (0, T ) × Γ b . (4.31)
Here, R i : R 2 → R, ξ → R i , i ∈ {1, 2}, are defined by

R i (ξ ) = Q S(∇ y φ ξ ) : ∇ y φ e i dy ∀ ξ ∈ R 2 , (4.32)
where for every ξ ∈ R 2 , ( φ ξ , π ξ ) is the solution of

               -div y ( S(∇ y φ ξ )) + ∇ y π ξ = 0 in T 2 × R + , div y φ ξ = 0 in T 2 × R + , φ ξ 3 (y , 0) = -∇Ψ(y ) • ξ on T 2 × {0}, S(∇ y φ ξ ) i,3 = 0, i = 1, 2, on T 2 × {0}. (4.33)
Remark 4.3.1. The roughness effect depends only on the non-Newtonian part S of the stress tensor, that describes the large-scale/small-scale interactions in Smagorinsky's LES model, and of the geometry of the riblets, specified by the function Ψ.

Remark 4.3.2. Theorem 4.3.3 can be easily adapted to cover the two other asymptotic regimes.

• In the subcritical case a ε /ε 5 3 → 0, the result holds replacing λ by 0 in conditions (4.31). In this regime, the roughness has no effect on the boundary condition satisfied by the homogenized velocity field w.

• In the supercritical case a ε /ε 5 3 → +∞, boundary conditions (4.31) are replaced by no-slip conditions w = 0 on (0, T ) × Γ b . In that case, the roughness effect is maximal.

The following proposition ensures the existence and uniqueness of the solution to system (4.33) and states some exponential decay properties in the y 3 -direction that will be useful to perform numerical simulations. We refer to [START_REF] Casado-Díaz | Exponential decay for the solutions of nonlinear elliptic systems posed in unbounded cylinders[END_REF] for a proof.

Proposition 4.3.4. Assume that C s -δ 2 γ > 0. For every ξ ∈ R 2 , there exists a unique solution ( φ ξ , π ξ ) ∈ (V/R) 2 × V × L 3 2 # ( Q, R) to system (4.33).
Moreover, there exist constants C, τ > 0 such that, for every z 3 > 0,

Q |∇ y φ ξ | 3 dy ≤ C |ξ | 3 T 2 |∇Ψ| 3 dy , (4.34)
T 2 ×(z 3 ,+∞)

|∇ y φ ξ | 3 dy ≤ C |ξ | 3 e -τ z 3 T 2 |∇Ψ| 3 dy . (4.35)
The particular case of riblets.

In a wide variety of applications to fluid mechanics, rough boundaries are composed of ribbed walls, which be modeled by Ψ ε (x ) = a ε Ψ(x 1 /ε), where Ψ is a positive one-variable function, periodic with period 1. In this particular setting, the description of the effective boundary condition can be simplified.

Corollary 4.3.5. Under the same hypotheses as in Theorem 4.3.3, if we assume in addition that Ψ ε (x ) = a ε Ψ(x 1 /ε), then boundary conditions (4.31) are replaced by

w 3 = 0 on (0, T ) × Γ b , -S(∇w) 1,3 + β w 1 + λ 3 w 2 1 R 1 (sgn(w 1 )e 1 ) = 0 on (0, T ) × Γ b , -S(∇w) 2,3 + β w 2 = 0 on (0, T ) × Γ b , (4.36)
where sgn is the sign function sgn(0) = 0, sgn(s) = s |s| if s ∈ R \ {0}, and R 1 (•) is given by (4.32). The proof of Theorem 4.3.1 is an adaptation of the original proof by Coletti. Its essential ingredients are the continuity, coerciveness and monotonicity of the operator S : R 3×3 → R 3×3 defined by (4.17), that we summarize here:

Elements of proof

∀ζ ∈ R 3×3 | S(ζ)| ≤ (C s + δ 2 2γ )|ζ| 2 , ∀ζ ∈ R 3×3 S(ζ) : ζ ≥ (C s - δ 2 2γ )|ζ| 3 , ∀ζ 1 , ζ 2 ∈ R 3×3 ( S(ζ 1 ) -S(ζ 2 )) : (∇ζ 1 -∇ζ 2 ) ≥ 0 . (4.37)

Proof of Theorem 4.3.3

In order to prove Theorem 4.3.3, we rely on an adaptation of the unfolding method (see e.g. [START_REF] Arbogast | Derivation of the double porosity model of single phase flow via homogenization theory[END_REF][START_REF] Casado-Díaz | Asymptotic behavior of a viscous fluid with slip boundary conditions on a slightly rough wall[END_REF][START_REF] Cioranescu | Periodic unfolding and homogenization[END_REF]), which is strongly related to the two-scale convergence method [START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF]. This approach requires some extra notation.

Extra notation related to the unfolding method

For k ∈ Z 2 and ρ > 0, we introduce the sets

C k ρ = ρk + ρY , Q k ρ = Ω ε ∩ (C k ρ × R) .
Then, we define κ : R 2 → Z 2 by

κ(x ) = k ⇔ x ∈ C k 1 .
κ is well defined up to a set of measure zero in R 2 . Moreover, for every ρ > 0 and k ∈ Z 2 , we have

κ x ρ = k ⇔ x ∈ C k ρ .
Finally, we set

K ε = k ∈ Z 2 : (0, 1) 2 ∩ C k ε = ∅ .
Localisation procedure To capture the behaviour of a sequence of functions, say, r ε defined on (0, T ) × Ω ε , in the vicinity of the oscillating wall Γ ε , the unfolding method relies on a localisation procedure that consists in defining the "unfolded" function r ε by

r ε (t, x , y) = r ε t, εκ( x ε ) + εy , εy 3 , (4.38)
for (t, x , y)

∈ (0, T )×T 2 × Y ε , with Y ε = y ∈ Y × R : -aε ε Ψ(y ) < y 3 < 1 ε .
Remark 4.4.1. In particular, for k ∈ Z 2 the restriction of r ε to (0, T )×C k ε × Y ε (i) does not depend on x ;

(ii) but as a function of y, is obtained from r ε by the change of variables

y = x -εk ε , y 3 = x 3 ε , (4.39)
which transforms

Q k ε into Y ε .
Asymptotic behaviour of the rescaled functions ( w ε , q ε ) We first study in Lemmas 4.4.1 and 4.4.2 the asymptotic behaviour of ( w ε , q ε ). These results represent a generalization to the evolutive case of Lemmas 4.5 and 5.6 in [START_REF] Suárez-Grau | Effective boundary condition for a quasi-Newtonian viscous fluid at a slightly rough boundary starting from a Navier condition[END_REF], originally stated for the stationary case.

Lemma 4.4.1. Let π ε be a bounded sequence in L ∞ (0, T ; L 3/2 (Ω ε , R)) and π ε be the associated unfolded sequence. Then there exists π ∈ L ∞ (0,

T ; L 3/2 (T 2 × Q)) such that ∀M > 0 ε 2/3 π ε π weakly- * in L ∞ (0, T ; L 3/2 (T 2 × Q M )) . (4.40)
Moreover, defining q ε = -∂ t π ε and q = -∂ t π, there holds for every ϕ ∈

C 1 c ((0, T ) × T 2 , C 1 # ( Q, R 3
)), with ∇ y ϕ(x , y) = 0 a.e. in {y 3 > M }, for some M > 0:

lim ε→0 ε -1/3 (0,T )×Ω q ε (t, x) (div y ϕ)(t, x , x ε ) dx dt = (0,T )×T 2 × Q M q div y ϕ dx dy dt . (4.41)
Lemma 4.4.2. We consider a sequence of functions w ε ∈ J ε , where J ε is defined by (4.18). Assume that there exists a constant C > 0, such that

w ε H 1 (0,T ;L 2 (Ωε,R 3 )) ≤ C , w ε L 3 (0,T ;W 1,3 (Ωε,R 3 )) ≤ C ,
and that w ε|Ω converges weakly to a function w in L 3 (0, T ; W 1,3 (Ω, R 3 )).

Then, w 3 vanishes on (0, T ) × Γ b . Moreover, the following result holds.

(i) If a ε ≈ ε 5 3 , with a ε /ε 5 3 → λ, 0 < λ < +∞, then there exists w ∈ L 3 (0, T ; L 3 (T 2 , V 3 )) with w 3 (t, x , y , 0) = -λ∇Ψ(y ) • w (t, x , 0), on (0, T ) × T 2 × Y , (4.42)
such that for every M > 0, the sequence w ε satisfies

ε -2 3 ∇ y w ε ∇ y w in L 3 ((0, T ) × T 2 × Q M , R 3×3 ) . (4.43)
Besides, if div w ε = 0 in (0, T ) × Ω ε , then

div y w = 0 in (0, T ) × T 2 × Q . (4.44) (ii) If a ε ε 5 3 , then w (t, x , 0) • ∇Ψ(y ) = 0 on (0, T ) × T 2 × Y . (4.45)
Proof of Lemma 4.4.1 In order to give a flavour of the arguments used in the unfolding method, we provide the details of the proof of Lemma 4.4.1.

The proof of Lemma 4.4.2 is longer and more technical so we choose to omit it.

Let M > 0 and ε > 0. Using Remark 4.4.1, there holds

ε T 2 × Q M | π ε (t)| 3/2 dx dy = ε T 2 × Q M |π ε (t, εκ( x ε ) + εy , εy 3 )| 3/2 dx dy = ε 3 k ∈Kε Q M |π ε (t, εk + εy , εy 3 )| 3/2 dy = k ∈Kε C k ε ×(0,εM ) |π ε (t, x , x 3 )| 3/2 dx ≤ Ωε |π ε (t, x)| 3/2 dx ,
which gives (4.40).

Then, for a regular ϕ with ∇ y ϕ(t, x , y) = 0 a.e. in {y 3 > M }, the limit (4.41) comes from the following computation:

lim ε→0 ε -1/3 (0,T )×Ω q ε (t, x) (div y ϕ)(t, x , x ε ) dx dt = lim ε→0 ε -1/3 (0,T )×Ω π ε (t, x) (div y (∂ t ϕ))(t, x , x ε ) dx dt = lim ε→0 (0,T )×T 2 × Q M ε 2/3 π ε (div y (∂ t ϕ)) dx dy dt = (0,T )×T 2 × Q M π div y (∂ t ϕ) dx dy dt .
Proof of Theorem 4.3.3

First, thanks to (4.21), convergences (4.29) hold. The rest of the proof is carried out by using suitable test functions ϕ ε in the variational formulation (4.26). We give below the main steps of the demonstration.

Step 1. Obtaining the limit system We consider

χ ∈ C 1 c ([0, T ), R), ϕ ∈ C 1 c (Λ, R) 3 with ϕ 3 (x , 0) = 0, ϕ ∈ C 1 c (T 2 , C 1 # ( Q, R 3
)) with ∇ y ϕ(x , y) = 0 a.e. in {y 3 > M }, for some M > 0, such that

ϕ(x) = ϕ(x , 0) if x 3 ≤ 0, ϕ 3 (x , 0) = 0, ϕ(x , y) = ϕ(x , y , 0) if y 3 ≤ 0, ϕ 3 (x , y , 0) = -λ∇Ψ(y ) • ϕ (x , 0). (4.46) Now, take a truncature function ζ ∈ C ∞ (R) s. t. ζ(s) = 1 if s < 1/3 and ζ(s) = 0 if s > 2/3, and R ε > 0 such that R ε → ∞, R ε a ε ε 5 3 -λ 3 + a ε ε 3 → 0 .
Then, we define

ϕ ε ∈ W 1,3 (Ω ε , R 3 ) by                      ϕ ε (x) = ϕ (x) + ε 2 3 ϕ (x , x ε )ζ(x 3 ) ϕ ε,3 (x) = ϕ 3 (x) + ε 2 3 ϕ 3 (x , x ε )ζ(x 3 ) +ζ( x 3 εR ε )∇Ψ( x ε ) • λ - a ε ε 5 3 ϕ (x , 0) - a ε ε ϕ (x , x ε , 0) .
By (4.46) and the fact that ϕ(x) and ϕ(x , y) equal zero for x outside a compact subset of

T 2 , v ε (t, x) = χ(t) ϕ ε (x) is an admissible test function in (4.26).
Developing ∇ϕ ε as

∇ϕ ε (x) = ∇ϕ(x) + ε -1 3 ∇ y ϕ(x , x ε ) + h ε (x),
where

Ωε |h ε | 3 dx ≤ O ε + CR ε a ε ε 5 3 -λ 3 + a ε ε one obtains T 0 χ(t) Ωε ∂ t w ε (t) + (w ε (t) • ∇)w ε (t) • ϕ dx dt + T 0 χ(t) Ωε S(∇w ε (t)) : ∇ϕ + ε -1 3 ∇ y ϕ(x , x ε ) dx dt - T 0 χ(t) Ωε q ε div ϕ + ε -1 3 div y ϕ(x , x ε ) dx dt +β T 0 χ(t) Γε w ε • ϕ dH 2 dt = T 0 χ(t) Ωε f • ϕ dx dt + O ε .
(4.47)

Passing to the limit in (4.47). We focus on the terms involving ϕ, which are specific to the unfolding method.

Taking into account that

Ωε\Ω ∇ϕ + ε -1 3 ∇ y ϕ(x , x ε ) 3 dx ≤ C ε |Ω ε \ Ω| ≤ C a ε ε = O ε , we write T 0 χ(t) Ωε S(∇w ε (t)) : ∇ϕ + ε -1 3 ∇ y ϕ(x , x ε ) dx dt = T 0 χ(t) Ω S(∇w ε (t)) : ∇ϕ dx dt + ε -1 3 T 0 χ(t) Ω S(∇w ε (t)) : ∇ y ϕ(x , x ε ) dx dt + O ε . (4.48) Since S(∇w ε ) is bounded in L 3 2 ((0, T ) × Ω ε , R 3×3 ), there exists ξ such that S(∇w ε ) |Ω ξ in L 3 2 ((0, T ) × Ω, R 3×3 ) . (4.49)
By Lemma 4.4.2, S(ε

-2 3 ∇ y w ε ) is bounded in L 3 2 ((0, T ) × T 2 × Q M , R 3×3 ∀M > 0, so there exists ξ such that S(ε -2 3 ∇ y w ε ) ξ in L 3 2 ((0, T ) × T 2 × Q, R 3×3 ) . (4.50)
Hence, an easy computation gives

ε -1 3 T 0 χ(t) Ω S(∇w ε ) : ∇ y ϕ(x , x ε ) dx dt = T 0 χ(t) T 2 × Q ξ : ∇ y ϕ dx dy dt + O ε .
Thus, (4.49) and the above equality imply

T 0 χ(t) Ωε S(∇w ε ) : ∇ϕ + ε -1 3 ∇ y ϕ(x , x ε ) dx dt = T 0 χ(t) Ω ξ : ∇ϕ dx dt + T 0 χ(t) T 2 × Q ξ : ∇ y ϕ dx dy dt + O ε . (4.51)
Arguing similarly for the pressure term, we can passing to the limit in (4.47) and obtain by a density argument on ϕ:

T 0 Ω ∂ t w + (w • ∇)w • ϕ + ξ : ∇ϕ -q div ϕ + T 0 T 2 × Q ξ : ∇ y ϕ -q div y ϕ + β T 0 Γ b w • ϕ = T 0 Ω f • ϕ , (4.52)
for every ϕ ∈ L 3 (0, T ; W 1,3 (Ω, R 3 )) such that ϕ(t) = 0 on Γ t and ϕ 3 (t) = 0 on Γ b , and every ϕ ∈ L 3 ((0, T ) × T 2 , V 3 ) satisfying the compatibility condition

ϕ 3 (t, x , y , 0) = -λ∇Ψ(y ) • ϕ (t, x , 0) on (0, T ) × T 2 × Y .
It remains to establish that ξ = S(∇w) and ξ = S(∇ w).

Step 2. Identification of ξ and ξ Using the monotonicity of S and the convergence results (4.29), (4.43), (4.49), (4.50), we get the new monotonicity inequality:

T 0 Ω (ξ-S(∇η)) : (∇w-∇η)+ T 0 T 2 × Q ( ξ-S(∇ y η)) : (∇ y w-∇ y η) ≥ 0 . (4.53) Here, (η, η) is a regular admissible couple, in the sense that η ∈ C 1 c ([0, T ] × Λ, R 3 ) and η ∈ C 1 c ([0, T ] × T 2 , C 1 # ( Q, R 3
)), with ∇ y η(t, x , y) = 0 a.e. in {y 3 > M }, for some M > 0.

At this point, we rely on an idea by Minty [START_REF] Minty | Monotone (nonlinear) operators in Hilbert space[END_REF]. Take η = w -τ ζ for an arbitrary τ > 0, with

ζ ∈ C 1 c ([0, T ] × Λ, R 3
) and η = w in (4.53). This gives

T 0 Ω ξ -S(∇w -τ ∇ζ) : ∇ζ ≥ 0 .
Observing that, by monotonicity of S, S(∇w -τ ∇ζ) : ∇ζ ≤ S(∇w) : ∇ζ, one can let τ → 0 and use the arbitrariness of ζ to get the equality ξ = S(∇w). A similar argument gives ξ = S(∇ y w).

Step 3. Expression of ( w, q) Replacing ξ by S(∇ y w) and letting ϕ = 0 in (4.52), we see that ( w, q) is solution of

               -div y S(∇ y w) + ∇ y q = 0 in Q, div y w = 0 in Q, w 3 (t, x , y , 0) = -λ∇Ψ(y ) • w (t, x , 0) on T 2 × {0}, S(∇ y w) i,3 = 0, i = 1, 2, on T 2 × {0}, (4.54)
a.e. in (0, T ) × T 2 . Defining ( φ ξ , π ξ ) by (4.33), we deduce the expressions

w(t, x , y) = λ φ w (t,x ,0) (y), q(t, x , y) = λ 2 π w (t,x ,0) (y) . (4.55)
Conclusion of the proof Finally, fixing a test function ϕ, we define ϕ by ϕ(t, x , y) = λ 2 i=1 ϕ i (t, x , 0) φ e i (y) and conclude using (4.52) that (w, q) is the variational solution to problem (4.30)-(4.31).

Numerical results

In this section, we investigate numerically the behaviour of the solution of the LES model (4.11)-(4.13), in the case where the lower plate is composed of riblets, defined by

Ψ ε (x ) = a ε Ψ(x 1 /ε), x ∈ T 2 , (4.56)
where Ψ is a given, positive 1-periodic and regular function of variable x 1 , and the sequence (a ε ) ε>0 satisfies assumptions (4.1). By Corollary 4.3.5, the effective velocity and pressure w, q satisfy the limit system (4.30), completed with boundary conditions (4.36).

In order to stress the influence of parameter λ, defined by (4.27), we simulate the stationary system We stress that, in the case of riblets, according to boundary conditions (4.36), one needs only to compute the coefficient R 1 (ξ ) for ξ = e 1 and ξ = -e 1 , and to have access to the sign of w 1 (x , 0), to define the boundary condition at point x ∈ T 2 . Let us first address the numerical approximation of coefficients R 1 (e 1 ), R 1 (-e 1 ).

                 (w • ∇)w -div (S(∇w)) + ∇q = f in Ω, div w = 0 in Ω, w = 0 on Γ t , w 3 = 0 on Γ b , -S(∇w) 1,3 + β w 1 + λ 3 w 2 1 R 1 (sgn(w 1 )e 1 ) = 0 on Γ b , -S(∇w) 2,3 + β w 2 = 0 on Γ b .

Numerical methods and algorithms

Resolution of the auxiliary system (4.33). In order to compute coefficients R 1 (e 1 ), R 1 (-e 1 ), one needs to solve system (4.33), with ξ ∈ {e 1 , -e 1 }. We fix a truncature parameter H > 0 and we consider the following approximated problem: find (

φ ξ H , π ξ H ) ∈ W 1,3 (T 2 × (0, H), R 3 ) × L 3/2 (T 2 × (0, H)) satisfying                        -div y ( S(∇ y φ ξ H )) + ∇ y π ξ H = 0 in T 2 × (0, H), div y φ ξ H = 0 in T 2 × (0, H), φ ξ H,3 (y , 0) = -∇Ψ(y ) • ξ on T 2 × {0}, S(∇ y φ ξ H ) i,3 = 0, i = 1, 2, on T 2 × {0}, S(∇ y φ ξ H ) i,3 = 0, i = 1, 2, 3, on T 2 × {H}. (4.58)
Then, the approximate values of R 1 (ξ ), for ξ ∈ {e 1 , -e 1 }, will be given by

R 1 (ξ ) ≈ T 2 ×(0,H) S(∇ y φ ξ H ) : ∇ y φ e 1 H dy.
Such an approximation can be justified by standard arguments, relying on Proposition 4.3.4 that states the exponential decay of the gradient of φ ξ , along the y 3 -direction. Since problem (4.58) is nonlinear (due to tensor S), we solve it using a standard fixed point algorithm.

Resolution of the effective problem (4.57). We proceed in an analogous way as in the treatment of system (4.58), by introducing a fixed point algorithm. In particular, we treat the nonlinear term R 1 (sgn(w 1 )e 1 ) explicitely, in order to utilize the approximate values of R 1 (e 1 ) and R 1 (-e 1 ) previously determined.

Spatial discretization by finite element methods. We rely on mixed formulations of problems (4.57) and (4.58), discretized by finite element methods using FreeFem++ [START_REF] Hecht | New Development in FreeFem++[END_REF].

Numerical simulations

Definition of the parameters

In order to illustrate the effect of coefficient λ, we consider a flow driven by an external force f = e 1 . We set the Reynolds number to Re = 10000, a value for which the use of Large Eddy Simulation is justified. We choose δ as the grid size of the tetrahedralization of domain Ω, and define C s = δ 2 and γ = 6, which appears to be standard for the numerical simulation of Smagorinsky's turbulent model (see [START_REF] Volker | Large eddy simulation of turbulent incompressible flows[END_REF]). Such definition of δ is coherent with the modelling of the turbulence developed in Sec. 4.2, since δ is a small averaging radius and the grid size is the smallest length scale accessible in the mesh.

To set the value of parameter β, we follow an argument developed by Berselli et al. in [START_REF] Berselli | Mathematics of Large Eddy Simulation of Turbulent Flows[END_REF], based on Maxwell's theory of gases [START_REF] Maxwell | On stresses in rarified gases arising from inequalities of temperature[END_REF], and define β = L/(Re δ), where L is a characteristic length of the domain (in our case, we take L = 1). This definition is consistent with Remark 4.2.1.

Numerical results and comments

In order to establish comparisons between the solutions to the effective system (4.57), let us introduce certain quantities that are classically associated with the stationary flow of a viscous fluid, in the vicinity of a solid wall where friction occurs:

• the kinetic energy 1 2 Ω |w| 2 dx ;

• the dissipated energy, obtained by summing two contributions:

1. the viscous dissipation

1 2 Ω 1 Re |∇w| 2 + C s |∇w| 3 - δ 2 2γ ∇w(∇w) T : ∇w dx 2. the frictional dissipation 1 2 Γ b β |w| 2 + λ 3 w 3 1 R 1 (sgn(w 1
)e 1 ) dH 2 ;

• the average streamwise velocity Ω w 1 dx.

We have plotted in Fig. 4.3 the kinetic and dissipated energies, the viscous and frictional dissipation and the average streamwise velocity, against the roughness parameter λ.

For every riblet geometry, increasing λ leads to decreasing the kinetic energy of the fluid, as well as the dissipated energy and the streamwise velocity. Dissipation by friction is the only physical quantity presenting a nonmonotonic behaviour, reaching a peak before decreasing for large values of λ. However, viscous dissipation dominates the frictional one, so that the dissipated energy behaves similarly to the viscous one with respect to λ.

Conclusion and perspectives

In this study, we have presented a model of practical interest in the simulation of turbulent flows, and identified an effective model accounting for the effect of wall roughness on the characteristics of the flow. We have also shown by numerical simulations that this effective model was usable in practice, in the case of a ribbed geometry which is widely studied in the experimental litterature on turbulent fluid mechanics. From a technical perspective, the asymptotic analysis that we have proposed relies on a careful use of the monotonicity properties of the operator associated with the Smagorinsky approximation of the Newtonian stress tensor.

Even though Smagorinsky model is still one of the most popular LES models, it possesses some well-known and documented drawbacks. For instance, the constant C s is an a priori input, which prevents the model to represent correctly various turbulent flows. Also, Smagorinky model introduces an excessive diffusion in the flow.

In the light of the results obtained in this chapter, we intend to extend our approach to more refined and up to date LES models, still enjoying nice mathematical properties but relying on more sophisticated approximations of the convolution by the Gaussian filter. One promising lead is the so-called Rational LES model proposed by Galdi and Layton in [START_REF] Galdi | Approximation of the Larger Eddies in Fluid Motions II: a Model for Space-Filtered Flow[END_REF].

Introduction

Microfluidics deals with the the experimental and theoretical study of lubricants that are geometrically constrained to a sub-millimeter scale. As previously introduced in Chapter 3, their behaviour can differ from "macrofluidic" behaviour since, at the microscale, factors such as surface tension, energy dissipation, and fluidic resistance start to dominate the system. As a result, one needs to understand and quantify very precisely the behaviour of the fluid near a solid wall in this context.

The mathematical models of lubrication usually result from the simplification of the lubricant film geometry, i.e. its thickness. Using this thickness as a small parameter ε, a simple asymptotic approximation can be easily derived providing a well-known Reynolds equation for the pressure of the fluid. Formal derivation goes back to the 19th century and the celebrated work of Reynolds [START_REF] Reynolds | On the theory of lubrication and its application to Mr. Beauchamp tower's experiments, including an experimental determination of the viscosity of olive oil[END_REF]. We refer to [START_REF] Bayada | The transition between the Stokes equations and the Reynolds equation: a mathematical proof[END_REF] for a rigorous justification of this approximation in the case of a Newtonian flow.

Nevertheless, most of the modern lubricants are no longer Newtonian fluids, since the use of additives has become a common practice to improve their performance. Therefore, several microcontinuum theories have been proposed to account for the effects of additives. In this chapter, we rely on Eringen micropolar fluid theory [START_REF] Eringen | Theory of micropolar fluids[END_REF]. In this model, every material volume contains undeformable microelements that can translate and rotate independently on the macrovolume motion. This leads to the introduction of a microrotation field w in the equations of motion of the fluid.

In this framework, a micropolar Reynolds equation was obtained in [START_REF] Bayada | On micropolar fluids in the theory of lubrication. Rigorous derivation of an analogue of the Reynolds equation[END_REF] for a micropolar flow in a thin film with a plain bottom assuming zero boundary conditions for microrotation. Other related results can be found in [START_REF] Prakash | Lubrication theory for micropolar fluids and its application to a journal bearing[END_REF]. However, more general boundary conditions for the microrotation were introduced to take into account the rotation of the microelements on the solid boundary. In the case where the boundary is flat, these conditions read α 2 (∇ × u) × ν = w × ν .

(5.1) Conditions (5.1) were effectively proved to be in good accordance with experiments [START_REF] Bessonov | Boundary viscosity conception in hydrodynamical theory of lubrication[END_REF][START_REF] Bessonov | A new generalization of the reynolds equation for a micropolar fluid and its application to bearing theory[END_REF]. The coefficient α describes the interaction between the given fluid and solid; it characterizes microrotation retardation on the solid surfaces. Nevertheless, it was proved impossible to simultaneously impose the boundary condition (5.1) and a no-slip condition for the velocity [START_REF] Bayada | New models in micropolar fluid and their application to lubrication[END_REF]. This would be like considering simultaneously, at the same boundary, a Neumann and a Dirichlet boundary condition. Hence, a compatible condition for the velocity must be introduced, allowing for slippage in the tangential direction, of the form (u -s)

× ν = δ 0 (∇ × w) × ν , u • ν = 0 (5.2)
where δ 0 is a parameter and s is the velocity of the wall. This condition appears as a new interpretation of the slippage observed in lubrication with micropolar fluids, expressed in terms of the microrotation field w.

In [START_REF] Bayada | New models in micropolar fluid and their application to lubrication[END_REF], using the nonzero boundary conditions (5.1)-(5.2), Bayada et al. derive rigourously a generalized version of the Reynolds equation in a 2d domain with flat boundary. In this chapter, we address the combined effect of such nonstandard boundary conditions with a rough wall geometry, in the particular case of riblets. Using the unfolding method (previously described in Chapter 4), we identify a critical scaling between the amplitude and the period of the roughness, leading to three different asymptotic regimes in the lubrication limit. In each case, we provide the corresponding micropolar Reynolds equation. Finally, we apply these results to carry out a numerical study of a model of squeeze-film bearing lubricated with a micropolar fluid.

A model for micropolar fluid lubrication

The fluid domain is given by

Ω ε = (x , x 3 ) ∈ ω × R, -Ψ ε x < x 3 < h ε , (5.3)
where ω is a smooth domain in R 2 , h > 0 is a dimensionless constant and Ψ ε is given by

Ψ ε (x ) = λ ε δ Ψ 1 ε x • e 1 (5.4)
for a regular, 1-periodic function Ψ that models the roughness profile on the lower surface. The lower, upper and lateral boundaries are respectively described by

Γ 0 ε = (x , -Ψ ε x ), x ∈ ω , Γ 1 ε = (x , h ε) x ∈ ω , Γ ε = ∂Ω ε -(Γ 0 ε ∪ Γ 1 ε ) .
We denote by u ε the velocity of the fluid, p ε the pressure and w ε the microrotation of the solid particles. Let µ be the viscosity of the fluid and µ r , c r be the micropolar viscosities. Setting N 2 = µr µ+µr and R M = cr µ+µr , the micropolar fluid flow is described by the following equations expressing the balance of momentum, mass and angular momentum:

-∆u ε + ∇p ε = 2N 2 (∇ × w ε ) in Ω ε , (5.5
)

div u ε = 0 in Ω ε , (5.6
)

-R M ∆w ε + 4N 2 w ε = 2N 2 (∇ × u ε ) in Ω ε . (5.7)
We consider the case of a squeeze-film motion: the upper plate Γ 1 ε moves downwards with a given velocity -V ε e 3 , while the lower plate Γ 0 ε is at rest. We denote by g ε the fluid velocity on the lateral boundaries of the domain. As discussed in the Introduction, we impose the following boundary condi-tions:

u ε = -V ε e 3 , w ε = 0 on Γ 1 ε , (5.8
)

u ε = g ε , w ε = 0 on Γ ε , (5.9
)

u ε • ν = 0, w ε • ν = 0 on Γ 0 ε , (5.10) α 2 [∇u ε ν ε ] tan = w ε × ν on Γ 0 ε , (5.11)
R M [∇w ε ν] tan = 2N 2 β u ε × ν on Γ 0 ε .
(5.12)

The incompressibility condition (5.6) imposes the following compatibility condition on the boundary data:

Γ ε g ε • ν dH 2 = V ε |ω| , (5.13)
where |ω| is the area of ω.

Let us point that condition (5.11) is an adaptation of (5.1) to the present case of a non flat boundary Γ 0 ε . This is justified by the fact that, for a flat boundary Γ and regular vector fields u, v satisfying u

• ν = v • ν = 0 on Γ, Γ [(∇ × u) × ν] • v dH 2 = Γ [∇u ν] tan • v dH 2 .
Interpretation of the model In [START_REF] Bessonov | Boundary viscosity conception in hydrodynamical theory of lubrication[END_REF], Bessonov proposed to interpret parameter α appearing in (5.11) as a microrotation retardation at the boundary and to connect it with the different viscosity coefficients, by a formula of the form

α = ν + ν r -ν b ν r (5.14)
where ν b is a boundary viscosity that differs from ν and ν r . Following [START_REF] Bessonov | Boundary viscosity conception in hydrodynamical theory of lubrication[END_REF], it is possible to give physical limits to ν b , inducing limits on α:

0 ≤ ν b ≤ ν + ν r ⇒ 0 ≤ α ≤ ν + ν r ν r . (5.15)
Condition α = 0 is equivalent to strong adhesion of the microparticles on the boundary surface, i.e. no rotation: w = 0. Thus, from now on, we consider α > 0 so that the stress tensor and the micro-rotation are coupled on the boundary. Finally, coefficient β ∈ R in (5.12) is a friction coefficient that controls the slippage of the fluid at the wall.

Weak formulation of problem (5.5)-(5.12) We define the functional spaces

V ε = {ϕ ∈ H 1 (Ω ε , R 3 ), ϕ |Γ 1 ε ∪Γ ε = 0, ϕ • ν = 0 on Γ 0 ε }, V 0 ε = {ϕ ∈ V ε , div ϕ = 0 in Ω ε } . Definition 5.2.1. We say that (u ε , w ε , p ε ) ∈ H 1 (Ω ε , R 3 ) 2 × L 2 (Ω ε )/R
(ϕ, ψ) ∈ V ε × V ε : Ωε ∇u ε : ∇ϕ dx - Ωε p ε div ϕ dx + R M Ωε ∇w ε : ∇ψ dx -2N 2 Ωε (∇ × u ε ) • ψ dx -2N 2 Ωε w ε • (∇ × ϕ) dx + 4N 2 Ωε w ε • ψ dx -2 1 α -N 2 Γ 0 ε (w ε × ν) • ϕ dH 2 -2N 2 β Γ 0 ε (u ε × ν) • ψ dH 2 = 0.
(5.16)

Main results

Asymptotic regimes on R M and V ε Different asymptotic behaviours of the flow may be deduced depending on the order of magnitude of the viscosity coefficients. We introduce a constant R c > 0 and consider here the asymptotic regime

R M = ε 2 R c , (5.17)
that leads to a strong coupling at main order of the velocity and microrotation.

Besides, we assume that the (vertical) velocity of the upper plate V ε is of order ε and we set

V ε = εS, (5.18)
where S is a positive constant.

Regularity of boundary data g ε We also need certain regularity on the boundary data g ε , and uniform estimates of relevant norms. A very general way of those properties is the following: there exists a sequence of lift functions

J ε ∈ H 1 (Ω ε , R 3 ) satisfying div J ε = 0 in Ω ε , the boundary conditions J ε = -V ε e 3 on Γ 1 ε , J ε = g ε on Γ ε , J ε • ν = 0 on Γ 0 ε , (5.19)
and the estimates

J ε L 2 (Ωε,R 3 ) ≤ Cε 1 2 , ∇J ε L 2 (Ωε,R 3×3 ) ≤ Cε -1 2 , J ε L 2 (Γ 0 ε ,R 3 ) ≤ C, (5.20
) where C > 0 is a universal constant.

Existence and uniqueness result Theorem 5.3.1. Assume that the coupling parameter N 2 satisfies the condition

N 2 ≤ 1 2 , (5.21)
and define the nonnegative parameter γ by γ = 1 α -N 2 -N 2 β. Assume that the asymptotic regimes (5.17) and (5.18) hold. Then, for any β such that

γ 2 < R c (1 -2N 2 ) h 2 , ( 5.22)
there exists ε 0 > 0 such that for any 0 < ε < ε 0 , system (5.5)-( 5.12) has a unique weak solution.

Effective system and Reynolds equation

To describe the asymptotic behaviour of the sequence (u ε , w ε , p ε ), we introduce a change of variables classically used for thin domains: the dilatation

y = x , y 3 = x 3 ε , (5.23)
which transforms Ω ε in the set Ω ε of height of order h, defined as follows:

Ω ε = (y , y 3 ) ∈ R 2 × R : y ∈ ω, -Ψ ε (y ) < y 3 < h , (5.24)
where Ψ ε (y ) = ε δ-1 Ψ 1 ε y • e 1 . Since δ > 1, it is clear that the sequence of domains Ω ε converges in usual senses to the limit domain Ω = ω × (0, h), of lower boundary Γ := ω × {0}. Accordingly, we define the rescaled functions

u ε , w ε ∈ H 1 ( Ω ε , R 3 ) and p ε ∈ L 2 ( Ω ε )/R by u ε (y) = u ε (y , εy 3 ), w ε (y) = w ε (y , εy 3 ), p ε (y) = p ε (y , εy 3 ). (5.25)
Depending on the relation between the amplitude parameter δ and the period parameter , we obtain three different asymptotic regimes, that we call critical, subcritical and supercritical. To lighten the presentation we only give the result in the critical case.

Generalized Reynolds equation

For simplicity, we consider here that the flow is independent on the y 2 coordinate, and that velocity and microrotation components u 2 and w 1 are zero. Thus, we address the following limit problem posed in Ω = (0, 1) × (0, h):

-∂ 2 y 3 u 1 + ∂ y 1 p + 2N 2 ∂ y 3 w 2 = 0 in Ω, -R c ∂ 2 y 3 w 2 + 4N 2 w 2 -2N 2 ∂ y 3 u 1 = 0 in Ω, (5.30)
with the boundary conditions

u 1 = 0, w 2 = 0 on (0, 1) × {h} , ∂ y 3 u 1 = 2 α w 2 + λ 2 E u 1 , on (0, 1) × {0} , R c ∂ y 3 w 2 = -2N 2 β u 1 on (0, 1) × {0} ,
and the incompressibility condition Sθ(y 1 ) dy 1 , ∀θ ∈ H 1 (0, 1), (5.32) where coefficient Θ λ depends explicitely on λ and other parameters N 2 , R c , α, E and h.

∂ y 1 h 0 u 1 (y 1 , y 3 ) dy 3 = S in (0, 1) . ( 5
Remark 5.3.1. We do not reproduce in this chapter the expression of Θ λ , which is rather long and complicated, but refer to [B3] for details.

Elements of proof 5.4.1 Existence and uniqueness of the solution

Thanks to the regularity assumption on the boundary data g ε , one can introduce a lift function J ε satisfying (5.19) and set

u ε = v ε + J ε . Prob- lem (5.5)-(5.12) is thus equivalent to finding a solution (v ε , w ε , p ε ) ∈ V 0 ε × V ε × L 2 (Ω ε )/R to the mixed formulation A ε ((v ε , w ε ), (ϕ, ψ)) + B ε ((ϕ, ψ), p ε ) = L ε (ϕ, ψ) ∀(ϕ, ψ) ∈ V ε × V ε , (5.33) B ε ((v ε , w ε ), q) = 0 ∀q ∈ L 2 (Ω ε ) , (5.34)
where

A ε ((v, w), (ϕ, ψ)) = Ωε ∇v : ∇ϕ dx - Ωε p ε div ϕ dx + ε 2 R c Ωε ∇w : ∇ψ dx -2N 2 Ωε (∇ × v) • ψ dx -2N 2 Ωε w • (∇ × ϕ) dx + 4N 2 Ωε w • ψ dx -2 1 α -N 2 Γ 0 ε (w × ν) • ϕ dH 2 -2N 2 β Γ 0 ε (v × ν) • ψ dH 2 , B ε ((v, w), q) = - Ωε q div v dx , L ε (ϕ, ψ) = - Ωε ∇J ε : ∇ϕ dx + 2N 2 Ωε (∇ × J ε ) • ψ dx + 2N 2 β Γ 0 ε (J ε × ν) • ψ dH 2 .
The existence and uniqueness of such (v ε , w ε , p ε ) is established in [START_REF] Bayada | New models in micropolar fluid and their application to lubrication[END_REF] for a flat lower boundary, and holds in that case for any N 2 ≤ 1. Our particular condition (5.21) results from the introduction of an oscillating boundary Γ 0 ε . Indeed, due to the presence of divergence and rotational terms, studying the coerciveness of A ε requires to relate the L 2 norm of the gradient of a given vector field, to the L 2 norm of its divergence and curl. However, we cannot rely on the well-known formula

Ωε |divϕ| 2 + |∇ × ϕ| 2 dx = Ωε |∇ϕ| 2 dx + Γ 0 ε ((ϕ • ∇)ν) • ϕ dH 2 (5.35)
(see, for instance, [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF], formula (IV.23)) since in our geometrical setting, ∇ν if of order ε δ-2 , hence diverging since ε δ-tends to zero. Instead, we use the elementary estimate

∇ × ϕ L 2 (Ωε) 3 ≤ √ 2 ∇ϕ L 2 (Ωε,R 3×3 ) ∀ϕ ∈ H 1 (Ω ε , R 3 ) ,
where the √ 2 constant explains the upper bound 1/2 in condition (5.21).

A priori estimates and convergences

Proposition 5.4.1. Assume that the asymptotic regimes (5.17) and (5.18) and conditions (5.21) and (5.22) hold. Then (u ε , w ε , p ε ) satisfies the estimates

u ε L 2 (Ωε,R 3 ) ≤ Cε 1 2 , ∇u ε L 2 (Ωε,R 3×3 ) ≤ Cε -1 2 , (5.36) w ε L 2 (Ωε,R 3 ) ≤ Cε -1 2 , ∇w ε L 2 (Ωε,R 3×3 ) ≤ Cε -3 2 , (5.37) p ε L 2 (Ωε) ≤ Cε -3 2 , ∇p ε H -1 (Ωε,R 3 ) ≤ Cε -1 2 .
(5.38)

with p 0 ε ∈ H 1 (ω) (it does not depend on the variable x 3 ), p 1 ε ∈ L 2 (Ω ε ), and

p 0 ε H 1 (ω) ≤ Cε -3 2 ∇p ε H -1 (Ωε) 3 , p 1 ε L 2 (Ωε) ≤ C ∇p ε H -1 (Ωε,R 3 ) . (5.46)
We denote by p 1 ε the rescaled function associated with p 1 ε , defined by p 1 ε (y) = p 1 ε (y , εy 3 ) for y ∈ Ω ε .

Corollary 5.4.4. There exist p ∈ H 1 (ω) and p 1 ∈ L 2 (Ω) such that

ε 2 p 0 ε p in H 1 (ω), ε p 1 ε|Ω p 1 in L 2 (Ω), (5.47)
and moreover ε 2 p ε|Ω → p in L 2 (Ω).

(5.48)

Proof. From (5.46) and (5.38) 2 , we get

p 0 ε H 1 (ω) ≤ Cε -2 , p 1 ε L 2 (Ωε) ≤ Cε -1 2 , (5.49)
and after rescaling p 1 ε , last inequality becomes p 1 ε L 2 ( Ωε) ≤ Cε -1 . Previous estimates and the fact that Ω ⊂ Ω ε and | Ω ε \Ω| → 0 imply (5.47). The strong convergence (5.48) for the complete pressure p ε is a direct consequence of (5.47) and the decomposition (5.45).

Effective system and Reynolds equation

The proof of Theorem (5.3.2) is based on the unfolding method. Since we have already addressed this method in Chapter 3, we skip this part for brevity.

The Reynolds equation is then obtained by explicitely solving the effective problem (5.30), and by plugging the corresponding formula for u 1 in (5.31).

Application to squeeze-film bearing

We consider in this section a squeeze-film bearing composed of two parallel plates separated by a micropolar fluid film. The lower surface is at rest and composed of a rough material, while the smooth upper surface is under normal squeeze motion. Hence, the distance between the plates is a decreasing function of time, which is expected to go to zero as the fluid is squeezed out of the gap. We assume that the motion of the upper plate is slow, so that the inertial effects can be neglected and the behaviour of the bearing can be captured by a quasi-static model.

Influence of parameters N, R c , νb , δ, E

To perform comparisons, we have chosen to unify the presentation of the numerical results by plotting the half-life time T half against N ∈ [0, 0.7] (which ensures that condition N 2 ≤ 1/2 is fulfilled), after normalization by its value for N = 0, for different values of R c ∈ {0.025, 0.05, 0.1, 0.2} and using various sets of parameters νb , δ, E.

Influence of νb , N, R c Fig. 5.1 suggest that parameter νb does not have much of an impact on the computed value of T half . Consequently, we have decided to fix it to νb = 0.1 in the rest of the simulations. On the opposite, T half is very sensitive to the couple of parameters (N, R c ).

Influence of E Fig. 5.2 shows that increasing the value of E produces a visible change in the model: the function N → T half becomes monotonic for the highest tested values of R c . This means that for certain micropolar fluids, the roughness of the lower plate may contribute to enhance the performance of the bearing, in the sense that T half becomes larger.

Influence of δ We observe on Fig. 5.3 that small values of δ favor an enhancement of the performance of the bearing, while large values of δ typically have a negative impact.

Conclusion and perspectives

In this chapter, we have derived a Reynolds equation describing the lubrication limit of a micropolar fluid model, taking into account recently introduced boundary conditions that finely describe the interaction between the lubricant and the surface of a bearing, in the context of a ribbed wall geometry. We have applied our results on the numerical study of a model of squeeze-film motion, bringing out the complex behaviour of the system with respect to physical characteristics of the fluid and geometric properties of the wall.

A comparable approach, combining asymptotic analysis and numerical experiments on the effective model, could be applied to other types of models of non-Newtonian fluids with practical use in modern lubrication. A first step in this direction is an ongoing work with M. Anguiano and F. J. Suárez-Grau, on the homogenization of a model of quasi-Newtonian flow through a porous medium. Chapter 6

Modeling and shape-optimization of aquaporins

This chapter presents results obtained in collaboration with F. Omnès and Y. Privat in [B8], on the numerical optimization via shape optimization techniques of a model of aquaporin proposed by soft matter physicists in [START_REF] Gravelle | Large permeabilities of hourglass nanopores: From hydrodynamics to single file transport[END_REF].

Introduction and modeling of the problem 6.1.1 Motivations

Aquaporins are proteins found in cell membrane in plants, bacteria and several organs of animals including humans. Along with membrane diffusion, aquaporins realize water filtration through the lipid bilayer, regulating the passage of water, ions and other solutes. From a technological perspective, finding artificial systems enjoying high energy performance has applications in domains where membrane filtration is involved, from water desalination, to industrial food processes and wastewater treatment.

Notwithstanding, numerical simulation of membranes is expected to help imagining new efficient structures. Since the pioneering work of Hummer et al. [START_REF] Hummer | Water conduction through the hydrophobic channel of a carbon nanotube[END_REF], many simulations of fluid transport at a nanometric scale have been developed. These simulations can be split up into two categories: • molecular dynamics (MD) models, that consider molecular interactions [START_REF] Karniadakis | Simple Fluids in Nanochannels[END_REF][START_REF] Park | Carbon nanofluidics of rapid water transport for energy applications[END_REF] ;

• methods based on continuum mechanics, that rely on the numerical resolution of the Navier-Stokes equations, using finite element methods (FEM) [START_REF] Gravelle | Large permeabilities of hourglass nanopores: From hydrodynamics to single file transport[END_REF] or spectral elements methods [START_REF] Karniadakis | Simple Fluids in Nanochannels[END_REF]Chapter 14].

To our best knowledge, MD have not yet been used to deal with shapeoptimization problems. However, testing on a wide range of parameters, MD and FEM have been shown to provide highly similar results for hydrodynamic resistance [START_REF] Gravelle | Large permeabilities of hourglass nanopores: From hydrodynamics to single file transport[END_REF]. This justifies the use of FEM-based methods for simulating the transport of fluids by aquaporins.

In [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF], Gravelle et al. investigate an hourglass model of aquaporin, imposing partial-slip boundary conditions on the channel walls. Varying the angle of the inlet and outlet cones, they show numerically the existence of an optimal angle minimizing the energy dissipation inside the water channel. The class of admissible shapes is extended in [START_REF] Belin | Optimal shape of entrances for a frictionless nanochannel[END_REF], where the shape of the inlet cone is optimized through three design parameters.

In this chapter, we present a generalization of the results in [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF][START_REF] Belin | Optimal shape of entrances for a frictionless nanochannel[END_REF], based on shape optimization techniques. In Section 6.1.2, we introduce the fluid model for aquaporins and gather some tools for the analysis of the resulting system of PDEs. In Section 6.1.3, the shape optimization problem aiming at improving the performances of aquaporins is introduced. We then analyze this problem in Section 6.2, investigating existence issues and determining a workable expression of the cost function shape derivative. Finally, Figure 6.2: The domain Ω (two reservoirs connected by an hourglass shaped channel). Section 6.3 is devoted to numerical issues, and describes the main features of the algorithm that we have developed to address the shape optimization problem in practice.

Geometry and fluid model

In what follows, we consider connected and bounded domains Ω in R 2 , describing possible geometries of aquaporins (see Fig. 6.2). The fluid domain Ω is made up of two reservoirs delimited by a lateral boundary Γ 0 , connected by a channel, and filled with a viscous fluid of viscosity µ. The central part of the channel is tubular and its lateral boundary is denoted Γ 2 . The inlet and outlet regions of the channel are conic, with lateral boundary Γ 1 . The upstream and downstream sections are resp. labeled Γ in and Γ out .

Fluid model and boundary conditions

Let u be the fluid velocity, p the pressure and σ(u, p) be the stress tensor σ(u, p) = 2µD(u) -p I , where I is the identity matrix in R 2×2 and µ > 0 is a constant viscosity coefficient. The fluid motion is described by Stokes equations -div(σ(u, p)) = 0 in Ω , div u = 0 in Ω . (6.1a) (6.1b)

The momentum (6.1a) and mass-conservation (6.1b) equations are completed with boundary conditions.

• On the upstream section Γ in , we fix a flow rate Q > 0 and impose

Γ in u • ν dH 1 = -Q . (6.2)
Following [START_REF] Formaggia | Numerical treatment of defective boundary conditions for the Navier-Stokes equations[END_REF], we interpret condition (6.2) as a (linear) constraint on the unknown u defined as a minimizer of an energy functional.

• In realistic applications, the conical regions of the channel are of nanometric size. At this scale, partial slip boundary conditions are considered relevant (see [START_REF] Bocquet | Flow boundary conditions from nanoto micro-scales[END_REF] and the Introduction to Chapter 3). Consequently, we set

u • ν = 0, [σ(u, p)ν] tan + βu tan = 0 on Γ 1 , (6.3)
where β > 0 is a constant friction parameter.

• Since the two reservoirs that are connected by the nanopore are typically of much larger size than the nanometric central channel, we impose the classical no-slip condition on their boundaries, that is, u = 0 on Γ 0 . (6.4)

• On the lateral boundary Γ 2 of the central tube, we neglect the frictional dissipation by assuming perfect slip boundary conditions u • ν = 0, [σ(u, p)ν] tan = 0 on Γ 2 . (6.5)

• Finally, on the downstream section Γ out of the domain, we consider free outflow conditions, modeled by Neumann boundary conditions σ(u, p)ν = 0 on Γ out . (6.6)

These boundary conditions are commonly used as passive conditions on artificial boundaries [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible navier-stokes equations[END_REF].

Mixed formulation of the Stokes problem with imposed inner flow through Γ in .

Let us introduce the functional space Hence, in the symmetric case, the flow is fully described by its restriction to the upper part Ω + of the domain, using the symmetry conditions 6.11.

V(Ω) = ϕ ∈ H 1 (Ω, R 2 ), ϕ |Γ 0 = 0, (ϕ • ν) |Γ 1 ∪Γ 2 = 0 .

The shape optimization problem

From a physical point of view, it is reasonable to look for a shape minimizing the energy dissipated by the fluid inside the aquaporin. Indeed, this criterion measures the irreversible conversion of mechanical energy into internal energy or heat.

The cost functional is then defined by

J(Ω) = 2µ Ω |D(u Ω,λ )| 2 dx + β Γ 1
|u Ω,λ | 2 dH 1 , (6.12)

where the triple (u Ω,λ , p Ω,λ , λ) ≡ (u λ , p λ , λ) ∈ V(Ω) × L 2 (Ω) × R is defined in Proposition 6.1.1.

Since our main objective is to improve the results obtained in [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF], we propose the following optimization strategy.

• First, we determine numerically the optimal angle of the inlet and outlet cones, minimizing the energy dissipation J(Ω) among all the domains Ω whose geometry is described by Figure 6.2. This step is a reproduction of the analysis performed in [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF], and leads to obtaining a reference domain Ω * , characterized by the position of the conical part Γ * 1 of its boundary.

• Our main contribution consists in developing a method for optimizing the shape of the conical regions of the domain, starting from the domain Ω * that was determined in the previous step. In this process, the only part of the boundary that is deformed is Γ * 1 , while the rest of the boundary of Ω * is fixed.

Analysis of the shape optimization problem

This section is devoted first to the statement of an existence result for the shape optimization problem (6.14), and second, to the writing of the first order necessary optimality conditions for this problem.

Existence issues

It can be noted that the admissible class O ad defined by (6.13) is not closed for usual domains topologies such as Hausdorff complementary topology or strong L 1 convergence of characteristic functions. To avoid the emergence of irregular shapes, we choose to impose geometrical constraints on the free boundary Γ 1 .

First, one wants to deal with (at least) Lipschitz domains since the definition of the functional space V(Ω) involves the outward pointing normal vector, and since such regularity is required for using standard tools in the analysis of variational problems in Fluids Mechanics. However, a minimizing sequence of Lipschitz domains may converge to a very irregular domain [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Henrot | Variation et optimisation de formes[END_REF].

A possible solution consists in restricting the class of admissible domains, by assuming some kind of uniform Lipschitz regularity. For that purpose, we recall the notion of ε-cone property, introduced in [START_REF] Chenais | On the existence of a solution in a domain identification problem[END_REF]. Definition 6.2.1. Let y be a point of R 2 , ξ a normalized vector and ε > 0. We denote by C(y, ξ, ε), the unpointed cone C(y, ξ, ε) = {z ∈ R 2 , z -y, ξ ≥ cos ε z -y and 0 < z -y < ε}.

We say that an open set Ω verifies the ε-cone property if ∀x ∈ ∂Ω, ∃ξ x ∈ S 1 , ∀y ∈ Ω ∩ B(x, ε), C(y, ξ x , ε) ⊂ Ω. Another geometrical constraint, which is standard in shape optimization, is to assume that all admissible shapes are contained in a compact set D to avoid the degeneracy of the free boundary. For this reason, we introduce an rectangular external box D, defined as the convex hull of Ω.

Finally, the shape optimization problem that we investigate reads inf{J(Ω), Ω ∈ O ad , Ω ⊂ D and Ω satisfies the ε-cone property} (6.15)

for some given parameter ε > 0.

One has the following existence result.

Theorem 6.2.1. The shape optimization problem (6.15) has a solution.

Computation of the shape derivative

We need to differentiate the solution u λ ∈ V(Ω) to system (6.7a)-(6.7c), with respect to deformations of the domain Ω acting on Γ 1 , but preserving the other parts of the boundary of Ω.

Let V be a regular vector field with compact support in R 2 , vanishing on ∂Ω \ Γ 1 . We introduce T > 0 and a regular mapping t ∈ (-T, T ) → Φ t = (Φ 1 t , Φ 2 t ), satisfying Φ 0 = Id and dΦt dt |t=0 = V. Moreover, we assume that Φ t (x) = x for every x ∈ ∂Ω \ Γ 1 and t ∈ (-T, T ). A typical choice is given by Φ t = Id + tV.

We define Ω t := Φ t (Ω) and denote by u λ,t ∈ V(Ω t ) the solution of system (6.7a)-(6.7b)-(6.7c) for Ω = Ω t . One has the following differentiability result. Proposition 6.2.2. Let Ω ∈ O ad . The mapping

t ∈ (-T, T ) → (u λ,t • Φ t , p λ,t • Φ t ) ∈ H 1 (Ω, R 2 ) × L 2 (Ω) is differentiable at t = 0.
In what follows, we will denote by dJ(Ω), V the shape derivative of J at Ω in the direction V, in other words dJ(Ω), V = lim t 0 J(Ω t ) -J(Ω) t .

From now on, we will assume more regularity on the domain Ω and on the vector field V in order to get a workable expression of the shape derivative, namely that ∂Ω is of class C 2 and V ∈ W 3,∞ (R 2 , R 2 ). These properties ensure that the boundary of the domain Ω t remains of class C 2 , provided that t is small enough (see e.g. [START_REF] Delfour | Shapes and geometries[END_REF]). Theorem 6.2.3. Assume that ∂Ω is C 2 . Let the triple (u λ , p λ , λ) be the unique solution to (6.7a)-(6.7b)-(6.7c). For every V ∈ W 3,∞ (R 2 , R 2 ) having a compact support that does not intersect ∂Ω\Γ 1 , there holds

dJ(Ω), V = Γ 1 j 1 (V • ν) + j 2 • ∇ Γ (V • ν) dH 1 (6.16)
where ∇ Γ stands for the tangential gradient on Γ 1 , and

j 1 = 2µ|D(u λ )| 2 + β ∂ n (|u λ | 2 ) + H|u λ | 2 -2(σ(u λ , p λ )ν • ν)∂ n (u λ • ν) , j 2 = 2(σ(u λ , p λ )ν • ν)[u λ ] tan .
In the definition of j 1 , H denotes the mean curvature on Γ 1 .

Numerical methods and algorithms

Relying on the tools developed in Section 6.2, we have developed an efficient numerical method for solving problem (6.14). The implementation of the method was performed using FreeFem++ finite element software [START_REF] Hecht | New Development in FreeFem++[END_REF].

In an attempt to improve the numerical results in [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF][START_REF] Belin | Optimal shape of entrances for a frictionless nanochannel[END_REF], we have enriched their approach by

• considering a wider class of admissible shapes for Γ 1 ,

• using numerical shape-optimization techniques based on the computation of the shape derivative.

We have restricted the shape optimization problem (6.14) to symmetric shapes, as stated in the last paragraph of Subsection 6.1. Our approach can be decomposed into two main steps:

Step 1. Following [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF], we recover the optimal inner angle between Γ 1 and Γ 2 .

Step 2. Starting from the resulting straight cone with optimal angle, we find a local minimizer for the shape optimization problem (6.14), restricted to symmetric domains.

Choice of parameters

In order to perform numerical tests, reasonable physical parameters µ, β, Q as well as geometric dimensions for the aquaporin are required. In what follows, we will concentrate on two relevant test-cases to present our numerical results (see Table 6.1).

Finding the optimal angle

As a validation of our model, we start by reproducing the analysis in [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF] and solve the optimization problem inf{J(Ω θ ), θ ∈ [0, θ max )} ,

where Ω θ stands for the domain Ω with angle θ between Γ 1 and Γ 2 . To this aim, we use a simple dichotomy procedure on the angle θ, obtaining θ * = 0.265 ± 0.001 for test-case 1 and θ * = 0.311 ± 0.001 for test-case 2 (see Fig. 6.4).

Optimizing the shape of Γ 1

In our numerical simulations, we consider polygonal shapes Ω, symmetric with respect to the hyperplane H = {x 2 = 0}, and build a conforming simplicial mesh T of their upper part Ω + = Ω ∩ {x 2 > 0} [START_REF] Frey | Mesh generation: application to finite elements[END_REF]. Let us now provide the skeleton of the algorithm. 6.1).

Shape optimization algorithm

• Initialization: Choose an initial admissible domain Ω, symmetric with respect to H. Define Ω 0 := Ω + as the upper part of Ω, and equip Ω 0 with a mesh T 0 .

• For n = 0, . . . , until convergence:

1. Compute the solution (u n λ , p n λ ) ∈ V(Ω n ) × L 2 (Ω n ) of Stokes equation (6.7a)-(6.7c), in Ω n , using the mesh T n .

2. Compute the shape derivative of J(Ω n ) (see Theorem 6.2.3) and infer a descent direction θ n for the optimization problem.

3. Choose an appropriate gradient step τ n and advect the shape Ω n into the new shape Ω n+1 := (Id + τ n θ n )(Ω n ); a mesh T n+1 of Ω n+1 is obtained.

Convergence is reached whenever

|J(Ω n+1 ) -J(Ω n )| τ n < ε stop .
At each iteration, τ n is initialized to a fixed value and divided by a constant factor until J((Id+τ n θ n )(Ω n )) < J(Ω n ). The mesh is finally advected point by point: the i th point x n i of mesh T n is moved to

x n+1 i := x n i + τ n θ n (x n i ) .
(6.17)

The mesh advection step will result in a valid mesh only if θ n is smooth enough and τ n is small enough. If not, self-intersections can appear. We address this difficulty by introducing an extension/regularization procedure.

Conclusion and perspectives

In this chapter, we have applied shape optimization techniques on a model of aquaporin initially studied by soft matter physicists in [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF][START_REF] Belin | Optimal shape of entrances for a frictionless nanochannel[END_REF], in the aim of identifying the geometry of a (local) minimizer of the energy dissipation. The model mixes boundary conditions of different types, such as Dirichlet and partial slip conditions, which are very common in recent models for nanofluidics but complicate the analysis of the problem. We have designed and implemented a numerical optimization method based on the construction of an elastic displacement of the mesh, obtained by exploiting the expression of the shape derivative to define a regular normal displacement of the boundary, and extend it to the whole domain.

In the continuity of this work, we have started to consider with Y. Privat different mechanical criteria to optimize in a similar context, such as the outgoing flow or the penalized dissipated energy, incorporating the penalization term responsible for the implementation of the impermeability condition. We are also working on the generalization of the method to 2daxisymmetric models, which are often used in nanofluidic literature since they permit to simulate cylindrical domains (such as carbon nanotubes) with a relatively low computational cost.

1 2r

 1 (h-β) 4r 2 hλ. Finally, letting β → 0 and λ → 1 we get the density estimate m(B(x, r)) ≥ 2r,

Figure 2 . 1 :Figure 2 . 2 :

 2122 Figure 2.1: Example with 3 points with the set of parameters P = , ε = 6/P , δ t = ε and λ ε = ε 2 . The first and second lines correspond respectively to the operators -ε∆ and -ε 3 ∆ 2 . From left to right: the phase field function u at different times t during the iterations. For the last picture, we plot the stationary solution obtained after dividing the ε parameter by 2.

Figure 2 . 3 :

 23 Figure 2.3: Numerical experiment with 16 points, with the set of parameters P = 2 8 , δ t = ε and λ ε = ε 2 . The first and second lines correspond respectively to the geodesic term R ε , R ε,max . Left to right : stationary solution u ε obtained respectively with different value of ε: ε = 6 √ 2/P , ε = 6/P , ε = 3 √ 2/P and ε = 3/P .

Figure 2 . 4 :

 24 Figure 2.4: Numerical experiments with a large number of points, with the set of parameters P = 2 8 , δ t = ε and λ ε = ε 2 . The first and second lines correspond respectively to the case of 50 and 100 points randomly distributed in the computation box Q. Left to right: stationary solution u ε obtained with different values of ε: ε = 6 √ 2/P , ε = 6/P , ε = 3 √ 2/P and ε = 3/P .

Figure 2 . 5 :

 25 Figure 2.5: Numerical experiments in 3d with the set of parameters P = 2 7 , δ t = ε 2 and λ ε = ε 2 . The first and second lines correspond respectively to the case of 3 and 10 points randomly distributed in the computation box Q. From left to right: initial red points a i , stationary solution u ε obtained with ε = 8/P , ε = 4/P and with ε = 2/P .

Figure 2 . 6 :

 26 Figure 2.6: Numerical experiments in 3d : the cube case with P = 2 6 , δ t = ε 2 and λ ε = ε 2 and ε = 5/P . Each pictures corresponds to the solution along the iterations until convergence to a stationary solution. Last picture picture corresponds to the stationary solution obtained with ε = 2.5/P .

  slip in micro-and nanofluidic models via homogenization This chapter gathers results obtained in collaboration with D. Gérard Varet and A. L. Dalibard in [B12] and with J. Olivier in [B1, B2].

e 1 :

 1 riblets parallel to the flow; -f = e 2 : riblets perpendicular to the flow. The issue is to derive a wall law for the system (3.1)-(3.2)-(3.3), i.e. to replace the mixed boundary condition (3.3) by a condition which does not depend on ε. We show that u ε behaves asymptotically like the solution ū of (3.1)-(3.2), endowed either with a Navier boundary conditionu 3 = 0, ∂ 3 u = M u on Γ b ,(3.6)or with a Dirichlet boundary conditionu = 0 on Γ b . (3.7)In (3.6), M is a 2 × 2 non-negative matrix, whose eigenvalues have the dimension of the inverse of a length. If M = λ I, the number λ -1 is called the "slip length". In the general case, the inverse of the eigenvalues provide the slip lengths in the directions of the eigenvectors. We note ūM the solution of (3.1)-(3.2)-(3.6) and ū0 in the special case M = ( 0 0 0 0 ). We write ū∞ the solution of (3.1)-(3.2)-(3.7).

Figure 3 . 1 :

 31 Figure 3.1: Patch configuration. For every k= (k 1 , k 2 ) ∈ [|0, . . . , ε -1 |] 2 , the intersection of the no-slip zone T ε with the cell [εk 1 , ε(k 1 +1))×[εk 2 , ε(k 2 +1)) is defined by εk + T ε = εk + a ε T .

Figure 3 . 2 :

 32 Figure 3.2: Riblet configuration. For k = (k 1 , k 2 ), the intersection of the no-slip zone T ε with the cell [εk 1 , ε(k 1 + 1)) × [εk 2 , ε(k 2 + 1)) is defined by εk + T ε = εk + (εT) × ( ε 2 + a ε I).

  riblets perpendicular to the main flow), then ūM rib 1

Figure 3 . 3 :

 33 Figure 3.3: Numerical value of the average slip length V ε,∞ e 1 • e 1 plotted against 1/ φ s ε , for circular and squared patches.

) 0 w

 0 (y , y 3 ) dy 3 • ∇ y η(y ) dy = 0 . (3.41) Boundary condition satisfied by ū on Γ b

  ) that results from the energy bound and Dirichlet boundary condition (3.25). Case 0 < α < 1: ū = 0 on Γ b By (3.32) and Korn inequality (3.29), there holds Bε |∇u ε | 2 ≤ C ε -α . Recall the trace inequality

(0, 1 ) 2 ūε, 3 ∂ 3

 1233 (y , εd(y )) η(y ) dyūε,3 (y , y 3 ) dy 3 η(y ) dy .

ε→0 2 ΩεD

 2 (u ε ) : D(φ) + 2ε α Bε D(u ε ) : D(φ) = 2 Ω D(ū) : D(φ) . • By Lebesgue dominated convergence theorem, one easily gets lim ε→0 Bε f • φ = 0 .

56 ) 2 d(y ) 0 ∂ y 3 ww

 56203 Using formula(3.48), the boundary conditions(3.52) and the coupling relation(3.44), we can express the integral ω ∂ y 3 w •∂ y 3 φ by computing (0,1) (y , y 3 ) • ∂ y 3 φ (y , y 3 ) dy 3 dy = (y , y 3 ) • φ (y , y 3 ) dy 3

4. 4 . 1

 41 Proof of existence theorem 4.3.1

(4. 57 )Figure 4 . 1 : 1 - 8 • 10 - 2 - 6 • 10 - 2 - 4 • 10 - 2 - 2 Figure 4 . 2 :

 57411810261024102242 Figure 4.1: Functions Ψ associated with V -shaped, U -shaped and blade riblets (from left to right).

Figure 4 . 3 :

 43 Figure 4.3: Kinetic energy, dissipated energy composed of viscous and frictional dissipations, and average streamwise velocity for different riblet geometries, against roughness parameter λ.

  is a weak solution to system (5.5)-(5.12) if (u ε , w ε ) satisfy boundary conditions (5.8)-(5.10), div u ε = 0 in Ω and the following relation holds for any

. 31 ) 5 . 3 . 3 . 1 0Θ

 315331 Theorem In the critical case δ = 3 2 -1 2 , the pressure p satisfies the Reynolds equation λ ∂ y 1 p(y 1 ) ∂ y 1 θ(y 1 ) dy 1 = 1 0

4 Figure 5 . 1 : 10 Figure 5 . 2 : 10 Figure 5 . 3 :

 45110521053 Figure 5.1: T half plotted against N , for δ = 1 and different values of R c ∈ {0.025, 0.05, 0.1, 0.2}, with νb ∈ {0.1, 0.2, 0.4}. First line: E = 0, second line: E = 10.

Figure 6 . 1 :

 61 Figure 6.1: Scheme of an aquaporin. Water molecules are transported in a single file across the cell membrane, through an hourglass-shaped nanopore. Image courtesy of François Detcheverry (ILM, Lyon 1).

Proposition 6 . 1 . 2 .

 612 denote by H the hyperplane H = {x 2 = 0}, and by Ref H the reflexion through H. We define H + = {x 2 > 0}, Ω + = Ω ∩ H + the upper part of the domain, and Γ sym = Ω ∩ H its lower boundary. Take Ω symmetric with respect to H, and assume that the solution (u, p) to the Stokes system (6.7a)-(6.7b)-(6.7c) belongs toH 2 (Ω, R 2 ) × H 1 (Ω). Then, (u, p) is symmetric in the sense that u = Ref H (u • Ref H ) and p = p • Ref H a.e. in Ω, (6.10) and as a consequence, [σ(u, p)ν] τ = 0 and u • ν = 0 on Γ sym . (6.11)

  Defining the class of admissible shapes byO ad = {Ω open connected with a Lipschitz boundary, Γ in ∪ Γ out ∪ Γ * 0 ∪ Γ 2 ⊂ ∂Ω} ,

Figure 6 . 3 :

 63 Figure 6.3: The external box D contains the domain Ω (in blue), completed with two lateral subregions (in pink).

Figure 6 . 4 :

 64 Figure 6.4: Criterion J(Ω) with respect to angle θ in test-cases 1 and 2 (see Table6.1).

Figure 6 . 5 :

 65 Figure 6.5: From top to bottom : mesh at iterations 0, 45 and 85 for Testcase 1.

Figure 6 . 6 :

 66 Figure 6.6: From top to bottom : mesh at iterations 0, 130 and 245 for Test-case 2.

  Minimization over pairs (u, γ) . . . . . . . . . . . . . . . . . . 2.2.1 Dissociation and regularization of functional F ε . . . . 2.2.2 Elements of proof . . . . . . . . . . . . . . . . . . . . . 2.3 Numerical approximation . . . . . . . . . . . . . . . . . . . . 2.3.1 Modified phase field models . . . . . . . . . . . . . . . 2.3.2 Numerical experiments . . . . . . . . . . . . . . . . . . 2.4 Conclusion and perspectives . . . . . . . . . . . . . . . . . . . 3 Slip in micro-and nanofluidic models 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Rough hydrophobic surface model . . . . . . . . . . . . . . .

	vi

1 Introduction 1.1 What is contained in this thesis . . . . . . . . . . . . . . . . . 1.2 What is not contained in this thesis . . . . . . . . . . . . . . 1.3 General notation and conventions . . . . . . . . . . . . . . . . 2 Approximation of Steiner tree problem 2.1 Phase-field approximation enforcing connectedness . . . . . . 2.1.1 Main result . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Elements of proof . . . . . . . . . . . . . . . . . . . . . 2.2

  + denotes the average of f t on [g t , hr] and m -denotes the average of f t on [-hr, g t ] (here we identified g t with its coordinate on the second axis). Next,

	r)	∂ ∂x 2	(P (ϕ ε )) dx ≥	λr -λr	Var(f t , [-hr, hr]) dt
	where f t := P (ϕ ε )| Lt with L t := t + Re 2 . On the other hand applying
	Lemma 2.1.4, for every t ∈ I λ,r,ε , we can write
	Var(f t , I t ) ≥ m + t + m -t -2f t (g t ),	(2.23)
	where m Var(f t , [-hr, hr]) ≥	1 r(h + β) It	f t (s) ds -2f t (g t ).
	Integrating over t ∈ I λ,r,ε = [-λr, λr] \ [t L , t R ] and applying Fubini's
	Theorem we get

Assume that spt µ is not reduced to {a 0 } and that H 1 (spt µ) < ∞. Assume also that (2.30) holds. Let ε k ↓ 0 and {u k

	) .	(2.30)
	Provided that H 1 (spt µ) < ∞, our second main result shows that sublevel
	sets of minimizers converge to a solution of the generalized Steiner problem
	min H 1 (K) : K ⊂ R 2 compact and connected, K ⊃ {a 0 }∪spt µ . (2.31)
	Note that for µ = N i=0 δ a i and some distinct points a i ∈ Ω, problem (2.31)
	coincides with the classical Steiner problem (2.1).
	Theorem 2.2.2.

  Proposition 3.3.2. Let f be in L 2 (Ω, R 3 ). Then for any 0 < α ≤ 2, the weak solution u ε to (3.30) satisfies the energy bound (3.32).

	Proof. Assuming 0 < α ≤ 2, combining Poincaré inequalities (3.26)-(3.27)
	with Korn inequalities (3.28)-(3.29), and testing against u ε in (3.31), we
	obtain
			.31)
	Definition 3.3.1 (Energy bound). We say that a family of solutions u ε ∈
	V D to (3.31) satisfies the energy bound if there exists a constant C > 0 such
	that
	|D(u ε )| 2 + ε α	|D(u ε )| 2 ≤ C .	(3.32)
	Ωε	Bε

Table 6 .

 6 1: Model and numerical parameters simplify notation, we present all the material on the full domain Ω, but the calculations remain valid on the symmetric problem, just by replacing Ω by Ω + , Γ 1 by Γ 1 ∩ H + , etc.

	3. In order to

2.2. MINIMIZATION OVER PAIRS (u, γ)

= O ε ,

Remerciements

Chapter 5

Generalized Reynolds equation in micropolar lubrication over a ribbed surface

This chapter presents results obtained in collaboration with I. Pažanin and F. J. Suárez-Grau in [B3].

Theorem 5.3.2. Assume that the asymptotic regimes (5.17) and (5.18) and conditions (5.21) and (5.22) hold. Assume that δ, satisfy the relation δ = 3 2 -1 2 (critical case). Let (u ε , w ε , p ε ) be a sequence of weak solutions of (5.5)- (5.12). Then, there exist u , w ∈ H 1 (0, h; L 2 (ω, R 2 )) and p ∈ H 1 (ω) ∩ L 2 (ω)/R such that the rescaled functions u ε , w ε , p ε satisfy u ε|Ω ( u , 0) and ε w ε|Ω ( w , 0) in H 1 (0, h; L 2 (ω, R 3 )),

(5.26)

The triplet ( u , w , p) is the unique solution of the following problem

Coefficients E, F ∈ R appearing in boundary conditions (5.29) are defined by

where ( φ i , q i ), i = 1, 2, are the respective solutions to

Proof of Proposition (5.4.1). Estimates (5.36)-(5.37) are obtained classically by testing against (ϕ, ψ) = (u ε -J ε , w ε ) in (5.33)- (5.34). The pressure estimate (5.38) is more specific to the thin domain configuration: it relies on an estimate of ∇p ε in H -1 and on the following inequality, proved in [START_REF] Casado-Díaz | Asymptotic behavior of the Navier-Stokes system in a thin domain with Navier condition on a slightly rough boundary[END_REF]Corollary 4.2]:

Convergence of the velocity and microrotation

As a consequence of the a priori estimates stated in Proposition 5.4.1, and the fact that Ω ⊂ Ω ε and | Ω ε \ Ω| → 0, we have the following convergence results for rescaled solutions u ε , w ε restricted to the limit domain Ω.

Lemma 5.4.2. Assume that the asymptotic regimes (5.17) and (5.18) and conditions (5.21) and (5.22) hold. There exist u , w in H 1 (0, h; L 2 (ω, R 2 )) with u (x , h) = w (x , h) = 0 for x ∈ ω, and

such that, up to extraction,

Since Lemma 5.4.2 is specific to the case of a thin domain developed in this chapter, we provide some details on its proof.

Proof. Rescaling estimates (5.36)-(5.37), one obtains

(5.42)

Hence, there exist u and w such that, up to extraction, (5.40) holds. The conditions u ε (x , h) = -εSe 3 and ε w ε (x , h) = 0 pass to the limit, yielding u(x , h) = w(x , h) = 0 for x ∈ ω.

The fact that u 3 ≡ 0 is a consequence of the incompressibility condition on u ε , expressed on the rescaled function u ε by

Indeed, multipying the previous equality by φ ∈ C ∞ c (Ω), and integrating by parts over Ω, we deduce

Using (5.40), we conclude that u 3 does not depend on y 3 , and since it vanishes on y 3 = h, it is identically null.

Next, we prove the divergence equation (5.39). Using condition (5.6), integration by parts, boundary conditions (5.8), (5.10) and the change of variables (5.23), we have for any

By the bound (5.41) and Hölder inequality, we can replace the integral over Ω ε by an integral over Ω, obtaining

Passing to the limit by the weak convergence (5.40), we deduce

which proves (5.39).

It remains to prove that w 3 = 0. Since w ε is not incompressible, we use the angular momentum equation. Let ψ ∈ C ∞ c (Ω) and define ψ ε = εψ(x , x 3 /ε)e 3 . Testing against (0, ψ ε ) in (5.33) and using (5.40), we get

(5.44)

By the no-penetration condition w ε • ν = 0 on Γ 0 ε , we also deduce that w 3 (y , 0) = 0 for x ∈ ω. Since w 3 vanishes on y 3 = h, we can take ψ = w 3 in (5.44) and easily deduce the thesis.

Convergence of the pressure

The convergence result for the rescaled pressure p ε is built on an accurate estimate for pressure p ε , obtained using the following decomposition result [START_REF] Casado-Díaz | Asymptotic behavior of the Navier-Stokes system in a thin domain with Navier condition on a slightly rough boundary[END_REF]Corollary 4.2]: Proposition 5.4.3. The pressure p ε ∈ L 2 (Ω ε )/R can be decomposed as

Derivation of the model

The fluid domain Ω ε given by (5.3) is now dependent on time since ε = ε(t), and parameters λ and δ used to describe the rough geometry also depend on time. However, since the shape of the rough plate is fixed, the relations

necessarily hold for some constants Λ, M which correspond respectively to the amplitude and period of the riblets.

Considering the critical regime δ(t) = 3 2 (t) -1 2 , the parameters (t), δ(t), λ(t) can thus be expressed as functions of ε(t), M, L, and in particular the parameter λ(t) simplifies to

Assuming that a constant load W is applied on the upper plate of the bearing, the rescaled pressure p ε satisfies ω p ε(t) (y , h) dy = W .

(5.52)

Since the system is independent on the x 2 -direction, we take ω = (0, 1) and apply Theorems 5.3.2 and 5.3.3 to approximate p ε(t) by p λ(t),S(t) /ε(t) 2 , where p λ(t),S(t) is the solution of Reynolds equation (5.32), with λ = λ(t) and S = S(t), S(t) being fixed consistently with relation (5.52). As is usual in the lubrication field, we impose Dirichlet boundary conditions for the pressure on x 1 ∈ {0, 1}, instead of the Neumann boundary conditions implicitly contained in the weak formulation (5.32). We obtain the relation

(5.53) Since Reynolds equation (5.32) is linear with respect to S, there holds p λ(t),S(t) = S(t) p λ(t),1 (where p λ(t),1 satisfies (5.32) with λ = λ(t) and S = 1) so that S(t) is given by

(5.54)

The normal velocity of the upper plate being equal to -ε(t)S(t), the motion of the bearing is driven by the differential equation

where λ(t) depends on ε(t) through relation (5.51). Finally, expressing p λ(t),1 by the formula p λ(t),1 (y 1 ) = y 1 (1-y 1 )

2Θ λ(t)

and setting h = 12W for simplicity, the ODE (5.55) can be rewritten

(5.56)

Equations of motion Denoting κ the dimensionless constant κ = Λ 2 M 3 , the model can be summarized by the system of equations Hence, in addition to the initial datum ε 0 , the solution t → ε(t) to the system (5.57) depends on the following set of parameters:

• parameters N, R c characterizing the mechanical properties of the micropolar fluid,

• νb , δ characterizing the fluid/surface interaction,

• E, κ related to the geometry of the rough pattern.

In the aim of estimating the influence of the parameters on the performance of the squeeze-film bearing, we solve the associated system (5.57) by a second-order Runge-Kutta method, and compute the "half-life time" T half , i.e. the first instant t such that ε(t) < ε 0 /2, which means that the width of the bearing has been divided by two. In our analysis, the configurations giving rise to the highest values of T half will be considered the most efficient.

For a given Q ∈ R \ {0}, we consider the following problem: find (u λ , p λ ) ∈ V(Ω) × L 2 (Ω) and λ ∈ R such that

Remark 6.1.1. Parameter λ appearing in (6.7a) can be regarded as the Lagrange multiplier associated to the constraint Γ in u λ • ν dH 1 = -Q, and expressed as

(u 1 , p 1 ) being the solution of (6.7a)-(6.7b) with λ = 1.

The next proposition interprets the solution u λ of the Stokes system as the unique minimizer of a certain energy. Proposition 6.1.1. For every Q ∈ R \ {0}, there exists a unique triple (u λ , p λ , λ) ∈ V(Ω) × L 2 (Ω) × R satisfying (6.7a)-(6.7b)-(6.7c). Moreover, the function u λ is the unique minimizer of the energy functional E Ω defined by

over the space

We end this section by pointing out some symmetry properties satisfied by the solution in the case where domain Ω is symmetric, that we use in the numerical simulations.

Mixed formulation of the Stokes problem with a symmetry condition.

In this paragraph, we summarize the symmetry properties enjoyed by the model if Ω is assumed symmetric with respect to the axis {x 2 = 0}. We

Construction of the descent direction θ n

This paragraph constitues the core of step 2. The method presented here relies on a H 1 -regularization step [START_REF] Dogan | Discrete gradient flows for shape optimization and applications[END_REF][START_REF] Burger | A framework for the construction of level set methods for shape optimization and reconstruction[END_REF][START_REF] Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF], followed by a linear elasticitybased extension. To our best knowledge, such extension is not standard, but offers the advantage to take into account Γ 1 j 2 • ∇ τ (θ • ν) dH 1 without assuming additional regularity on factor j 2 .

Substep 1: regularization procedure Let φ ∈ H 1 0 (Γ 1 ) be the solution of the following variational formulation:

for all ψ ∈ H 1 0 (Γ 1 ) . (6.18) Setting ψ = φ yields the inequality

Substep 2: extension to the whole domain We look for a smooth vector field θ defined on Ω and satisfying the boundary condition

Combining (6.19) with definition (6.16) of the shape derivative, the previous condition ensures that θ is a descent direction.

To construct this extension, we adopt a widespread strategy in meshing, based on the intuitive idea that elastic deplacements tend to reduce compression (see for instance [START_REF] Baker | Mesh movement and metamorphosis[END_REF][START_REF] Dapogny | Geometrical shape optimization in fluid mechanics using FreeFem++[END_REF]). More specifically we choose θ as the unique solution in H 1 (Ω, R 2 ) of the linear elasticity problem

where σ e (θ) stands for the elasticity tensor given by σ e (θ) = µ e D(θ) + divθ.

Numerical results

The evolution of the mesh during Step 2 of the algorithm is displayed in Figs. [START_REF] Allaire | Homogenization and two-scale convergence[END_REF].5 (test-case 1) and 6.6 (test-case 2). This step appears highly beneficial, leading to a 35% decrease of J in case 1, and 40% in case 2 respectively. Left: evolution of the cost functional J, and its subparts J d , J f over iterations, relatively to the initial value J 0 of J. Right: evolution of the stopping criterion.

To take the analysis a step further, let us investigate which term in J contributes most to J between • the viscous dissipation J v (Ω) := 2µ Ω |D(u Ω,λ )| 2 dx,

• the dissipation by friction J f (Ω) := β Γ 1 |u Ω,λ | 2 dH 1 .

In Figs. 6.7-6.8, we observe that the reduction of the total dissipated energy achieved by step 2, results from an important decrease of the viscous dissipation J d , which is the main contributor to the cost functional J. The frictional dissipation term J f accounts for less than 10 percent of the total dissipation, and is slightly increased during the process.