Alain Bastide

Pierre-Henri Cocquet

Delphine Ramalingom

PENALIZATION MODEL FOR NAVIER-STOKES-DARCY EQUATION WITH APPLICATION TO POROSITY-ORIENTED TOPOLOGY OPTIMIZATION

Keywords: Navier-Stokes-Darcy model, Porous media, PDE-constrained optimization, Topology optimization AMS subject classifications. 35Q30, 76S05, 49J20, 76D55

HAL is

1. Introduction. Many practical problem in applied science and engineering aim at finding the location of a solid to improve a desired physical behavior. Such kind of problem, usually referred as topology optimization, can be written as a constrained optimization problem where the parameter we look for is the location of the solid, the constraint are given by the governing equation of the fluid and the function one wants to minimize describe a physical effect.

Several techniques have been proposed to use the location of the porous media as optimization parameter. For instance, the topological asymptotic expansion [START_REF] Amstutz | The topological asymptotic for the Navier-Stokes equations[END_REF][START_REF] Amstutz | Topological derivatives for a class of quasilinear elliptic equations[END_REF] consider the solid as holes in the computational domain. The topological gradient is then the first order term of the asymptotic expansion of the cost functional as the size of the hole vanishes. One of the advantage of this method is that the solid is clearly taken into account in the optimization process. However, this method requires the computation of an asymptotic expansion which can be a hard task that depend on the shape of the hole, the dimension of space and the problem under study. For these reasons, the penalization method that has been introduced in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] has been favored over the past few years.

The method introduced by T. Borrvall and J. Petersson [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] is based on a penalization model [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] that define the location of the porous media with a friction term α(ρ)u where u is the velocity of the fluid and α(ρ) = ψ(ρ -1 -1) is (physically) an inverse of a permeability by a kinematic viscosity with ψ a dimensioned coefficient. The fluid part of the optimal design is then defined when ρ = 1 while ρ = 0 gives the solid part since, at least formally, the velocity of the fluid goes to zero in these zones. As pointed out by A. Evgrarov [START_REF] Evgrafov | The limits of porous materials in the topology optimization of Stokes flows[END_REF], the parameter ρ belong to a non-convex nor weakly * closed set thus both the mathematical study and the numerical approximation of the topology optimization problem is difficult. It has then been proposed (see e.g. [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF][START_REF] Wiker | Topology optimization of regions of Darcy and Stokes flow[END_REF]) to work with an interpolation of α hence to use ρ ∈ [0, 1] as design parameter. This interpolation technique has then been used extensively in several work [START_REF] Wiker | Topology optimization of regions of Darcy and Stokes flow[END_REF][START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF][START_REF] Evgrafov | On Chebyshevs method for topology optimization of Stokes flows[END_REF][START_REF] Evgrafov | The limits of porous materials in the topology optimization of Stokes flows[END_REF][START_REF] Andreasen | Topology optimization of microfluidic mixers[END_REF][START_REF] Gersborg | Topology optimization of flow problems[END_REF][START_REF] Othmer | A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows[END_REF] (see also the review paper [START_REF] Dbouk | A review about the engineering design of optimal heat transfer systems using topology optimization[END_REF]) and extended to other physical modeling like heat transfert in fluid [START_REF] Qian | Topology optimization of a coupled thermal-fluid system under a tangential thermal gradient constraint[END_REF][START_REF] Kontoleontos | Adjoint-based constrained topology optimization for viscous flows, including heat transfer[END_REF][START_REF] Ramalingom | Optimisation topologique des échanges thermiques dans un canal vertical asymétriquement chauffé[END_REF] or natural convection [START_REF] Alexandersen | Topology optimisation for natural convection problems[END_REF].

The major drawback of the technique introduced by T. Borrvall & J. Petersson is that the penalization model is relevant only for perfect solid that is a material that does not allow the fluid to enter. In this paper, we want to propose an alternative method to solve topology optimization problem that allows the optimal design to be characterized by a porous medium and not only a perfect solid. Following the same idea as the penalization method for which we have a solid as α goes to infinity, we would like to use the porosity α as optimization parameter and get the Navier-Stokes equations when α = 0 while, for positive α, we wish to have the Darcy law for porous medium with porosity α. Our new penalized model is going to be given by a singular perturbation of the Navier-Stokes problem whose weak limit is the Navier-Stokes-Darcy model. The techniques used to build our penalized model are based on some vanishing-viscosity method similar to those developed in [START_REF] Gastaldi | On the coupling of two dimensional hyperbolic and elliptic equations: analytical and numerical approach[END_REF][START_REF] Gander | Viscous problems with inviscid approximations in subregions: a new approach based on operator factorization[END_REF] for advection-diffusion problem.

This paper is organized as follow: We first derive the singular perturbation model used latter as constraint in our topology optimization problem. We show next that its weak-limit is a fluid-porous media interface problem given by a Navier-Stokes-Darcy problem. We then use this penalization model as constraint for a topology optimization problem and give necessary condition for the existence of optimal design as well as first order necessary optimality condition. This paper ends with some numerical simulation for Stokes flow to show the potential application of the approach developed in the present paper.

2. The penalized Navier-Stokes-Darcy model. Let Ω be a bounded open set of R d , d = 3 2 , with Lipschitz boundary Γ = Γ 1 ∪ Γ 2 (Γ 1 ∩ Γ 2 = ∅) whose outward unitary normal is denoted by n. We assume that Ω can be written as Ω = Ω p ∪ Ω f , where the subscript f denotes the fluid part and p the porous part of the computational domain. We assume that one has the incompressible Navier-Stokes equation in Ω f and the Darcy law for porous medium with porosity α in Ω p . This reads

       -ν∆u f + 1 ρ f ∇p f + (u f • ∇)u f = f in Ω f , div u f = 0 in Ω f αu p + 1 ρ f ∇p p = f in Ω p , div u p = 0 in Ω p , (2.1)
with suitable interface condition on (∂Ω f ∩ ∂Ω p ) \ Γ and some boundary condition on Γ that will be specified latter. In (2.1), ρ f is the constant fluid volumetric mass and is equal to 1 in the sequel of the paper. Our first goal is to obtain Equation (2.1) as the limit of a model involving a single equation giving the Navier-Stokes equation for α = 0 and, at the limit, the Darcy law when α > 0. The latter can be achieved using function of the form exp(-τ α) which converges, as τ → +∞, toward the indicator function of the set {x ∈ Ω | α(x) = 0} := Ω f . Therefore, multiplying the second order term and the non-linear term of Equation (2.1) by exp(-τ α) is going, at least formally, to do the job. However, this method yields a problem whose mathematical analysis can be difficult and this is the reason why we propose below another approach based on the following formula

(u • ∇)u = 1 2 ∇|u| 2 -u × (∇ × u), (2.2)
Using (2.2), one can write the incompressible Navier-Stokes equations as follow

-ν∆u + ∇p + 1 2 ∇|u| 2 -u × (∇ × u) = f in Ω f .
We then consider the following penalized incompressible Navier-Stokes system with Darcy law for porous medium

       -div νe -τ α ∇u -pI 3 -1 2 e -2τ α |u| 2 I 3 -e -2τ α u × (∇ × u) + αu = f in Ω, div(u) = 0 in Ω, u = 0 on Γ 1 , νe -τ α ∂ n u -(p + 1 2 e -2τ α |u| 2 )n = ϕ on Γ 2 , (2.3)
2 The case Ω ⊂ R 2 is discussed along the paper and does not yield additional difficulties.

where div denotes alternatively the divergence of a tensor or a vector field, α ∈ L ∞ (Ω) is a kinematic viscosity divided by a permeability, f ∈ L 2 (Ω) 3 is a density of forces and ϕ ∈ L 2 (Γ 2 , R 3 ) is a surface source term. We now derive formally the limiting problem as τ → +∞. Note that we recover the Navier-Stokes equations on Ω f since α vanishes. Also taking the formal limit as τ → +∞ of the solution to Problem (2.3), we get that it satisfies

div u = 0, αu + ∇p = f in Ω p ,
which is the Darcy law for porous media. Problem (2.3) is therefore a penalization model for an interface fluid-porous media problem. Remark 2.1 (Two-dimensional case). In a two-dimensional setting, one only has to use that

(u • ∇)u - 1 2 ∇|u| 2 = u 2 ∂ 2 u 1 -u 2 ∂ 1 u 2 u 1 ∂ 1 u 2 -u 1 ∂ 2 u 1 .
Therefore, we simply replace ∇× by the above formula in Problem (2.3) to get the penalized Navier-Stokes-Darcy equation in R 2 .

We start by studying the existence of solution to Problem (2.3). Let X be the subspace of H 1 (Ω) 3 defined by

X := ψ ∈ H 1 (Ω) 3 | ψ = 0 on Γ 1 .
It is a Hilbert space for the norm . X := . H 1 (Ω) d . We recall the following Poincaré inequality that holds if

|Γ 1 | > 0 ∃C P > 0 such that ∀ϕ = (ϕ j ) ∈ X, ϕ j L 2 (Ω) ≤ C P ∇ϕ j L 2 (Ω) 3 , (2.4)
where the constant C P only depends on Ω. In the sequel, we use the following notation ∀u ∈ X, ∇u

2 L 2 (Ω) d := d j=1 ∇u j 2 L 2 (Ω) d ,
which is actually the norm of the Jacobian matrix of the vector field x ∈ Ω → u(x).

A variational formulation for Problem (2.3) is given below

F ind (u, p) ∈ X × L 2 (Ω) such that ∀v ∈ X, a(u, v) + N (u, u, v) + b(v, p) = Ω f (x) • v(x)dx + Γ2 ϕ • vdσ, ∀q ∈ L 2 (Ω), b(u, q) = 0, (2.5)
where, using subscript j to denote the components of a vector, the multilinear forms are defined as follow

a(u, v) = ν d j=1 Ω e -τ α ∇u j • ∇v j dx + Ω αu • vdx, N (u, v, w) = - Ω e -2τ α 2 u • vdiv(w) + e -2τ α u × (∇ × v) • wdx, (2.6) b(v, p) = - Ω divvp dx.
The equivalence between the weak-formulation (2.5) and Problem (2.3) is classical and can be found in [START_REF] Bernardi | Finite element discretization of the Stokes and Navier-Stokes equations with boundary conditions on the pressure[END_REF]. We show below the continuity of the trilinear form associated to the non-linear term N .

Lemma 2.2. The trilinear form N is continuous on

X 3 |N (u, v, w)| ≤ C N ∇u L 2 (Ω) 3 ∇v L 2 (Ω) 3 ∇w L 2 (Ω) 3 ,
where C N > 0 depends only on Ω.

Proof. We recall the following compact Sobolev embedding theorem (see e.g [START_REF] Galdi | An introduction to the mathematical theory of the Navier-Stokes equations: Steady-state problems[END_REF][START_REF] Temam | Navier-Stokes equations: theory and numerical analysis[END_REF])

H 1 (Ω) → → L 4 (Ω). (2.7)
Let (u, v, w) ∈ X 3 . Using (2.7), that e -τ α ≤ 1 and Holder inequality then give

|N (u, v, w)| ≤ 2 u L 4 (Ω) 3 ∇v L 2 (Ω) 3 w L 4 (Ω) 3 ≤ C N ∇u L 2 (Ω) 3 ∇v L 2 (Ω) 3 ∇w L 2 (Ω) 3
which proves that N is continuous. We first study the existence and uniqueness of (u τ , p τ ) satisfying Problem (2.3) for finite 0 < τ . Theorem 2.3. Suppose that the following set of assumptions hold 1. α ∈ L ∞ (Ω) is almost everywhere positive.

2. The domain Ω can be written as

Ω = Ω p ∪ Ω f where Ω f = {x ∈ Ω | α(x) = 0} satisfies Ω f ∩ Γ 1 = ∅.
3. There exists Ω p ⊂ Ω and α min > 0 such that a.e x ∈ Ω p , α min ≤ α(x).

Then for any 0 < τ ≤ M , there exists at least one (u τ , p τ ) satisfying Problem (2.5) and the following a priori estimates

u τ H 1 (Ω) 3 ≤ C 1 e M α L ∞ (Ω) ν f L 2 (Ω) 3 + ϕ L 2 (Γ2) 3 , p τ L 2 (Ω) ≤ β -1 C 2 1 + α L ∞ (Ω) C 1 e M α L ∞ (Ω) ν f L 2 (Ω) d + ϕ L 2 (Γ2) 3 +β -1 C N C 2 1 e M α L ∞ (Ω) ν 2 f L 2 (Ω) 3 + ϕ L 2 (Γ2) 3 2 ,
where C 1 , C 2 are positive constant that does not depend on ν nor on M . Proof. We first prove the upper bound for the solution to Problem (2.3). Note that the non-linear term satisfy

N (u, u, u) = - Ω e -2τ α |u| 2 div(u) + u × (∇ × u) • u dx = Ω e -2τ α (∇ × u) • (u × u)dx = 0.
Taking v = u in the weak formulation (2.5), using Cauchy-Schwarz, Poincaré (2.4) and trace inequality and that e -τ α ≥ e -M α L ∞ (Ω) then show that

u H 1 (Ω) 3 ≤ C 1 e M α L ∞ (Ω) ν f L 2 (Ω) 3 + ϕ L 2 (Γ2) 3 , (2.8)
where C 1 is a constant that does not depend on ν nor M . To get the estimate on the pressure, one needs the next inf-sup condition whose proof can be found in [START_REF] Bernardi | Spectral approximations of the Stokes equations with boundary conditions on the pressure[END_REF] Theorem 2.1, p. 338.

∃β > 0 such that ∀q ∈ L 2 (Ω), sup v∈X b(v, q) v X ≥ β q L 2 (Ω) . (2.9) Since p ∈ L 2 (Ω) satisfies Problem (2.5), we obtain b(v, p) = Ω f (x) • v(x)dx - Γ2 ϕv • ndσ -a(u, v) -N (u, u, v), ∀v ∈ X.
Using the inf-sup condition (2.9) together with Lemma 2.2 then give

p τ L 2 (Ω) ≤ β -1 sup v∈X Ω f (x) • v(x)dx -Γ2 ϕv • ndσ -a(u τ , v) -N (u τ , u τ , v) v X , ≤ β -1 f L 2 (Ω) 3 + ϕ L 2 (Γ2) 3 + (1 + α L ∞ (Ω) ) u τ H 1 (Ω) 3 +β -1 C N u τ 2 H 1 (Ω) 3
which yields the desired result thanks to (2.8).

We now prove the existence of at least one solution to Problem (2.3). First we introduce the kernel of the linear form q ∈ L 2 (Ω) → b(., q) ∈ X namely

V := {v ∈ X | div(v) = 0 in Ω} .
It is worth noting that u ∈ V. We introduce the non-linear mapping Φ : V → V defined by

Φ(u), v := a(u, v) + N (u, u, v) - Ω f (x) • v(x)dx - Γ2 ϕ • vdσ.
We then need to prove that there exists u ∈ V such that Φ(u) = 0 to get the existence of u satisfying Problem (2.5). Since N (u, u, u) = 0, we infer

Φ(u), u ≥ C 1 u 2 H 1 (Ω) 3 -C 2 ( f L 2 (Ω) 3 + ϕ L 2 (Γ2) 3 ) u H 1 (Ω) 3 ,
where C j > 0 does not depend on ν. One thus gets that Φ(u), u is non-negative on the sphere with radius r 0 := C 2 ( f L 2 (Ω) 3 + ϕ L 2 (Γ2) 3 )/C 1 . Since V is separable (as a closed subspace of the Banach space H 1 (Ω) 3 ), there exists an increasing sequence of finite dimensional subspaces V n such that ∪ n V n is dense in V. From Brower's fixed point Theorem, there exists u n ∈ V n such that u n H 1 (Ω) 3 ≤ r 0 and

∀v n ∈ V n , Φ(u n ), v n = 0.
The sequence V n being increasing, one has 3 is compact (see (2.7)), there exists a subsequence (still denoted (u n ) n ) that converges weakly in X and strongly in L 4 (Ω) 3 to some u ∈ N. The definition of the weak convergence then yield

∀m ≥ n, ∀v m ∈ V m , Φ(u n ), v m = 0. Since (u n ) n ⊂ V is bounded in X and the embedding X ⊂ L 4 (Ω)
lim n→∞ a(u n , v m ) = a(u, v m ). Since u n u in X and u n → u strongly in L 4 (Ω) 3 , one gets that u n • v m → u • v m strongly in L 2 (Ω) 3 and thus lim n→∞ N (u n , u n , v m ) = N (u, u, v m ).
From the above limits, we infer

∀m ∈ N, ∀v m ∈ V n , Φ(u), v m = 0.
The density of ∪ n V n into V then gives Φ(u), v = 0 for any v ∈ V and thus the existence of a velocity field for Problem (2.3) is proven. The existence of p ∈ L 2 (Ω) then follow from the inf-sup condition (2.9) and Lemma 4.1 p.38 from [START_REF] Girault | Finite element methods for Navier-Stokes equations: theory and algorithms[END_REF].

Remark 2.4. The existence proof also works for f ∈ X * and ϕ ∈ H -1/2 (Γ 2 , R 3 ) and one only has to change the norm in the upper bound of the solution.

We do not detail here the proof of the unicity since this can be done using standard technique from [START_REF] Temam | Navier-Stokes equations: theory and numerical analysis[END_REF][START_REF] Galdi | An introduction to the mathematical theory of the Navier-Stokes equations: Steady-state problems[END_REF][START_REF] Bernardi | Finite element discretization of the Stokes and Navier-Stokes equations with boundary conditions on the pressure[END_REF] and we are mainly interested in existence when doing PDE-constrained optimization problems. Nevertheless, this demand ν to be large enough.

3. Derivation of the limiting problem. We are going to characterize the limit of (u τ , p τ ) satisfying Problem (2.3) as τ → +∞. We assume from now that the data satisfies f ∈ L 2 (Ω) 3 and ϕ = 0. We require f ∈ L 2 (Ω) 3 because this regularity is useful to give sense to the interface condition as shown in the proof of Theorem 3.2. We also consider vanishing ϕ because this greatly simplify the derivation of some bound proved in the sequel. Finally, we assume that the fluid part Ω f ⊂ Ω and the porous part Ω p ⊂ Ω are bounded open set with Lipschitz boundary since we are going to write (2.3) as an interface problem.

Using that Ω = Ω f ∪ Ω p and omitting the τ -dependance to lighten the expression, Problem (2.3) can be written as the following transmission problem

                               -div νe -τ α ∇u p -p τ,p I 3 -1 2 e -2τ α |u p | 2 I 3 -e -2τ α u p × (∇ × u p ) + αu p = f in Ω p , div(u p ) = 0 in Ω p , -div ν∇u f -p f I 3 -1 2 |u f | 2 I 3 -u f × (∇ × u f ) = f in Ω f , div(u f ) = 0 in Ω f , u f = 0 on Γ 1 ∩ ∂Ω f , u p = 0 on Γ 1 ∩ ∂Ω p , ν∂ n u f -(p f + 1 2 |u f | 2 )n = 0 on Γ 2 ∩ ∂Ω f , e -τ α ν∂ n u p -(p p + 1 2 e -2τ α |u p | 2 )n = 0 on Γ 2 ∩ ∂Ω p , u f -u p = 0 on Γ p,f , T τ [(u p , p p ), (u f , p f )] = 0 on Γ p,f , (3.1)
where Γ p,f = (∂Ω p ∩ ∂Ω f ) \ Γ is the fluid-porous media interface. The application T τ defines the flux transmission conditions. Chosing n pf as unitary normal to Γ p,f , for instance inward to Ω p and outward to Ω f , the latter is given by

T τ [(u p , p p ), (u f , p p )] = (νe -τ α ∇u p -p p I 3 -1 2 e -2τ α |u p | 2 I 3 )n pf -ν∇u f -p f I 3 -1 2 |u f | 2 I 3 n pf . (3.2)
Note that Theorem 2.3 ensures the existence of (u s,τ , p τ,p ,

u f,τ , p τ,f ) ∈ H 1 (Ω p ) 3 ×L 2 (Ω p )×H 1 (Ω f ) 3 ×L 2 (Ω f ) 3
that satisfy the transmission Problem (3.1-3.2). Before giving the variational formulation of Problem (3.1-3.2), we introduce the following functional spaces

X f := ψ ∈ H 1 (Ω f ) 3 | ψ = 0 on Γ 1 ∩ ∂Ω f . X p := ψ ∈ H 1 (Ω p ) 3 | ψ = 0 on Γ 1 ∩ ∂Ω p .
The variational formulation then reads

F ind (u τ,f , p τ,f , u p , p τ,p ) ∈ X f × L 2 (Ω f ) × X p × L 2 (Ω p ) such that u f -u p = 0 on Γ p,f and a f (u τ,f , v f ) + N f (u τ,f , u τ,f , v f ) + b f (v f , p τ,f ) + a p (u τ,p , v p ) + N p (u τ,p , u τ,p , v p ) +b(v p , p τ,p ) - Γ p,f νe -τ α ∇u τ,p -p τ,p I 3 - 1 2 e -2τ α |u τ,p | 2 I 3 n pf • (v p -v f )dσ = Ω f f (x) • v f (x)dx + Ωp f (x) • v p (x)dx b f (u τ,f , q f ) + b p (u τ,p , q p ) = 0, ∀(v f , v p ) ∈ X f × X p ∀(q f , q p ) ∈ L 2 (Ω f ) × L 2 (Ω p ). (3.3)
Above the a f , a p , N f , N p , b f , b p are defined as in (2.6) with Ω replaced by Ω f and Ω p according to their subscript.

3.1. First estimates. We now prove some estimates that are sharper than those from Theorem 2.3 and are going to be useful to pass to the limit in Problem (3.1-3.2).

Theorem 3.1. Assume that the assumptions of Theorem 2.3 hold and that α| Ωp ∈ Lip(Ω p ). Then the solution to Problem (3.1) satisfy the following bounds

u f,τ L 2 (Ω f ) 3 + u p,τ L 2 (Ωp) 3 ≤ C γ f L 2 (Ω) 3 , ∇u f,τ L 2 (Ω f ) 3 ≤ C γ f L 2 (Ω) 3 , e -τ α ∇u p,τ L 2 (Ωp) 3 ≤ e -τ αmin/2 C √ νγ f L 2 (Ω) 3 , p p,τ L 2 (Ωp) + p f,τ L 2 (Ω f ) 3 ≤ B,
where γ = min(C -1 P ν, α min ), C > 0 is a generic constant that does not depend on τ nor on ν and B > 0 is a constant that does not depend on τ .

Proof. We omit the τ dependence to lighten the overall expressions. We now choose (v f , v p ) ∈ X f × X p defined as v f = u f , v p = u p in the variational formulation (3.3). Using then that N t (w t , w t , w t ) = 0 for any w t ∈ X f with subscript t ∈ {f, p} referring to the fluid or the porous part of Ω, one gets

3 j=1 Ωp νe -τ α |∇u s,j | 2 dx + Ω f ν|∇u f,j | 2 dx + Ωp α|u p | 2 dx = Ωp f • u p dx + Ω f f • u f dx
Using Cauchy-Schwarz, trace and Poincaré inequalities (2.4), we obtain

α min u p 2 L 2 (Ωp) 3 + C -1 P ν u f 2 L 2 (Ω f ) 3 ≤ f L 2 (Ω) 3 u p L 2 (Ωp) 3 + u p L 2 (Ωp) 3 .
Since 2(a 2 + b 2 ) ≥ (a + b) 2 for any a, b ∈ R + , we get the existence of a constant C > 0 that does not depend on ν nor τ such that

u f L 2 (Ω f ) 3 + u p L 2 (Ωp) 3 ≤ C min(C -1 P ν, α min ) f L 2 (Ω) 3 , (3.4)
which yields the first estimate. The second estimates is easily obtained thanks to Poincaré inequality (2.4).

For the last estimate, we infer from a Cauchy-Schwarz inequality and (3.4) that

3 j=1 Ωp νe -τ α |∇u s,j | 2 dx ≤ C min(C -1 P ν, α min ) ( f L 2 (Ω) 3 ) 2 . (3.5)
Multiplying (3.5) by e -τ αmin and using that α min ≤ α(x) for any x ∈ Ω p , we obtain

ν e -τ α ∇u p 2 L 2 (Ωp) 3 = ν 3 j=1 Ωp e -2τ α |∇u p,j | 2 dx ≤ νe -τ αmin Ωp e -τ α |∇u p,j | 2 dx (3.6) ≤ ν Ce -τ αmin νmin(C -1 P ν, α min ) f 2 L 2 (Ω) 3
and the bounds for the velocity are therefore complete. To get the bound for the pressure (p p , p f ) ∈ L 2 (Ω p ) × L 2 (Ω f ), we use the following formula which comes from the weak formulation (2.5) together with the transmission condition (3.2)

b(v, p) = - Ωp p p divv dx - Ω f p f divv dx = Ω f (x) • v(x)dx -a p (u p , v p ) -N p (u p , u p , v p ) -a f (u f , v f ) -N f (u f , u f , v f ), where v ∈ X and (v p , v f ) = (v| Ωp , v| Ω f ).
Following the proof of Lemma 2.2 and using that e -2τ α = (e -τ α ) 2 we obtain the next refined estimates

|N p (u p , u p , v p )| ≤ C N (e -τ α u p ) 2 H 1 (Ωp) 3 v p H 1 (Ωp) 3 , |a p (u p , v p )| ≤ e -τ α ∇u p L 2 (Ωp) 3 ∇v p L 2 (Ωp) (3.7) + α L ∞ (Ω) u p L 2 (Ωp) 3 ) ∇v p L 2 (Ωp) 3 .
Since α| Ωp ∈ Lip(Ω p ), it is almost everywhere differentiable with (essentially) bounded derivatives thanks to Rademacher Theorem. One thus have ∇e -τ α = -τ e -τ α ∇α on Ω p . Using now that 0 ≤ e -τ α(x) ≤ e -τ λmin for almost every x ∈ Ω, we obtain

(e -τ α u p ) H 1 (Ωp) 3 = ∇(e -τ α u p ) L 2 (Ωp) 3 + e -τ α u p L 2 (Ωp) 3 (3.8) ≤ e -τ α ∇u p L 2 (Ωp) 3 + 1 + τ ∇α L ∞ (Ωp) e -τ α u p L 2 (Ωp) 3 ≤ e -τ α ∇u p L 2 (Ωp) 3 + e -τ λmin 1 + τ ∇α L ∞ (Ωp) u p L 2 (Ωp) 3 .
Using now bounds (3.6-3.4) and that τ > 0 is large enough, we infer

|N p (u p , u p , v p )| ≤ C e -τ λmin/2 1 + τ ∇α L ∞ (Ωp) 2 f L 2 (Ω) 3 2 v p H 1 (Ωp) 3 ≤ Cτ 2 e -τ λmin f L 2 (Ω) 3 2 v p H 1 (Ωp) 3 (3.9)
where C > 0 is a generic constant that does not depend on τ nor on f . The upper bound (3.9), the inf-sup condition (2.9) together with (3.6-3.4-3.7) then give the existence of B > 0 that does not depend on τ such that

p p L 2 (Ωp) + p f L 2 (Ω f ) = p L 2 (Ω) ≤ β -1 sup v∈X b(v, q) v X ≤ B,
where (p| Ωp , p| Ω f ) = (p p , p f ).

3.2. Characterization of the limiting problem for (u τ , p τ ). From Theorem 3.1, we infer that the sequence ((u p,τ , p p,τ , u f,τ , p f,τ )) τ is bounded, uniformly with respect to τ in

L 2 (Ω p ) 3 × L 2 (Ω p ) × H 1 (Ω f ) 3 × L 2 (Ω f ) and thus there exists subsequence and (u p , p p , u f , p f ) ∈ L 2 (Ω p ) 3 × L 2 (Ω p ) × H 1 (Ω f ) × L 2 (Ω f ) such that u s,τ → u p weakly in L 2 (Ω p ) 3 , u f,τ → u f weakly in H 1 (Ω f ) 3 and strongly in L 2 (Ω f ) 3 , (3.10) p p,τ → p p weakly in L 2 (Ω p ) and p f,τ → p f weakly in L 2 (Ω f ), e -τ α ∇u s,τ → 0 strongly in L 2 (Ω p ).
The limiting problem satisfied by (u p , p p , u f , p f ) is given in the next result. Theorem 3.2. The distributions

(u p , p p , u f , p f ) ∈ L 2 (Ω p ) 3 × L 2 (Ω f ) × H 1 (Ω f ) 3 × L 2 (Ω f
) defined as the weak limits of (u τ , p τ ) satisfying Problem (3.1,3.2) satisfy the following equation

                               ∇p p + αu p = f in D (Ω p , R 3 ), div(u p ) = 0 in D (Ω p ), -div (ν∇u f -p f I 3 ) + (u f • ∇)u f = f in D (Ω f , R 3 ), div(u f ) = 0 in D (Ω f ), u f = 0 in H 1/2 (Γ 1 ∩ ∂Ω f , R 3 ), n • u p = 0 in H -1/2 (Γ 1 ∩ ∂Ω p ) u f • n pf = u p • n pf in H -1/2 (Γ p,f ), ν∂ n u f -(p f + 1 2 |u f | 2 )n = 0 in H -1/2 (Γ 2 ∩ ∂Ω f , R 3 ), -p p = 0 in H 1/2 (Γ 2 ∩ ∂Ω p ), ν∇u f -p f I 3 -1 2 |u f | 2 I 3 n pf = -p p n pf in W 1-1/q,q (Γ p,f , R 3 ) , (3.11
)

where q = 2p/(2 -p) for any p ∈ [1, 3/2].
Proof. The four first equations of (3.11) follow from (3.10) by passing to the limit in the variational formulation (3.3). To pass to the limit in the non-linear term N f , we use the continuity of N f together with strong L 2 and weak H 1 convergence of (u τ,f ) τ . For the non-linear term N p , it goes to zero thanks to (3.9). The Dirichlet boundary condition on Γ 1 is obtained thanks to the compactness of the trace operator. The boundary condition on Γ 2 ∩∂Ω f and on Γ 2 ∩∂Ω p are easily recovered by integration by part in the variational formulation.

To recover the transmissions condition on Γ p,f and the boundary conditions on Γ 1 ∩∂Ω p , we first introduce the reflexive Banach space L p div (Ω, Hom(R 3 )) which is defined by

L p div (Ω, Hom(R 3 )) = G ∈ L p (Ω, Hom(R 3 )) | div G ∈ L p (Ω) 3 .
Note that any v ∈ W 1,q (Ω, R 3 ), for q = 1/(1 -p -1 ), has a trace v| ∂Ω that belong to W 1-1/q,q (Ω). Therefore, for any (G, v) ∈ L p div (Ω, Hom(R 3 )) × W 1,q (Ω, R 3 ), the following Green's identity holds

Ω v • div Gdx + 3 j=1 Ω G j • ∇v j dx = 3 j=1 G j • n, v , (3.12)
where ., . is the duality product between W 1-1/q,q (Ω) and W 1-1/q,q (Ω). To get the matching of the normal trace on Γ p,f , one only has to remark that u τ ∈ L 2 div (Ω, R 3 ) with uniform bound with respect to τ since div u τ = 0. Therefore, one has u τ,f • n pf = u τ,s • n pf and, using the Green formula, we then get that

u p • n pf = u f • n pf in H -1/2 (Γ 2 ). Similarly, since u p,τ = 0 in H 1/2 Γ 1 ∩ ∂Ω p , R 3 , we end up with u p • n = 0 in H -1/2 (Γ 1 ∩ ∂Ω p ).
We now consider the distribution F τ defined below

F τ = νe -τ α ∇u τ -p τ I 3 - 1 2 e -2τ α |u τ | 2 I 3 .
Theorem 3.1 and the Sobolev embedding (2.7) ensure that F τ ∈ L 2 (Ω, Hom(R 3 )) and that there exists a positive constant C > 0 that does not depend on τ such that

F τ L 2 (Ω,Hom(R 3 )) ≤ C.
We now turn our attention on proving that div(F τ ) is also bounded uniformly with respect to τ . Starting from Equation (2.3), Theorem 3.1 and Holder inequality, we infer that

divF τ L p (Ω) 3 = -f + αu τ -e -2τ α u τ × (∇ × u τ ) L p (Ω) 3 ≤ f L 2 (Ω) 3 + α L ∞ (Ω) u τ L 2 (Ω) 3 + C e -τ α u τ L q (Ω) 3 e -τ α ∇u τ L 2 (Ω) 3 ,
where C > 0 does not depend on τ , q = 2p/(2 -p) and 1 ≤ p ≤ p max for some p max that is going to be specified below. Theorem 3.1 shows that e -τ α ∇u τ L 2 (Ω) 3 is bounded uniformly with respect to τ . It only remains to show that e -τ α u τ L q (Ω) 3 is also bounded independently of τ . To prove this, we recall the continuous embedding H

e -τ α u τ L q (Ω) ≤ u f,τ L q (Ω f ) 3 + e -τ α u p,τ L q (Ωp) 3 ≤ C ∇u f,τ H 1 (Ω f ) 3 + ∇ e -τ α u p,τ H 1 (Ωp) 3 ≤ M,
where C, M > 0 does not depend on τ . Note that, thanks to the Sobolev embedding, we require q ≤ 6 which yield p max = 3/2. This finally prove that F τ is uniformly bounded in L p div (Ω, Hom(R 3 )) and thus there exists a subsequence (still denoted F τ ) which converges weakly toward some F ∈ L p div (Ω, Hom(R 3 )). From (3.10), we get that F τ converge weakly in L 2 (Ω, R 3 ) toward

F = ν∇u f -p f I 3 -1 2 |u f | 2 I 3 in Ω f ∇p p in Ω p .
The uniqueness of the weak limit and the Green's formula (3.12) then shows that 3 j=1 F j n weakly converge to 3 j=1 F j n in W 1-1/q,q (Γ p,f , R 3 ) which complete the proof. Remark 3.3 (Two dimensional case). The penalized model in a two-dimensional setting is given in Remark 2.1. One can see that all the proofs of the Theorem 2.3, 3.1, 3.2 remain the same using that

u 2 ∂ 2 u 1 -u 2 ∂ 1 u 2 u 1 ∂ 1 u 2 -u 1 ∂ 2 u 1 • u = 0.
Also, the limiting Problem (3.11) is exactly the same.

Note that the first equation of Problem (3.11) gives p p ∈ H 1 (Ω p ) since u p ∈ L 2 (Ω p ). Also, since u f ∈ X f , the second equation is valid in L 2 (Ω f ). Now, using that u p = -α -1 ∇p, we obtain the same interface condition on Γ p,f that have been used in [START_REF] Girault | DG approximation of coupled Navier-Stokes and Darcy equations by Beaver-Joseph-Saffman interface condition[END_REF]. Keeping both the velocity and the pressure in the porous medium, they read:

u f • n pf = u p • n pf , ν∇u f -p f I 3 - 1 2 |u f | 2 I 3 n pf = -p p n pf .
Let (t 1 , t 2 ) be any orthonormal basis of the tangent plane to Γ p,f so that one can write v = v n pf n pf + v 1 t 1 + v 2 t 2 for any vector field. Projecting the interface condition according to this decomposition, one gets:

ν∇u f n pf • n pf -p f I 3 -1 2 |u f | 2 I 3 = -p p on Γ p,f (ν∇u f n pf • t j ) t j = 0, j = 1, 2, on Γ p,f . (3.13)
The variational formulation to the limit problem can be directly obtained from (3.3) by letting τ → +∞ and reads

F ind (u f , p f , u p , p p ) ∈ X f × L 2 (Ω f ) × X p × L 2 (Ω p ) such that a f (u f , v f ) + b f (v f , p f ) + Ωp αu p • v p + N f (u f , u f , v f ) + b p (v p , p p ) - Γ p,f (p p -p f ) n pf • (v p -v f )dσ = Ω f f (x) • v f (x)dx + Ωp f (x) • v p (x)dx b f (u f , q f ) + b p (u p , q p ) = 0, ∀(v, v p ) ∈ X f × X p ∀(q f , q p ) ∈ L 2 (Ω f ) × L 2 (Ω p ). (3.14)
The latter can be used to get some a priori estimate as well as the existence and uniqueness of (u f , p f , u p , p p ). Note also that the well-posedness of Problem (3.11)-(3.13) have been studied in [START_REF] Discacciati | Navier-Stokes/Darcy coupling: modeling, analysis, and numerical approximation[END_REF].

Remark 3.4 (Link with Beaver-Joseph-Saffman condition). The transmission condition on the tangential part of ν∇u f n pf is mentioned in [START_REF] Discacciati | Navier-Stokes/Darcy coupling: modeling, analysis, and numerical approximation[END_REF] as a simplified form of the so-called Beavers-Joseph-Saffman interface condition which is used when dealing with fluid-porous media interface problem (see e.g. [START_REF] Cao | Coupled Stokes-Darcy model with Beavers-Joseph interface boundary condition[END_REF][START_REF] Girault | DG approximation of coupled Navier-Stokes and Darcy equations by Beaver-Joseph-Saffman interface condition[END_REF]). The latter reads

∂ n pf u f • t j = C BJ √ K (u f -u p ) • t j on Γ p,f , (3.15)
where C BJ is a dimensionless constant depending on the structure of the porous media and K is the permeability of the porous media. Condition (3.15) has been mathematically justified in [START_REF] Jäger | Modeling effective interface laws for transport phenomena between an unconfined fluid and a porous medium using homogenization[END_REF] using two-scale convergence.

In this paper, we derive the interface fluid-porous media problem through a singular perturbation approach where the porous media is modeled by adding the term αu in the Navier-Stokes equation, also called a penalization model [START_REF] Cimolin | Navier-Stokes/Forchheimer models for filtration through porous media[END_REF]. This approach allow to recover the interface fluid-porous media problem when α = K/η, as η goes to zero, which explain why we only obtained the first order approximation of the Beavers-Joseph-Saffman condition (3.15) in the limiting problem (3.11).

Remark 3.5 (On the computation of the order of convergence). The strongly penalized Darcy-Navier-Stokes problem (2.3) is a singular perturbation problem. Therefore, we cannot expect to have a regular expansion of (u τ , p τ ) as τ → +∞ that is valid everywhere on Ω since boundary layers occurs at the interface between the fluid and the porous medium. In the context of matched asymptotic expansion (see e.g. [START_REF] Cousteix | Asymptotic analysis and boundary layers[END_REF]) the solution to the limit problem (3.11) can be view as the zeroth order exterior expansion and is used to approximate the solution outside the boundary layer. Similar singular perturbation problem involving Stokes-Brinkman interface model has been studied in [START_REF] Angot | Asymptotic study for StokesBrinkman model with jump embedded transmission conditions[END_REF] where the author compute boundary layer corrector as well as the complete asymptotic expansion. Nevertheless, this paper aims at using the model (2.3) for topology optimization to be able to look for optimal design that are porous medium. The complete asymptotic expansion of (u τ , p τ ), as τ → +∞, is thus not needed in the present paper.

Direct extensions.

We present below several direct extension that can be drawn, without additional proof, from the above analysis.

Tensorial porosity. One can easily consider α ∈ Lip(Ω p , Hom(R d )) for d = 2, 3. The coercivity assumption over the porous part of the domain now reads a.e. x ∈ Ω p , α(x) t = α(x), α min I d α(x), for α min > 0. The latter means that the permeability tensor is symmetric and comparable to the identity over the solid part of the domain. Note that the limiting problem (3.11) is exactly the same.

Penalization model using Cauchy-Stress tensor. Since divu = 0, one has ∆u = div(∇u + ∇u T ) = div(2D(u)) where D(u) is the symmetric part of the Jacobian matrix of vector field x ∈ Ω → u(x) ∈ R d . As a result, considering the trace of 2D(u)n is equivalent to considering ∂ n u either in the boundary conditions (see e.g. [START_REF] Court | A fictitious domain approach for the Stokes problem based on the extended finite element method[END_REF] p.4 Remark 1) or in the interface condition and we could have use div (2νe -τ α D(u)) instead of div (νe -τ α ∇u). It is worth noting that the above analysis remains the same thanks to the following formula

Ω div 2νe -τ α D(u) • vdx = ∂Ω 2νe -τ α D(u)n • vdσ -2ν Ω D(u) : D(v)dx,
       -div (νe -τ α ∇u -pI 3 ) + αu = f in Ω, div(u) = 0 in Ω, u = 0 on Γ 1 , νe -τ α ∂ n u -= 0 on Γ 2 .
(3.16) Since we do not have to pass to the limit in the non-linear term, we obtain the following fluid-porous media interface limit problem

                               ∇p p + αu p = f in D (Ω p , R 3 ), div(u p ) = 0 in D (Ω p ), -div (ν∇u f -p f I 3 ) = f in D (Ω f , R 3 ), div(u f ) = 0 in D (Ω f ), u f = 0 in H 1/2 (Γ 1 ∩ ∂Ω f , R 3 ), n • u p = 0 in H -1/2 (Γ 1 ∩ ∂Ω p ), u f • n pf = u p • n pf in H -1/2 (Γ p,f ), ν∂ n u f -p f n = 0 in H -1/2 (Γ 2 ∩ ∂Ω f , R 3 ), -p p = 0 in H 1/2 (Γ 2 ∩ ∂Ω p ), (ν∇u f -p f I 3 ) n pf = -p p n pf in H -1/2 (Γ p,f , R 3 ), (3.17)
where, as seen before, the above system can be written in a stronger sense than in the sense of distribution.

4. Application to topology optimization. We now work with fixed τ > 0 and we consider the following general PDE-constrained optimization problem (P) inf

(u,p,α)∈X×L 2 (Ω)×U J (u, p, α), such that (u, p) satisfy (2.5), α ∈ U ad , (4.1)
where J is a real-valued cost functional. The set U ad ⊂ U is a set of admissible design defined thank to the following sets

Ω f (α) = {x ∈ Ω |α(x) = 0} , Ω p (α) = {x ∈ Ω | α(x) ≥ α min },
where α min > 0. Note that both Ω f (α) and Ω p (α) are Lebesgue mesurable set if α ∈ L ∞ (Ω) since such functions are defined almost everywhere. The set of control is then

U := α ∈ L ∞ (Ω) | Ω = Ω f (α) ∪ Ω p (α) and α ∈ Lip(Ω p ) .
Theorem 2.3 ensures the existence of solution to Problem (2.5) for any α ∈ U and thus the set of constraints is not empty. Nevertheless, Theorem 3.2 which gives the limiting problem demands, for any α ∈ U , that Ω f (α) and Ω p (α) are bounded open set with Lipschitz boundary because the derivation of the limiting problem involves the study of the interface problem (3.1). We emphasize that such a regularity on the level set of α is hard to impose and that U as defined above is clearly not convex. To bypass these difficulties, we propose to work with a fixed decomposition of Ω hence assuming that the location of the fluid/porous media is given. Therefore, the set U we work with is now

U = {α ∈ L ∞ (Ω) | α(x) = 0 a.e x ∈ Ω f , α(x) ≥ α min a.e x ∈ Ω p , α ∈ Lip(Ω p )} ,
where Ω = Ω p ∪ Ω f . Note that the above set is a convex subset of L ∞ (Ω). It is also closed with respect to the strong topology of L ∞ and thus weakly closed. This section is now devoted to prove some existence results for the optimization problem (4.1) and to give next first order optimality conditions. 4.1. Existence of optimal control. We show here the existence of optimal design for Problem (4.1). Recall that for any non-linear function h on a reflexive Banach space and α n a weak-convergent sequence, one has h(α n ) h(α 0 ) if and only if the function is affine. Also, it is classical in PDE-constrained optimization problem to require the state equation to be weakly continuous with respect to the each of its variable [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF][START_REF] Hinze | Optimization with PDE constraints[END_REF]. We then need to consider U ad to be a strong compact of L ∞ (Ω) since we have to pass to the limit in terms of the form e -τ α . We give after the proof some example of set U satisfying this assumption.

Theorem 4.1. Assume that i) J is lower-semicontinous with respect to the (weak, weak, weak * ) topology of

H 1 (Ω) d ×L 2 (Ω)×L ∞ (Ω). ii) ∀(u, p, α) ∈ X × L 2 (Ω) × L ∞ (Ω) × U , one has J (u, p, α) > -∞.
iii) The set U ad is a closed convex bounded subset of U . It is also a strong compact of U for the L ∞ norm. Then the optimization Problem (4.1) has at least one optimal solution.

Proof. The proof is based on standard weak compactness argument for minimizing sequence and on the continuity of the mapping α ∈ U ad → (u, p) ∈ X × L 2 (Ω) where (u, p) is a solution to Problem (2.5). Let (u n , p n , α n ) be a minimizing sequence for Problem (4.1). Since α n ∈ U ad , Theorem 2.3 show that there exists some (u n , p n ) ∈ X × L 2 (Ω) satisfying Equation (2.5) such that u n H 1 (Ω) d and p n L 2 (Ω) are both uniformly bounded with respect to n. One then has a subsequence of (u n , p n ) which converges weakly in H 1 (Ω) d × L 2 (Ω) toward some (u 0 , p 0 ) and, thanks to (2.7), u n → u 0 strongly in L p (Ω) d for 1 ≤ p ≤ 4. From iii), we infer that α n converges strongly toward α 0 ∈ U ad in L ∞ (Ω). Owning to this strong convergence in L ∞ (Ω), we get that e -τ αn → e -τ α0 strongly in L ∞ (Ω) and, passing to the limit in the variational formulation (2.5), we obtain that (u 0 , p 0 ) is a solution to Problem (2.5) with α = α 0 . Assumptions i), ii) then give

J (u 0 , p 0 , α 0 ) ≤ lim inf n→+∞ J (u n , p n , α n ) = inf (P) ,
which shows that (u 0 , p 0 , α 0 ) is a solution to Problem (4.1).

We now give some examples where the control set U ad is a strong compact of L ∞ . We recall that we assumed a fixed location of fluid/porous media and now prove that iii) holds. Indeed, if α n ∈ U ad is bounded, it is also bounded in Lip(Ω p ) and vanishes on Ω f , Ascoli's theorem then ensures that there exists a subsequence (still denoted α n ) and α 0 ∈ Lip(Ω p ) such that

α n → α 0 strongly in L ∞ (Ω p ), α n → 0 strongly in L ∞ (Ω f ).
This gives the strong convergence of α n toward some ᾱ in L ∞ (Ω) and the proof of Theorem 4.1 remains the same.

Another way to get assumption iii) satisfied is to consider α that belong to some higher order Sobolev spaces that are compactly embedded into C(Ω). Such set U are given (see e.g. [START_REF] Galdi | An introduction to the mathematical theory of the Navier-Stokes equations: Steady-state problems[END_REF]) by the two following Banach spaces

U = W 1,q 0 (O) → → C(O) if q > n W k,q (O) → C r,γ (O) → → C(O) for γ ∈ (0, 1), k -r -γ = d/q, (4.2)
where O ⊂ R d is a bounded open set with Lipschitz boundary and the compactness of the last embedding follows from Ascoli's Theorem. As a result, one can get the required regularity by adding suitable Sobolev norms of α in the objective functional. These additional terms can be seen as regularization terms and read as follow

J (u, p, α) = J 0 (u, p) +        c 1 q α q L q (Ωp) + c 2 q ∇α q L q (Ωp) + c 3 2 ∇α 2 L 2 (Ωp) , c 1 2 α 2 L 2 (Ωp) + c 2 2 ∇α 2 L 2 (Ωp) + c 3 2 d i,j=1 ∂ i ∂ j α 2 L 2 (Ωp) , (4.3)
where c 1 , c 2 , c 3 > 0 and we used r = 0, q = 2, γ = 2 -d/2 for the parameters of the second embedding to have a Sobolev space involving integer derivatives. The main steps of the existence proof then remains the same since, assuming that J 0 (u, p) > -∞, ensure that α is bounded in U and thus has a subsequence that vanishes in Ω f and converges strongly in C(Ω p ).

4.2.

First order optimality condition. The next result gives some general assumptions yielding first order optimality condition for a PDE-constrained optimization problem and can be found in [START_REF] Hinze | Optimization with PDE constraints[END_REF] p. 64, Section 5.2.

Theorem 4.2. Let U, Y, Z be Banach spaces and consider the following general nonlinear problem

min (y,u)∈Y ×U f (y, u) subject to E(y, u) = 0, u ∈ U ad .
Assume that the following hold 1. U ad ⊂ U is non-empty and convex.

2. f :

Y × U → R and E : Y × U → Z are continuously Fréchet differentiable. 3. For all u ∈ V in a neighborhood V ⊂ U of U ad ,
∂ u f (ū), u -ū U * ,U ≥ 0, ∀u ∈ U ad .
The reduced problem associated to (4.1) reduces to find the minimum of J 1 (α) = J (u(α), p(α), α) over the closed convex set U ad . We show below that the mapping α → S(α) = (u(α), p(α)) is C 1 with respect to the topology of the strong convergence and thus, assuming that J is differentiable, one get that the optimum (ū, p, ᾱ) satisfies the following first order optimality condition

J 1 (ᾱ)[δα -ᾱ] ≥ 0, ∀δα ∈ U ad .
Since J 1 (α) = J (S(α), α), the chain rules gives

J 1 (α)[δα] = ∇ u,p J (S(α), α), S (α)[δα] + ∇ α J (S(α), α), δα , (4.4)
where the brackets stand for the duality product. We emphasize that, among the assumptions of Theorem 4.2, only the inversibility of the linearized problem has to be proved.

Derivability of the design to flow mapping. We study first the linearization of the penalized Navier-Stokes-Darcy system. Computing the Fréchet derivative of the operator involved in Problem (2.5) at (u, p) ∈ X × L 2 (Ω), we get that (w, q) is a solution to the linearization of (2.5) if it satisfies

F ind (w, q) ∈ X × L 2 (Ω) such that ∀v ∈ X, a(w, v) + N (u, w, v) + b(v, q) = Ω f (x) • v(x)dx + Γ2 ϕ • vdσ, ∀π ∈ L 2 (Ω), b(w, π) = 0, (4.5)
where

N (u, w, v) = - Ω e -2τ α (u • wdiv(v) + u × (∇ × w) • v + w × (∇ × u) • v) dx.
Integrating by parts then gives the following strong formulation of Problem (4.5)

-div νe -τ α ∇w -qI 3 -e -2τ α u • wI 3 -e -2τ α (u × (∇ × w) + w × (∇ × u)) + αw= f in Ω, div(w)= 0 in Ω, (4.6) w= 0 on Γ 1 , νe -τ α ∂ n w -(q + e -2τ α u • w)n= ϕ on Γ 2 .
The well-posedness of the linear Problem (4.5) can be obtained with standard method and its precise statement depends on the regularity of (u, p) as one can see in the next result.

Theorem 4.3.

Let (u, p) ∈ X × L 2 (Ω) and (f , ϕ) ∈ X × H -1/2 (Γ 2 , R 3 ) be given. We have the following i) Assume that (u, p) ∈ (X ∩ L ∞ (Ω))×L 2 (Ω)
and consider the unbounded operator of X , (L, V) defined by

Lw, v = Aw + N w, v = a(w, v) + N (u, w, v), ∀v ∈ V.
Then the spectrum of (L, V) is discrete with no accumulation point and if 0 / ∈ σ(L), then Problem (4.6) is well-posed for any ν > 0. ii) If ν > 0 is only large enough, then Problem (4.6) is well-posed. Proof. From the inf-sup condition (2.9), the pressure is determined as soon as the velocity is. Also, the unicity can be obtained using standard technique since the underlying problem is linear. Note that the velocity is a zero of the linear mapping Φ : V → V defined by

Φ(w), v := a(w, v) + N (u, w, v) - Ω f (x) • v(x)dx - Γ2 ϕ • vdσ. We detail below both cases. i) Assume that (u, p) ∈ (X ∩ L ∞ (Ω)) × L 2 (Ω).
We first prove that the spectrum of the unbounded operator (L, V) is discrete. Using 2ab ≤ a 2 t + b 2 t -1 for any t > 0, we obtain Lw, w = a(w, w) + N (u, w, w)

≥ νC 1 (α min , Ω) u 2 H 1 (Ω) - Ω e -2τ α u × (∇ × w) • wdx ≥ (νC 1 (α min , Ω) -u L ∞ (Ω) t 2 ) w 2 H 1 (Ω) -u L ∞ (Ω) 1 2t w 2 L 2 (Ω) 3 ,
where t > 0 is small enough to ensure that

(νC 1 (α min , Ω) -u L ∞ (Ω) t/2) > νC 1 (α min , Ω)/2.
Therefore, using Lax-Milgram lemma, the operator (L + λI) is invertible with a continuous inverse if λ > 0 is large enough. This shows that the resolvent set of (L, V) is non-empty. Since the embedding V → X is compact, one has that the resolvent of (L, V) is compact and thus its spectrum is discrete with no accumulation point. Remark now that if 0 / ∈ σ(L) then, for all ν > 0, there exists a unique w ∈ V such that Lw = f for any f ∈ X . This gives Φ(w), v = 0 for all v ∈ V and thus the existence and uniqueness of w.

ii) Due to the lack of regularity of (u, p), we cannot get a Garding inequality for the operator L as above. Instead, to get a bound on N (u, w, w), we split Ω using that Ω = Ω f ∪ Ω p and apply Holder inequality

Ω e -2τ α u × (∇ × w) • wdx ≤ C(Ω) e -αmin ∇w p 2 L 2 (Ωp) e -τ α ∇u p L 2 (Ωp) + C(Ω) ∇w f 2 L 2 (Ω f ) ∇u f L 2 (Ω f ) .
Using Young inequality, we then infer that Φ satisfies a Garding inequality if ν > 0 is large enough owning to Theorem 3.1. The existence is then proved as in ii) and the unicity holds for large enough ν > 0.

Remark 4.4. Theorem 4.3 giving the well-posedness of the linearized Problem (4.5) can be extended for non vanishing divergence as well (see e.g. [START_REF] Amrouche | Lp-theory for vector potentials and Sobolev's inequalities for vector fields: Application to the stokes equations with pressure boundary conditions[END_REF][START_REF] Raymond | Stokes and NavierStokes equations with nonhomogeneous boundary conditions[END_REF] for Stokes and Navier-Stokes equations). To handle such case, consider the following orthogonal decomposition

X = V ⊕ V ⊥ , (4.7)
and the associated projectors P 0 : X → V, P ⊥ : X → V ⊥ . Projecting Problem (4.5) according to the decomposition (4.7) then gives the following equivalent system

b(P ⊥ w, q) = (g, q) L 2 (Ω) , ∀q ∈ L 2 (Ω), a(P 0 w, v 0 ) + N (u, P 0 w, v 0 ) = f , v 0 , ∀v 0 ∈ V, b(v, p) + a(P 0 w + P ⊥ w, v) + N (u, P 0 w + P ⊥ w, v) = f , v , ∀v ∈ X, (4.8)
Note that the second equation of (4.8) is well-posed as soon as any assumption of Theorem 4.3 hold. This gives the existence and uniqueness of P 0 w. Now, let B : X → L 2 (Ω) be the operator defined by Bu, q = b(u, q). The inf-sup condition (2.9) with Lemma 4.1 p.58 [START_REF] Girault | Finite element methods for Navier-Stokes equations: theory and algorithms[END_REF] then ensures that the first and last equations of (4.8) are well-posed and thus the linearized problem (4.5) with non-homogeneous divergence is well-posed under the same assumptions as those of Theorem 4.3.

We can now apply the implicit function theorem to get some smoothness of the design-to-flow mapping. First order optimality condition and adjoint model. We now use Theorems 4.2,4.3,4.5 to get the first order optimality condition involving the adjoint model. This technique is standard and can be found for instance in [START_REF] Hinze | Optimization with PDE constraints[END_REF] Corollary 5.2.4 p. 66.

Recall first that S(α) = (u(α), p(α)) satisfying Problem (2.5) are also solution to F (u, p, α) = 0 where

F : X × L 2 (Ω) × U ad → X × L 2 (Ω) is defined by F (u, p, α), (v, q) = a(u, v) + N (u, u, v) + b(v, p) + b(u, q) - Ω f (x) • v(x)dx - Γ2 ϕ • vdσ.
The implicit function theorem then gives that

S (α)[δα] = -∂ (u,p) F (u(α), p(α), α) -1 • ∂ α F (u(α), p(α), α)[δα], (4.9)
where ∂ . is the Fréchet derivative with respect to the subscripted variable. Using now (4.4), we obtain the derivative of J 1 (α) = J (S(α), α) at any α ∈ U ad where the . - * denotes the adjoint of the inverse of the linear operator ∂ (u,p) F (S(α), α). Let (ū, p, ᾱ) ∈ X × L 2 (Ω) × U ad be a solution to the optimal control problem (4.1). Since U ad is convex, one has that optimal solution satisfies

J 1 (ᾱ)[δα -ᾱ] ≥ 0, ∀δα ∈ U ad , (4.11)
where the derivative of J 1 is given in (4.10) and involves the resolution of the so-called adjoint problem

F ind (w, q) ∈ X × L 2 (Ω) such that ∀v ∈ X, a(v, w) + N (u, v, w) + b(v, q) = Ω f (x) • v(x)dx + Γ2 ϕ • vdσ, ∀π ∈ L 2 (Ω), b(w, π) = -(g a , π) L 2 (Ω) , (4.12)
which is the variational formulation to the Problem

-div νe -τ α ∇w -qI 3 + ∇ × e -2τ α u × w -e -2τ α (udiv w + ∇ × u × w) + αw= f a in Ω, div(w)= g a in Ω, (4.13) w= 0 on Γ 1 , νe -τ α ∂ n w -(qn + e -2τ α n × (u × w))= ϕ a on Γ 2 .
The well-posedness of the adjoint equation (4.12) can be proved similarly as Theorem 4.3. For nonhomogeneous divergence, this can be done as in Remark 4. [START_REF] Andreasen | Topology optimization of microfluidic mixers[END_REF]. Now assume that the derivative of the cost functional can be written as

∂ u J (u, p, α), δu = ∂ u J Ω (u, p, α), δu X * ×X + ∂ u J Γ (u, p, α), δu H -1/2 ×H 1/2 , ∂ p J (u, p, α), δp = (∂ p J (u, p, α), δp) L 2 (Ω) ,
where (∂ u J Ω (ū, p, ᾱ), ∂ u J Γ (ū, p, ᾱ)) ∈ X * × H -1/2 (Γ 2 , R 3 ). The first order optimality condition (4.11) can finally be recast as follow (ū, p) ∈ X × L 2 (Ω) satisfy (2.5), ( w, q) ∈ X × L 2 (Ω) satisfy (4.13) with div w = g a ,

where g a = ∂ p J (ū, p, ᾱ) and (f a , ϕ a ) = (∂ u J Ω (ū, p, ᾱ), ∂ u J Γ (ū, p, ᾱ)) ,

(∂ α J (ū, p, ᾱ) -∂ α a(ū, w) -∂ α N (u, u, w)) [ᾱ -δα] ≥ 0, ∀δα ∈ U ad , with ∂ α a(u, w)[δα] = ν d j=1 Ω (-τ δα)e -τ α ∇u j • ∇v j dx + Ω δαu • vdx ∂ α N (u, u, w)[δα] = τ Ω e -2τ α |u| 2 div(w) + 2e -2τ α u × (∇ × u) • w δαdx.
We finally emphasize that the first order optimality conditions (4.14) correspond to the standard adjoint method (see e.g. [START_REF] Hinze | Optimization with PDE constraints[END_REF] Corollary 5.2.4). Remark 4.6 (Derivatives of J with regularization terms). We give here the derivative of the cost functional (4.3) involving regularization terms which ensures that α belong to a strong compact of L ∞ (Ω). In the case where J 1 (α) = c1 q α q L q (Ωp) + c2 q ∇α q L q (Ωp) + c3 2 ∇α 2 L 2 (Ωp) for α ∈ W 1,q 0 (Ω p ) with q > d, one has ∀δα ∈ W 1,q 0 (Ω p ),

∂ α J 1 (α)[δα] = Ωp c 1 |α| q-2 αδα + c 2 |∇α| q-2 ∇α • ∇δα + c 3 ∇α • ∇δαdx = - Ωp div((c 3 + c 2 |∇α| q-2 )∇α) -c 1 |α| q-2 α δα dx.
Therefore we end up with solving a nonlinear coercive elliptic problem involving the q-Laplacian.

In the case

J 1 (α) = c1 2 α 2 L 2 (Ωp) + c2 2 ∇α 2 L 2 (Ωp) + c3 2 d i,j=1 ∂ i ∂ j α 2 L 2 (Ωp) with α ∈ H 2 0 (Ω p ) we get ∀δα ∈ H 2 0 (Ω p ), ∂ α J 1 (α)[δα] = Ωp c 1 α -c 2 ∆α + c 3 ∆ 2 α δα.
We then need to solve a fourth order problem.

5. Numerical illustration: Velocity tracking. In this section, we wish to give an illustration of the penalized Navier-Stokes-Darcy model applied to topology optimization. We would like to emphasize that the goal of this section is to show some numerical illustration and thus we do not intend to prove that all the algorithms and discretization used converges. Nevertheless, the convergence of finite element method applied to the penalized Stokes-Darcy model can be deduced from [START_REF] Girault | Finite element methods for Navier-Stokes equations: theory and algorithms[END_REF] since both problems have a coercive principal part. Having this in mind, we are going to use FreeFEM ++ [START_REF] Hecht | New development in FreeFem++[END_REF].

Let Ω = [0, 1] × [0, 2] be an enclosure where air flows comes from the left of the room with a given velocity u 0 = (1, 0). The velocity and pressure of air in the room then satisfy a steady state Stokes equation with inhomogeneous Dirichlet boundary condition. We are interested in finding the location of a Darcy porous medium with porosity α such that the following velocity-tracking cost functional is minimized

J (α) = 1 2 Ω |u -u d | 2 dx
where u d = (1, 0), u b = (1, 0). We chose to take into account the porous media with the following penalized Stokes-Darcy model

           -div (νe -τ α ∇u -pI 3 ) + αu = 0 in Ω, div(u) = 0 in Ω, u = u 0 on {0} × [0, 1], u = 0 on [0, 2] × ({0} ∪ {1}) , νe -τ α ∂ n u -pn = 0 on {2} × [0, 1].
(5.1)

The velocity-tracking optimisation problem is now given as min (u,p,α)∈X×L 2 (Ω)×U J (α), (u, p) satisfies (5.1), α ∈ U ad ,

where the set U ad is defined below

U ad = {α ∈ L ∞ (Ω) | α(x) = 0 or α(x) ≥ α min } .
We are going to need the adjoint equation associated to problem (5.2). The latter can be obtained from (4.13) and reads here

           -div (νe -τ α ∇v -qI 3 ) + αv = f in Ω, div(v) = 0 in Ω, v = 0 on {0} × [0, 1], v = 0 on [0, 2] × ({0} ∪ {1}) , νe -τ α ∂ n v -qn = 0 on {2} × [0, 1], (5.3) where f = -(u -u d ).
The variational formulation of the direct and adjoint problems can be respectively found in Equations (2.5) and (4.12) setting N = N = 0. The finite element discretisation is done with a regular mesh composed of triangles build from a uniform n × n grid. Both direct and adjoint states are computed with P 2 finite element for the velocity fields and P 1 for pressure. The porosity fields α is sought in P 0 and is thus a piecewise constant function.

The optimization problem (5.2) is now solved with a steepest descent algorithm. The derivative of J with respect to α can be computed thanks to (4.10) (see also (4.14)) which is actually the standard adjoint algorithm for PDE-constrained optimization problem [START_REF] Herzog | Algorithms for PDE-constrained optimization[END_REF][START_REF] Petra | Model variational inverse problems governed by partial differential equations[END_REF].

Algorithm 1:

1. α = 0 ∈ P 0 is given. 2. Compute (u h , p h ) ∈ P 2 × P 1 satisfying Problem (5.1).

3. Compute (v h , q h ) ∈ P 2 × P 1 satisfying the adjoint problem (5.3). 4. Update the porosity thanks to the optimality condition (4.14)

β = max    0, α -λΠ h   -ντ e -τ α d j=1 (∇u j,h • ∇v j,h ) + u h • v h      , α = max {α min , β} χ (β>0)
where Π h : L 2 (Ω) → P 0 is the finite element projector and the max has to be understood on each cell of the mesh. Also, χ (β>0) is defined on each cell of the mesh and return zero if β = 0 and the value of β if β > 0. 5. Return to Step 2. Remark 5.1. We use the Crout solver to solve both direct and adjoint problem. The latter needs every sub-matrices to be invertible and we thus add the term επ, π = p, q, in the incompressibility condition of both direct and adjoint problems. The convergence, as ε → 0, of the solution to the discrete perturbed problems to the solution of the discrete original one is of the order of ε aas proved in [START_REF] Roberts | Mixed and Hybrid Methods, Handbook of Numerical Anaysis[END_REF]. The analysis can be extended for our penalized Navier-Stokes-Darcy model since both problem share the same mathematical structure.

The following parameters are used in the numerical simulations α min = 10, ν = 1, τ = 10, n = 100, λ = 10, n iter = 1500, ε = 0.000001.

We chose ν = 1 since we work with Stokes equation that can be used as a model for fluid flow at low Reynolds number. We emphasize that we run Algorithm 1 with 1500 iterations since, as we show below in Figure 5.1 the values of the cost functional are stationary, when rounding up to four digits, before reaching the end of iterates. Values of the cost function without porous medium embedded in Ω is J = 0.1801 and after 1500 iterations is J = 0.0043 hence reducing the initial value of the cost function by a factor 42. We represent in Figure 5.2 the velocity field when no porosity is embedded in Ω and in Figure 5.3 after 1500 iteration of the steepest descent algorithm. The optimal porosity design can be found in Figure 5.4. Note that most of the porous media is added at the outlet of the room. This can be explained by the fact that we impose u 0 = (1, 0) = u d at the entrance. Finally, we represent the term exp(-τ α) in Figure 5.4. This clearly shows where the the porous media is located. Note that, when α > 0, exp(-τ α) numerically vanished and so the penalized model (5.1) is actually the Darcy model for porous media. As a result, we have the Stokes model inside the blue zone and the Darcy model inside the red zones.

6. Conclusions and Perspectives. This paper proposed a penalization model for the Navier-Stokes-Darcy system and its application to topology optimization problems. The main advantage of this approach is that one can use the porosity as optimization parameter and thus find optimal design that are porous medium. We first prove that we recover the standard Navier-Stokes-Darcy model as the weak-limit of our penalization model. We then gave necessary condition on both the control set and the cost functional to have at least one optimal solution. We also derive the first order optimality condition and finally end this paper with some numerical simulations, on Stokes flow, to show the interest of the approach.

A first perspective is to perform more intensive numerical simulation of the penalized Navier-Stokes-Darcy model applied to topology optimization. Indeed, we used in this paper one of the simplest optimization algorithm (steepest descent) and this could be worth using, for instance, Newton method. Also, we emphasize that, at the continuous level, one needs the porosity to belong to a strong compact of L ∞ to get at least on optimal design and we thus propose to add some regularization terms (see Remark 4.6) that were not considered in the numerical simulation. The impact of those terms on the optimal design could be studied as well. Another perspective could be to consider other methods for solving both primal and adjoint problem. For instance, domain decomposition methods should be more suitable to handle the heterogeneous nature of the fluid-porous media transmission problem. Error estimates for the numerical approximation of the optimality condition could also be investigated for finite element or finite volume method. Such analysis has already been carried out for finite element approximation of distributed control problem for Navier-Stokes equation [START_REF] Casas | Error estimates for the numerical approximation of a distributed control problem for the steady-state Navier-Stokes equations[END_REF][START_REF] Gunzburger | Analysis and finite element approximation of optimal control problems for the stationary Navier-Stokes equations with distributed and Neumann controls[END_REF]. Convergence of finite volume method has been proved in [START_REF] Evgrafov | Convergence of cell based finite volume discretizations for problems of control in the conduction coefficients[END_REF] for a control problem in linear conduction coefficient. All of these previous works can serve as a basis to handle the case of topology optimization with our penalized model. Another possible future work is linked to the topological asymptotic expansion [START_REF] Amstutz | The topological asymptotic for the Navier-Stokes equations[END_REF][START_REF] Amstutz | Topological derivatives for a class of quasilinear elliptic equations[END_REF] which compute the so-called topological derivative as the first order term in the asymptotic expansion of the solution to PDE in domain with hole with vanishing size. An idea could be to replace the hole by a porous medium with a given permeability and compute the asymptotic expansion of the solution to the fluid-porous medium transmission problem as the characteristic size of the porous medium goes to zero. A last direction for future work is related to the derivation of similar penalization model for other modeling of porous media (see e.g. [START_REF] Ehrhardt | An Introduction to Fluid-Porous Interface Coupling[END_REF]). For instance the Brinkman model [START_REF] Angot | Asymptotic study for StokesBrinkman model with jump embedded transmission conditions[END_REF] involves a second order term, with an effective viscosity, in the porous medium and can be used to account high porosity of the porous media. Another example is the Forchheimer model [START_REF] Cimolin | Navier-Stokes/Forchheimer models for filtration through porous media[END_REF] whose model involves an additional non-linear term inside the porous media to deal with microscopic inertial effects. The case of time-dependent model is also an interesting continuation of the present work.

  where A : B = d i,j=1 A ij B ij . Extension to Stokes-Darcy model. The above analysis extends without any change to the case of Stokes-Darcy model. The resulting penalized model is given by

  the state equation E(y, u) = 0 has a unique solution y = y(u) ∈ Y . 4. ∂ y E(y(u), u) ∈ L(Y, Z) has a bounded inverse for all u ∈ U ad . Then, denoting f (u) = f (y(u), u), any local solution ū of the reduced problem min u∈U f (u) such that u ∈ U ad , satisfies the variational inequality

Theorem 4 . 5 .

 45 Let β ∈ U and (u(β), p(β)) ∈ X × L 2 (Ω) be a solution to Problem (2.5) with α = β. Assume that Problem (4.6) is well-posed (see Theorem 4.3). Then there exists an open neighborhood O(u) × O(p) ⊂ X × L 2 (Ω) , an open neighborhood O(β) ⊂ U and a C 1 mappingS : α ∈ O(β) → (u(α), p(α)) ∈ O(u) × O(p),where (u(α), p(α)) satisfy Problem (2.5).

J 1 (

 1 α)[δα] = ∇ (u,p) J (S(α), α), S (α)[δα] + ∇ α J (S(α), α), δα = ∇ α J (S(α), α), δα -(4.10) ∂ (u,p) F (S(α), α) - * ∇ (u,p) J (S(α), α), ∂ α F (S(α), α)[δα] ,

Fig. 5 . 1 .

 51 Fig. 5.1. Value of cost functions over iteration.

Fig. 5 . 2 .

 52 Fig. 5.2. Velocity fields without porosity. Left: u 1 , Right: u 2 .

Fig. 5 . 3 .

 53 Fig. 5.3. Velocity fields after optimization. Left: u 1 , Right: u 2 .

Fig. 5 .

 5 Fig. 5.4. α after 1500 iteration. Left: α, Right: exp(-τ α).

  1 (O) ⊂ L 6 (O) (see e.g. [44] Chap. II p.159) for any locally Lipchitz set O ⊂ R d for d = 2, 3. Using this together with Theorem 3.1 for the H 1 bound on Ω f and (3.8) for the H 1 bound on Ω p , we obtain

Affiliation: Physique et Ingénierie Mathématique pour l' Énergie et l'Environnement (PIMENT), Université de la Réunion,

rue Joseph Wetzell, 97490 Sainte-Clotilde. E-mail: {Alain.Bastide, Pierre -Henri.Cocquet, Delphine.Ramalingom}@univ-reunion.fr. 1