Olivier Heen
email: olivier.heen@inria.fr

Thomas Genet
email: thomas.genet@inria.fr

On the Unobservability of a Trust Relation in Mobile Ad Hoc Networks

 where anonymity and unobservability contribute to a better privacy.

Introduction

These last years we see an explosion of the number of mobile devices featuring wireless capabilities. Such devices can self-organized as the nodes of a Mobile Ad hoc NETwork, a MANET. Some nodes can join and leave the MANET at any time according to their moves and communication range.

The notion of trust naturally appears in MANET: some nodes can establish long term trust relations in order to better communicate the next time they meet. Typically, trusted nodes may exchange private data while the other nodes only exchange public data. The so-called resurrecting duckling [START_REF] Stajano | The resurrecting duckling: Security issues for ad-hoc wireless networks[END_REF] is a well known-solution for the establishment and detection of trust relations between nodes. Other solutions like [START_REF] Prigent | Secure long term communities in ad hoc networks[END_REF] are particularly suitable for home networks. According to many solutions the user contributes to the establishment of the trust relations, and later on the nodes detect their trust status as soon as they start to communicate without any user intervention. By trust status we mean the fact of sharing a trust relation or not.

In the same time, privacy becomes a serious concern in a variety of domains including many typical application fields of MANET: home networks, Personal Area Networks (PAN), vehicular networks (VANET) [START_REF] Toh | Research challenges in intelligent transportation networks[END_REF][START_REF] Dötzer | Privacy Issues in Vehicular Ad Hoc Network[END_REF], etc. The greater the ability of devices to communicate, the greater the risk of disclosing private information. In the field of VANET for instance, malicious observers should not trace the moves of a vehicle for a long period of time, simply based on its network communications [START_REF] Schneier | Blog article, tracking automobiles through their tires[END_REF].

In some cases, even the existence of a trust relation between two nodes must remain private. This can be the case when the nodes belong to a same group of interest, friends, family etc. In the field of VANET, this can be critical for unmarked police cars: they should not be disclosed by their trust relation.

In this paper, we address the following question: is it possible to detect trust relations in MANET in an anonymous and unobservable manner? In particular, can two communicating nodes detect their trust status and still be sure that no untrusted node has detected it either?

We mainly focus on the existence of a solution to the problem, in presence of passive or active attackers. To do so, we use a binary notion or trust: a node does trust another node, or does not. Note that some more graduated notions of trust exist, with for instance a parameter for the intensity of trust or taking in account the number of positive vs. negative experiences with other nodes. . . We do not formally consider these notions. Instead we consider that any non-null trust is a maximal trust, and bring unobservability in this worst case.

The section 2 provides the necessary definitions as well as the constraints we will respect in the design of the solution. The section 3 provides the appropriate notations and gives a full description of our solution. The section 4 gives arguments for the proof of the security, the anonymity and the unobservability properties. A significant part of the verification is performed using the AVISPA [START_REF] Armando | The AVISPA Tool for the automated validation of internet security protocols and applications[END_REF] framework together with the animation tool SPAN [START_REF] Boichut | Using Animation to Improve Formal Specifications of Security Protocols[END_REF].

Preliminaries

Definitions

Firstly, we give the definition related to the nodes and the trust relation.

Node: anodeA is a communicating device with a pair of keys denoted (K A ,K -1 A ). No key in the pair will be public. The node set is noted N . We define the trust relation for which we want to achieve unobservability against some attackers.

Trust relation: an o d eB trusts another node A, denoted B>A,i f B knows the key K A . We can see that B>A does not imply that A>B.

Trust community: the trust community of a node B, denoted T B , is the set of nodes trusted by B.

T B = {A 2N,B>A}
We give in the following the most common definition of the anonymity and some related concepts. These definitions are taken from [START_REF] Pfitzmann | Anonymity, unobservability, and pseudonymity -a proposal for terminology[END_REF][START_REF] Pfitzmann | Anonymity, unobservability, and pseudonymity: A consolidated proposal for terminology[END_REF].

Anonymity set: an anonymity set is a set of nodes having the same attributes and capable to perform the same actions.

Anonymity: a node remains anonymous if we cannot identify it in an anonymity set. Sender (resp. receiver) anonymity is provided when a message analysis do not allow determining its sender (resp. receiver) in the anonymity set.

Unobservability: there is unobservability of an event when an attacker cannot deduce the existence of this event. There is unobservability of the sender (resp. receiver) when an attacker cannot deduce that a message was sent (resp. received).

We will first provide a solution for the unobservability of a trust relation against a passive attacker.

Passive attacker: a passive attacker can dump every message sent on the network. She does not know, ap r i o r i , any cryptographic key.

Then, we will refine this solution to thwart an active attacker.

Active attacker: an active attacker can dump, add, remove and modify every message on the network. We say that "the attacker is the network". Note that this kind of attacker is also called a Dolev-Yao attacker.

Trust relation

Numerous trust models exist, coming from the simple secret key sharing to certificate chains. For further information on the different trust models, read [START_REF] Balfanz | Talking to strangers: Authentication in adhoc wireless networks[END_REF] or [START_REF] Prigent | Secure long term communities in ad hoc networks[END_REF]. The model we use is particularly adapted to the MANET. We give here the main characteristics and some use cases.

Establishment of trust: the trust relation B>A is established by transmitting K A to B over a secure channel. Note that if A changes its key or if B forgets K A it will be necessary to rebuild the trust relation. If B does not trust A anymore it just has to remove K A from its key database.

Extension of trust (a.k.a recommendation): in the case B>A, B can decide to recommend A to another node X; it just has to send K A to X over a secure channel. In the example of vehicular network, a truck B can trust road equipment A (base station, toll, etc. ) and transmit this trust to another truck X,evenifA is not here during the key transmission. Our solution is based on three significant techniques: renewable pseudonyms, message encryption, periodical broadcast.

Our solution

The general format of a message in our protocol is:

{i.0.db} K -1 A .{j.0.db} K -1 B .g i a.m where {i.0.db} K -1 A is a pseudonym of the sender, {j.0.db} K -1 B
is a pseudonym of the receiver, g i a Diffie-Hellman value, m the encrypted payload.

Note that since an observer should not distinguish anything, all the messages must respect the same format.

Our solution builds a Diffie-Hellman key and tries to authenticate this key by using the information contained in the pseudonyms. The authentication will succeed only when a trust relation exists, otherwise the nodes know that they do not share any trust relation (refer to [START_REF] Blake-Wilson | Authenticated Diffie-Hellman key agreement protocols[END_REF] for more information about authenticated Diffie-Hellman exchanges).

We now provide the sequence of message according to our solution in various cases.

Case B>A (B knows K A ) and A starts:

1. A ! All : {i.0.db} K -1 A .R.g i .m 1 (with R some padding) 2. B checks that {{i.0.db} K -1 A } K A = i 0 .0.db and g i 0 = g i 3. B ! All : {i.j.db} K A .{i.0.db} K -1 A .g j .m 2 4. A checks that {{i.j.db} K A } K -1 A = i 0 .j 0 .db and i 0 = i and g j 0 = g j 5. A ! All : {i.0.db} K -1 A .{i.j.db} K A .g k .{Public} g ij 6. B ! All : {i.j.db} K A .{i.0.db} K -1 A .g l .{Private} g ij At the step 2, B detects that {i.0.db} K -1
A is a pseudonym of a trusted node. A simple manner for detecting this is to sequentially3 try all the keys in T B until the key K A successfully decrypts the pseudonym of A.

At the step 4, A detects that {i.j.db} K A is a pseudonym of a trusted node, nevertheless A does not know that this node is B.

The nodes A and B may then continue to use there pseudonym in the further exchanges. They can send as many encrypted messages as they need, with A sending encrypted Public messages and B sending Public or Private messages:

A ! All : {i.0.db} K -1 A .{i.j.db} K A .g l .{Public} g ij B ! All : {i.j.db} K A .{i.0.db} K -1 A .g m .{Private} g ij
If one node changes its pseudonym, a new detection of trust will happen and re-establish the communication.

Case without trust relation: We describe in this section the sequence of exchanged message between two nodes A and C that do not share any trust relation. Nevertheless, A and B want to communicate with each other.

1. A ! All : {i.0.db} K -1 A .R.g i .m 1 (with R some padding) 2. B ! All : {j.0.db} K -1 B .{i.0.db} K -1 A .g j .m 2 3. A ! All : {i.0.db} K -1 A .{j.0.db} K -1 B .g k .{Public} g ij
In this case, A and B do not share any key and then they cannot collaborate for checking the Diffie-Hellman values. The key used to encrypt the payload is the Diffie-Hellman key g ij generated with the third parts of the messages 1 and 2. Note that a man-in-the-middle attack is possible here against the Diffie-Hellman part of the solution; this point is further discussed in 4.4.

Case B>A and B starts: this case is not managed in a particular way : since messages are periodically broadcasted one cannot predict if A or B starts the protocol. Here, B starts as in the case without trust. Since A does not trust B, A continues as in the case without trust. Only then B has the possibility to detect the trust relation. At this point the situation is exactly the same as when A started the communication. In particular B will drop its old pseudonym {j.0.db} K -1

B

for choosing a pseudonym of the form{j.k.db} K A very much like in step 3 when A starts.

Case A>B and B>A: this case is not managed in a particular way: depending on which node starts the case A>B or the case B>A is resolved first. Nevertheless, once the anonymous secure channel is built between A and B, it is always possible for a node to ask the authentication of the other node, for instance by signing a random value. For instance, if the case B>A is resolved first, A can authenticate itself to B by sending:

A ! All : {i.0.db} K -1 A .{i.j.db} K A .g l .{n.{n} K -1 A } g ij , B checks {n} K -1 A = n 0 and n 0 = n
4 Analysis of the solution

Basic security properties

Our solution must provide confidentiality of communications between nodes A and B when B>A.

We use the AVISPA framework to prove this property. We first write a full specification of the protocol (given in appendix), then we specify the secrecy property and we run the AVISPA detection tools. The specification corresponds to the case B>A and we check the two situations: when A starts the communication and when B starts the communication.

No attack was found. The case A>B is verified by symmetry of the roles of A and B.

Anonymity properties

We first remark that all the participants X to a communication are using pseudonym {i.0.db} K -1

X

. In order to keep the long term secret K X undisclosed, the cryptographic algorithm {} K -1 X must reveal nothing about its key. This is one basic property of asymmetric encryption algorithms.

We also remark that the nodes are regularly updating their pseudonyms. In particular, they choose a new pseudonym each time they want to establish a trusted connection with other nodes. Thus, one single node can use many pseudonyms at a same time: some for untrusted relations, some other for trusted relations. Depending of whether there is trust or not, the way pseudonyms are built varies but neither form of the pseudonyms reveals anything about long term secrets.

Since no permanent secret is revealed, the effective anonymity only relies on the size of the anonymity set. If the anonimity set is restricted to one single device, of course the anonimity does not hold. But in the case of VANET for instance, the typical anonymity set is the set of all pseudonyms used by all vehicles communicating in the attacker's vicinity during the observation period.

Unobservability against passive attackers

Regarding unobservability, our solution exhibits three properties:

1. All communications are broadcasted. 2. All nodes are regularly sending messages. 3. All message components are encrypted or have the form g i .

According to [START_REF] Pfitzmann | Anonymity, unobservability, and pseudonymity: A consolidated proposal for terminology[END_REF] §8, properties 1 and 3 imply receiver unobservability. Properties 2 and 3 imply sender unobservability. These to properties imply message unobservability, which in turns implies the unobservability of the trust relation by the argument hereafter: Ad absurdum, we assume message unobservability but not trust relation unobservability. We consider the shortest message sequence that leads the attacker to observe the trust relation. We consider the last message of this sequence: it is responsible for the detection of the trust relation (otherwise this message can be removed from the sequence, which contradicts the assumption that this is the shortest). We remark that attacker is not able to detect the trust relation without observing this last message (otherwise this also contradicts the assumption that the sequence is the shortest). Thus, reading just the last message, the attacker detects the trust relation. In other terms, the attacker is able to observe the last message, which contradicts message unobservability.

In order to better illustrate the unobservability property, we provide the two tables hereafter. The first table shows the case B>A and what a passive attacker can deduce. The second table shows the case without trust and what a passive attacker can deduce. α, β, γ are observed pseudonyms that the attacker cannot decrypt. The g x are observed g i values that the attacker cannot reduce (she does not know i). The m x are observed encrypted payloads.

Note that some tools that we used are proven complete [START_REF] Turuani | Security of Cryptographic Protocols: Decidability and Complexity[END_REF]: when no attack is found this means that there is no attack involving an arbitrary number of intruder operations. This attack is not serious per se since it only discloses information that nodes are willing to send even without trust relation. But the attack has consequences over the unobservability of the trust relation between the nodes with pseudonyms {i.0.db} K -1 A and {j.0.db} K -1 B : if the active attacker is able to decrypt the Public payload, this means that there is no trust relation between these nodes. Otherwise the attacker cannot decrypt the payload and then deduces that the node with pseudonym {i.0.db} K -1 A detected g i 6 = g. This case only happens if a trust relation exists.

This attack only defeats the unobservability of the trust relation. The anonymity is still guaranteed, since there is still no way for the attacker to infer A from {i.0.db} K -1 A nor to infer identity information from the Public payload.

Countermeasure

The countermeasure repairs unobservability when a message of the form {j.0.db} K -1

B

.{i.0.db} K -1 A .g k .m is received and when the abnormal situation j 6 = k is detected (in the above description of the attack, we had k = 1). In this case, the detecting node can suspect an attack. We modify its behavior so that it falls back to a non-trust behavior instead of continuing to enforce trust. In particular it will compute the Diffie-Hellman key g jk and send public data, like {Public} g jk . Of course, the attacker will be able to decrypt this data. But, since decryption will always work, she will not learn anything about the trust relation.

The trust relation will be successfully detected and enforced only if there is not attack, that is in the nominal case: {j.0.db} K -1

B

.{i.0.db} K -1

A

.g k .m and j = k. It might be argued that the protocol now silently fails when a node detects an attack, thus giving the attacker an additional way to perform a denial-of-service attack. This is true but in our attacker model the attacker already have the possibility to block all messages (as often in presence of mobile wireless communications).

Conclusion

Privacy issues in network communications become more and more important. In particular in MANET or VANET users take care about their privacy and do not want to reveal anything about their activities and personal or professional travels. In certain cases, just detecting that Bob has meet Alice on some place may reveal partial information about industrial secrets or vendor strategies. In this kind of networks, we have provided a solution for detecting a particular trust relation between two nodes in an anonymous and unobservable way. We believe that these two properties will be of first importance in a near future in the design of security protocols, for instance in the RFID research field [START_REF] Van Deursen | Untraceability of rfid protocols[END_REF][START_REF] Pietro | Information confinement, privacy and security in rfid systems[END_REF]. In future work, we plan to address some complexity issues of our solutions: decreasing the complexity of trust detection algorithm, reduce the use of asymmetric cryptography and add appropriate cryptographic puzzles for mitigating the exhaustion of computation resources. It is also suitable to study different trust models, not necessarily based on asymmetric cryptography.

Figure 1

 1 Figure1gives all the notations used in this paper.

If B trusts a lot of nodes the set TB is large and the detection can take a lot of time. There exist many ways to improve the efficiency, like trying to decrypt with the most often used key first.

Fig. 2. Verification of an active attack against unobservability, as found by AVISPA and displayed by SPAN.

Acknowledgment

The author would like to thank Ciarán Bryce for valuable discussions and appreciation. The author would like to thank the WISTP 2009 reviewers for many in depth comments and useful remarks.

Communication when there is no trust Observations of a passive attacker

We see that the passive attacker observes the same information in both tables. In fact, the general expression of what the passive attacker observes at step n in both cases is obs(n):

Active attack against unobservability

We show here an attack found by the OFMC tool of AVISPA. The countermeasure is given in 4.5 but we find it profitable to precisely explain the attack as it is illustrative for:

the concrete benefit of using automated verification tools; -one way to defeat unobservability while keeping anonymity; -the power of active attackers and the practical conditions of the attack.

If there is no trust between A and B (see 4.4) the secrecy of their subsequent communications only depends on the Diffie-Hellman key agreement performed within the third part of each message. Thus the secret only holds against a passive attacker but not against an active attacker because she can perform a man-in-the-middle attack. The practical conditions of this attack may be very complicated but one variant is much simpler. This variant is automatically found by OFMC tool of AVISPA: the figure 4.4 shows the corresponding message sequence chart as redrawn by the tool SPAN [START_REF] Boichut | Using Animation to Improve Formal Specifications of Security Protocols[END_REF]. The attack works as follows:

1. The attacker D captures a pseudonym, for instance {j.0.db} K -1 B .

The node A normally starts a communication by sending:

A ! All : {i.0.db} K -1

A

.R.g i .m 1 (with R some padding)

3. The active attacker D regularly tries to answer A by sending: D ! All : {j.0.db} K -1

B

.{i.0.db} K -1

A

.g.m. Note that the third part of the message is g (in fact any low power of g will work, like g 2 ,g 3 ...). 4. If A continues the protocol normally, it will send a message with an encrypted payload:

.g k .{Public} g i . Note that the encryption key g i is known by any node having received the first message.

We provide here the formal specification of the protocol, as used for the verification of security properties within the AVISPA framework.


Fabien Caubet

Marc Dambrine

Stability of critical shapes for the drag minimization problem in Stokes flow

Keywords: stability of critical shape, drag minimization, optimal profiles, shape calculus, shape Hessian, Stokes equations AMS Classification: 49Q10, 49K40, 35Q93, (76D55)

We study the stability of some critical (or equilibrium) shapes in the minimization problem of the energy dissipated by a fluid (i.e. the drag minimization problem) governed by the Stokes equations. We first compute the shape derivative up to the second order, then provide a sufficient condition for the shape Hessian of the energy functional to be coercive at a critical shape. Under this condition, the existence of such a local strict minimum is then proved using a precise upper bound for the variations of the second order shape derivative of the functional with respect to the coercivity and differentiability norms. Finally, for smooth domains, a lower bound of the variations of the drag is obtained in terms of the measure of the symmetric difference of domains.

Introduction

In fluid mechanics, the study of the minimization of the drag of a body in a fluid (i.e. the computation of optimal profiles) is a very popular problem. A typical application is the study of the geometry of blunt bodies in flow at low Reynolds numbers (see [START_REF] Schäfer | Evaluation of a CFD benchmark for laminar flows[END_REF]). In this work, the parameter is the shape of the body immersed in the fluid and we are interested in studying the stability of this shape optimization problem.

In an optimization problem, the study of the stability consists to know if a critical shape of a functional is a local strict minimizer. This point is particularly crucial in order to make numerical simulations: if the problem is unstable, regularization is required in the numerical minimization of the functional (see for example [START_REF] Afraites | On second order shape optimization methods for electrical impedance tomography[END_REF][START_REF] Badra | Detecting an obstacle immersed in a fluid by shape optimization methods[END_REF]). The question of stability is addressed with the second order derivative of the functional at a critical point. The notion of derivative used in this work is the now classic shape calculus (Hadamard's approach) presented for example in [START_REF] Henrot | Variation et optimisation de formes[END_REF][START_REF] Simon | Differentiation with respect to the domain in boundary value problems[END_REF][START_REF] Simon | Second variations for domain optimization problems[END_REF][START_REF] Sokołowski | Introduction to shape optimization[END_REF].

The strategy we follow is the following: a first step is to prove the existence of second order shape derivatives. It is classically obtained through an implicit function theorem: this is due to Simon [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF] for Stokes equations and to Bello et al. [START_REF] Bello | Optimal shape design for Navier-Stokes flow[END_REF][START_REF] Bello | The variation of the drag with respect to the domain in Navier-Stokes flow[END_REF] for the Navier-Stokes equations). In a second step, we obtain the Euler-Lagrange equation, then we compute the shape Hessian. These shape derivatives of the drag (at least at the first order) was computed by some authors (see for example [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF] for the Stokes equations and [START_REF] Bello | The variation of the drag with respect to the domain in Navier-Stokes flow[END_REF][START_REF] Bello | The differentiability of the drag with respect to the variations of a Lipschitz domain in a Navier-Stokes flow[END_REF] for the Navier-Stokes equations). However, in this work, we need a different expression of the shape Hessian (in [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF]). Then, we use the Euler-Lagrange equations and some computations to derive a sufficient condition of positivity of the shape Hessian computed at a critical shape.

In a third step, we prove the stability of the minimization. In finite dimension, the knowledge of the sign of the Hessian permits to fully answer the question of stability using Taylor-Young expansion. However, in shape optimization or in infinite dimension, the two norms discrepancy problem can occur: the coercivity norm is often weaker than the differentiability norm. Let us refer to the paper [START_REF] Descloux | On the two dimensional magnetic shaping problem without surface tension[END_REF] by Descloux for a concrete example of such a situation known as the magnetic shaping problem: on this example, the coercivity at a critical point holds in the H 1/2 norm while the differentiability holds in C 2 topology. Since the quantity o( • C 2 ) is not smaller than C • 2 H 1/2 , the classical argument using the Taylor-Young formula does not insure that this critical point is a local strict minimum. A method to overcome this problem is given by Dambrine et al. in [START_REF] Dambrine | About stability of equilibrium shapes[END_REF][START_REF] Dambrine | On variations of the shape Hessian and sufficient conditions for the stability of critical shapes[END_REF]. The key is a precise estimate of the variations of the second order shape derivative of the functional with respect to the coercivity and differentiability norms around a critical shape. In [START_REF] Dambrine | About stability of equilibrium shapes[END_REF][START_REF] Dambrine | On variations of the shape Hessian and sufficient conditions for the stability of critical shapes[END_REF], the authors study Poisson's equation or the case of a strictly and uniformly elliptic operator. Here, we adapt these methods to the drag minimization problem for a Stokes flow. Finally, in a fourth and last step, we derive a precise version of the minimality inequality: we provide a lower bound of the variations of the drag in term of a geometrical quantity following new ideas introduced by Fusco et al. [START_REF] Fusco | Equilibrium configurations of epitaxially strained elastic films: second order minimality conditions and qualitative properties of solutions[END_REF][START_REF] Acerbi | Minimality via second variation for a nonlocal isoperimetric problem[END_REF].

The paper is organized as follows. We first define the studied problem: we define the shape functional under consideration and the used domain perturbations. In particular, we recall some geometrical constructions from [START_REF] Dambrine | About stability of equilibrium shapes[END_REF][START_REF] Dambrine | On variations of the shape Hessian and sufficient conditions for the stability of critical shapes[END_REF]. In Section 3, following Hadamard's approach, we use shape calculus to prove the coercivity of the shape Hessian of the considered functional at an equilibrium shape in the H 1/2 norm. Then, in Section 4, we state the main result of this work which is a stability result for the drag minimization problem. It claims that a critical shape can be a local strict minimum if a given criterion is satisfied. In Section 5, we translate our result in a more explicit estimation in purely geometrical quantities. Finally, we conclude the paper in Section 6 mentioning some possible extensions and prospects concerning the problem we study here.

The problem setting

General notations

Let us introduce the notations that we adopt in this paper. For an open set Ω ⊂ R d (d = 2 or 3), we denote by L p (Ω), W m,p (Ω) and H s (Ω) the usual Lebesgue and Sobolev spaces. We note in bold the vectorial functions and spaces: L p (Ω), W m,p (Ω), H s (Ω), etc. Moreover, for k ∈ N and α ∈ (0, 1), the space C k (Ω) is defined as the set of functions having continuous derivatives up to order k in Ω and C k,α (Ω) denotes the usual Hölder space. We denote by |Ω| the measure of Ω. Moreover, n represents the external unit normal to ∂Ω, and for a smooth enough function u, we note ∂ n u the normal derivative of u. Finally, we define the space

L 2 0 (Ω) := p ∈ L 2 (Ω), Ω p = 0 .

Tubular coordinates and normal representation of a perturbed domain

Let us recall some basic facts from differential geometry and fix the notations. We follow in this paragraph the constructions and proof of [START_REF] Dambrine | About stability of equilibrium shapes[END_REF]. Let Ω be a smooth bounded open set and let n denote the outer unit normal field to ∂Ω. There is a non negative real ρ(Ω), such that the application T ∂Ω defined by

T ∂Ω : ∂Ω × R → R d (M, h) → M + h n(M ) is a local diffeomorphism from ∂Ω × (-ρ(Ω), ρ(Ω)
) on an open tubular neighborhood of ∂Ω we will denote T ∂Ω . This property expresses the fact that any point x in T ∂Ω has a unique orthogonal projection p ∂Ω (x) on ∂Ω and that the relation

x = p ∂Ω (x) + h(x) n(p ∂Ω (x)) where |h(x)| = x -p ∂Ω (x) ,
holds. Notice that h is uniquely defined. Let us recall the definition of the Micheletti's distance on the set S of C 2,α bounded domains of R d that are diffeomorphic to Ω. The idea is to use the Banach imbedding of diffeomorphisms in

C 2,α (R d , R d ). For Ω 1 and Ω 2 in S, set d 2,α (Ω 1 , Ω 2 ) := inf Θ -I 2,α + Θ -1 -I 2,α
where the infimum is taken over all the diffeomorphisms

Θ ∈ C 2,α (R d ) mapping Ω 1 onto Ω 2 .
In particular, if Θ is a perturbation of the identity such that

Θ -I 2,α + Θ -1 -I 2,α < ρ(Ω),
any point of the boundary ∂Ω per of the perturbed domain Ω per := Θ(Ω) lays in the tubular neighborhood T ∂Ω of ∂Ω so that it can be described in terms of normal deformations: for all

x ∈ ∂Ω Θ(x) = R Θ (x) + d ∂Ω (Θ(x)) n(R Θ (x)),
where R Θ is a diffeomorphism from ∂Ω into itself and d ∂Ω is the signed distance to ∂Ω. Note that n := ∇d ∂Ω is a unitary extension of the normal field that coincides with n • p ∂Ω in the tubular neighborhood T ∂Ω . As a consequence, the flow Φ t of the vector field X Θ defined in T ∂Ω by

X Θ (x) = h(x) n(p ∂Ω (x)) (2.1)
defines a path t ∈ [0, 1] → Ω t within domains connecting Ω 0 = Ω to Ω 1 = Ω per . However, this flow does not in general preserve the volume of the domain even if the original and perturbed domains share the same volume. Therefore, an alternative field Y Θ defined in tubular coordinates (M, h) by

Y Θ (M, h) = - d ∂Ω (Θ•R -1 Θ (M )) 0 det DT ∂Ω (M, s) ds det DT ∂Ω (M, h) n(M ) (2.2)
was build in [9, Subsection 2.1] in order to deal with volume constraint. Indeed, the divergence of Y Θ cancels in the tubular neighborhood T ∂Ω and the flow of Y Θ preserves the measure of Ω for t in [0, 1].

Let us remark that the field X Θ involves the normal field n to ∂Ω and that Y Θ involves first order derivatives of n. Hence their flow is a C 2,α -diffeomorphism that is to say let stable the admissible class of domain under additional regularity assumptions. Let us assume from now on that the initial shape is in fact C 4,α . For time t ∈ [0, 1], the flow Φ t of both fields X Θ and Y Θ remains a perturbation of the identity: there is a constant C such that, for all t ∈ [0, 1]:

Φ t -I 2,α ≤ C Θ -I 2,α . (2.3)
Noting that the outer normal field n t to Ω t can be written as

n t (y) = ∇(d ∂Ω • Φ -1 t ) ∇(d ∂Ω • Φ -1 t ) (y),
J(t, .) -1 C 1 (∂Ω) ≤ C Φ t -I 2,α , ∀t ∈ [0, 1]; (2.4)
• the normal field n t to ∂Ω t satisfies

n t (Φ t (.)) -n(.) C 1 (∂Ω) ≤ C Φ t -I 2,α , ∀t ∈ [0, 1].
(2.5)

• Set m V := V • n for V = X Θ or Y Θ , then for all t ∈ [0, 1]: m V (Φ t (.)) -m V L 2 (∂Ω) ≤ C m V L 2 (∂Ω) Φ t -I 2,α , (2.6)
m V (Φ t (.)) -m V H 1/2 (∂Ω) ≤ C m V H 1/2 (∂Ω) Φ t -I 2,α . (2.7)
Remark 2.2. The fact that we have to impose a C 4,α regularity of the boundary of the initial shape (whereas we can work with shapes with a C 2,1 boundary in order to have the twice differentiability with respect to the domain) comes to the fact that we consider normal perturbations (so we loose one rank of regularity) which are divergence free (which imposes the lost of an additional derivative by construction): see [9, Section 2.1].

The general notations and the functional

The drag functional. Let Ω be a bounded and connected open subset of R d (with d = 2 or d = 3) containing a Newtonian and incompressible fluid with coefficient of kinematic viscosity ν > 0. We assume that Ω is smooth (at least with a C 2,1 boundary). We also assume that ν is constant. Let δ > 0 fixed (large). We define the set of admissible shapes by

O δ := ω ⊂ Ω open set with a C 2,1 boundary such that d(x, ∂Ω) > δ, ∀x ∈ ω
and such that Ω\ω is connected .

Let us consider g ∈ C 2,α (∂Ω), with α ∈ (0, 1), satisfying the compatibility condition ∂Ω g • n = 0 and the unique solution

(u, p) ∈ C 2,α (Ω\ω) × C 1,α (Ω\ω) ∩ L 2 0 (Ω\ω) of        -ν∆u + ∇p = 0 in Ω\ω div u = 0 in Ω\ω u = g on ∂Ω u = 0 on ∂ω.
The existence and uniqueness of the solution of such a problem is classical. We refer for example to the book of Galdi [14, Theorems IV.7.1 and IV.7.2]. The energy dissipated by the fluid is given by

J(ω) := 1 2 Ω\ω ν |D(u)| 2 ,
where

D(u) := ∇u + t ∇u = (∂ i u j + ∂ j u i ) i,j , i, j = 1, . . . , d.
The drag minimization problem is to minimize J over all subdomains ω of Ω with a given measure M :

ω * := argmin J(ω); ω ∈ O δ such that L d (ω) = M ,
where L d denotes the Lebesgue measure on R d . We also define the following stress tensor:

σ(u, p) := ν ∇u + t ∇u -p I.
Admissible deformations. Let ω ∈ O δ . Let us define some admissible deformations of the domain ω we will use in this paper. Since we want to perturb only ω (and not Ω which is fixed), we define Ω δ an open set with a C ∞ boundary and such that

{x ∈ Ω ; d(x, ∂Ω) > δ/2} ⊂ Ω δ ⊂ {x ∈ Ω ; d(x, ∂Ω) > δ/3} .
We then consider a diffeomorphism

Θ ∈ U := θ ∈ C 2,1 (R d ), θ ≡ I in R d \Ω δ .
The principle of Hadamard's approach is to consider the flow Φ Θ,t of an adequate autonomous vector field

V Θ ∈ U := θ ∈ C 2,1 (R d ); Supp θ ⊂ Ω δ , i.e. the solution of ∂ t Φ = V Θ (Φ) Φ(0, x) = x. Notice that V Θ represents one of the previous fields X Θ or Y Θ defined in Section 2.2. It defines a regular path (ω t := Φ Θ,t (ω)) t∈[0,1] in O δ .
For the rest of the paper, we use a subscript " t " to indicate that the quantity is defined on the time t dependent domain. For instance, n t is the external unit normal of Ω\ω t . Moreover, in order to simplify the notations, we omit the dependency with respect to Θ: in particular, we use the notation Φ t instead of Φ Θ,t .

Let us now construct some vector fields whose flows connects the original shape ω * to the perturbed domain ω 1 = Θ(ω).

Computation of the shape Hessian

In this section, we prove that the shape Hessian of the functional J can be coercive for a critical shape ω * ∈ O δ with a C 4,α boundary in the H 1/2 (∂ω * ) norm. We assume the existence of such a critical shape. We first define, for t ∈ [0, 1],

j(t) := J(ω t ) := 1 2 Ω\ωt ν |D(u t )| 2 ,
where

(u t , p t ) ∈ C 2,α (Ω\ω t ) × C 1,α (Ω\ω t ) ∩ L 2 0 (Ω\ω t )
is the solution of the following perturbed problem:

       -ν∆u t + ∇p t = 0 in Ω\ω t div u t = 0 in Ω\ω t u t = g on ∂Ω u t = 0 on ∂ω t , (3.1)
with ω t := Φ t (ω * ). The existence of the second order shape derivatives is proved in some previous works (see for example [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF][START_REF] Badra | Detecting an obstacle immersed in a fluid by shape optimization methods[END_REF][START_REF] Caubet | A Kohn-Vogelius formulation to detect an obstacle immersed in a fluid[END_REF] for the Stokes case and [START_REF] Bello | Optimal shape design for Navier-Stokes flow[END_REF][START_REF] Caubet | Instability of an inverse problem for the stationary Navier-Stokes equations[END_REF] for the Navier-Stokes case). Let V ∈ U be divergence free. Classically, the shape derivative

(u ′ t , p ′ t ) which belongs to C 2,α (Ω\ω t ) × C 1,α (Ω\ω t ) ∩ L 2 0 (Ω\ω t ) is characterized as the solution of the following Stokes problem:        -ν∆u ′ t + ∇p ′ t = 0 in Ω\ω t div u ′ t = 0 in Ω\ω t u ′ t = 0 on ∂Ω u ′ t = -∂ nt u t (V • n t ) on ∂ω t . (3.2)

Shape gradient of the functional and Euler-Lagrange equation

Simon proves in [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF]Theorem 3] that

DJ(ω * ) • V = - 1 2 ∂ω * ν |∂ n u| 2 (V • n).
Since we work under the constraint of constant volume and since ω * is a critical point, there exists

Λ 0 ∈ R such that DJ(ω * ) • V + Λ 0 DV(ω * ) • V = 0 for any perturbation V (where V(ω * ) is the volume of ω * ). Hence, for all V ∈ U , - 1 2 ∂ω * ν |∂ n u| 2 (V • n) + Λ 0 ∂ω * V • n = 0,
and then we obtain the Euler-Lagrange equation satisfied at the critical shape: there exists

Λ ∈ R such that |∂ n u| 2 = 2 ν Λ 0 =: Λ. (3.3)
Notice that this result is also proved in [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF]Theorem 7] in a different way.

Shape Hessian of the functional

Let T be a fixed small positive real number. Simon proves in [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF]Theorem 4] that, for t ∈ [0, T ),

j ′ (t) = - 1 2 ∂ωt ν |D(u t )| 2 (V • n).
In [START_REF] Simon | Domain variation for drag in Stokes flow[END_REF]Theorem 8], the second variation j ′′ (0) is computed. However, we shall use a simplified expression of j ′′ (0) which, in particular, does not depend on the second order shape derivative of the state u. Let us obtain this new expression. Since V = 0 on ∂Ω, we use Green's formula to get

j ′ (t) = - 1 2 Ω\ωt νdiv |D(u t )| 2 V .
Thus, we use Hadamard's formula (see [START_REF] Henrot | Variation et optimisation de formes[END_REF]Corollary 5.2.5]) to derive this function. For t ∈ [0, T ), it holds

j ′′ (t) = - Ω\ωt ν div (D(u t ) : D(u ′ t ))V - 1 2 ∂ωt ν div |D(u t )| 2 V (V • n t ) = - 1 2 ∂ωt ν 2(D(u t ) : D(u ′ t )) + ∇(|D(u t )| 2 ) • V + |D(u t )| 2 div V (V • n t ) = - 1 2 ∂ωt ν 2(D(u t ) : D(u ′ t )) + ∇(|D(u t )| 2 ) • V (V • n t ), (3.4)
since we work with a divergence free vector field V . Moreover,

∂ωt ν D(u t ) : D(u ′ t ) (V • n t ) = 2 ∂ωt ν ∇u t : D(u ′ t )(V • n t ) = 2 ∂ωt ∇u t : (ν D(u ′ t ) -p ′ t I)(V • n t ).
Here we used that the solution u t of Stokes problem is divergence free in Ω\ω t to get

∂ωt ∇u t : (p ′ t I)(V • n t ) = ∂ωt p ′ t div u t (V • n t ) = 0.
Moreover, u t satisfies the Dirichlet boundary condition u t = 0 on ∂ω t , and then it follows that ∇u t = ∂ nt u t ⊗ n t on ∂ω t . Thus,

∂ωt νD(u t ) : D(u ′ t ) (V • n t ) = 2 ∂ωt ∇u t : σ(u ′ t , p ′ t )(V • n t ) = 2 ∂ωt ∂ nt u t • σ(u ′ t , p ′ t ) n t (V • n t ). (3.5)
Then, combining Equalities (3.4) and (3.5),

j ′′ (t) = -2 ∂ωt ∂ nt u t • (σ(u ′ t , p ′ t )n t )(V • n t ) - 1 2 ∂ωt ν∇(|D(u t )| 2 ) • V (V • n t ).
Finally, since the shape derivative u ′ t satisfies (3.2) and in particular

u ′ t = -∂ nt u t (V • n t ) on ∂ω t , j ′′ (t) = 2 ∂ωt u ′ t • (σ(u ′ t , p ′ t )n t ) - 1 2 ∂ωt ν∇(|D(u t )| 2 ) • V (V • n t ). (3.6)

At a critical shape

We recall that ω * is a critical shape. Now, we study the expression of j ′′ (0) and then use the fact that ω * is a critical shape. The previous characterization of j ′′ (t) gives, for t = 0, the expression

j ′′ (0) = 2 ∂ω * u ′ • (σ(u ′ , p ′ )n) =:J 1 - 1 2 ∂ω * ν∇(|D(u)| 2 ) • V (V • n) =:J 2 . (3.7)
that we split into two terms we will study separately.

D N : H 1/2 (∂ω) -→ H -1/2 (∂ω) h -→ σ(U , P )n,
where

(U , P ) ∈ H 1 (Ω\ω) × L 2 (Ω\ω) solves        -div (σ(U , P )) = 0 in Ω\ω div U = 0 in Ω\ω U = 0 on ∂Ω U = h on ∂ω. (3.8)
Then, there exists a constant c depending of Ω and ω such that

D N (h) , h H -1/2 (∂ω),H 1/2 (∂ω) ≥ c h 2 H 1/2 (∂ω) . (3.9)
Proof. Using U as a test function in Problem (3.8), we obtain

D N (h) , h H -1/2 (∂ω),H 1/2 (∂ω) = 1 2 ν D(U ) 2 L 2 (Ω\ω) .
Then, since U = 0 on ∂Ω, Korn's inequality (see for example [16, eq. (2.14) page 19]) leads

D N (h) , h H -1/2 (∂ω),H 1/2 (∂ω) ≥ c U 2 H 1 (Ω\ω) ≥ c h 2 H 1/2 (∂ω) .
Consider Problem (3.2) (for t = 0) solved by the shape derivative (u ′ , p ′ ). Then, noticing that, since div u ′ = 0,

-ν∆u ′ + ∇p ′ = -div (νD(u ′ )) + ∇p ′ = -div (σ(u ′ , p ′ )),
and that u ′ = -∂ n u (V • n) on ∂ω * , we get ∂ω * u ′ • (σ(u ′ , p ′ )n) = ∂ω * ∂ n u(V • n) D N (∂ n u(V • n)).
Hence, there exists a constant c > 0 such that

J 1 ≥ c ∂ n u (V • n) 2 H 1/2 (∂ω * ) = c Λ V • n 2 H 1/2 (∂ω * ) . (3.10)
We used the coercivity of the operator D N ((3.9) in Lemma 3.1) and the Euler-Lagrange equation |∂ n u| 2 = Λ according to (3.3)

Second term J 2 . Let us now study the second term of the characterization of j ′′ (0) given by (3.7). We recall that, according to subsection 2.2, only normal perturbation directions are considered. Thus, we focus on

∂ω * ν∇(|D(u)| 2 ) • n (V • n) 2 .
We have to compute ν∇(|D(u)| 2 ) • n. This is the object of next Lemma. • F is the second fundamental form of the surface ∂ω * ;

• ∇ τ denotes the tangential gradient.

Proof. By using local coordinates as in [START_REF] Descloux | On the two dimensional magnetic shaping problem without surface tension[END_REF] or [START_REF] Dambrine | About stability of equilibrium shapes[END_REF][START_REF] Dambrine | On variations of the shape Hessian and sufficient conditions for the stability of critical shapes[END_REF], we can assume that ω * is locally above the graph of a function f :

(-ε, ε) d-1 → R (of class C 4,α ) which is such that f (0) = ∂ i f (0) = 0 for i = 1, . . . , d -1. For all x ∈ (-ε, ε) d-1 , noticing n = (n 1 , . . . , n d ), ∇ |D(u)(x, f (x))| 2 • n(x, f (x)) = ∇   d i,j=1 (∂ i u j (x, f (x)) + ∂ j u i (x, f (x))) 2   • n(x, f (x)) = d k=1 ∂ k   d i,j=1 (∂ i u j (x, f (x)) + ∂ j u i (x, f (x))) 2   n k (x, f (x)) = 2 d i,j,k=1 (∂ 2 ki u j (x, f (x)) + ∂ 2 kj u i (x, f (x)))(∂ i u j (x, f (x)) + ∂ j u i (x, f (x)))n k (x, f (x)).
Then, for x = 0, since n(0, f (0)) = (0, . . . , 0, 1) in our system of coordinates,

∇ |D(u)(0, f (0))| 2 • n(0, f (0)) = 2 d i,j=1 (∂ 2 di u j (0, f (0)) + ∂ 2 dj u i (0, f (0)))(∂ i u j (0, f (0)) + ∂ j u i (0, f (0))).
Moreover, since u = 0 on ∂ω * , we have ∂ i u j (0, f (0)) = 0 for all i = 1, . . . , d -1 and j = 1, . . . , d. In addition, since div u = 0 in Ω\ω, it follows that

∂ d u d (0, f (0)) = 0. Then, ν∇ |D(u)(0, f (0))| 2 • n(0, f (0)) = 4 d-1 j=1 ν ∂ 2 dd u j (0, f (0)) + ∂ 2 dj u d (0, f (0)) ∂ d u j (0, f (0)).
(3.12) Let us now interpret the boundary condition on ∂ω * , that is

u(x, f (x)) = 0 for x ∈ (-ε, ε) d-1 .
By differentiation, for i, j = 1, . . . , d -1 and k = 1, . . . , d,

∂ j u k (x, f (x)) + ∂ j f (x)∂ d u k (x, f (x)) = 0 and ∂ 2 ij u k (x, f (x)) + ∂ i f (x)∂ 2 dj u k (x, f (x)) + ∂ 2 ij f (x)∂ d u k (x, f (x)) + ∂ j f (x) ∂ 2 id u k (x, f (x)) + ∂ i f (x)∂ 2 dd u k (x, f (x)) = 0. For x = 0, since f (0) = ∂ i f (0) = 0 for i = 1, . . . , d -1, ∂ 2 ij u k (0, f (0)) + ∂ 2 ij f (0)∂ d u k (0, f (0)) = 0 ∀i, j = 1, . . . , d -1, ∀k = 1, . . . , d. (3.13)
Notice that -ν∆u+∇p = 0 in Ω\ω, i.e., since div u = 0 in Ω\ω, -νdiv (D(u))+∇p = 0 in Ω\ω. Hence, for j = 1, . . . , d and any

x ∈ (-ε, ε) d-1 , ν ∂ 2 dd u j (x, f (x))+∂ 2 dj u d (x, f (x)) = - d-1 i=1 ν ∂ 2 ii u j (x, f (x))+∂ 2 ij u i (x, f (x)) +∂ j p(x, f (x)) . (3.14)
In Equality (3.13) with j = i, we sum on i to get, since

d-1 i=1 ∂ 2 ii f (0) = (d -1)H(0)
, where H is the mean curvature of ∂ω * (see for example the comment after Definition 5.4.7 in [START_REF] Henrot | Variation et optimisation de formes[END_REF]), that for all k = 1, . . . , d

d-1 i=1 ∂ 2 ii u k (0, f (0)) + (d -1)H(0)∂ d u k (0, f (0)) = 0. (3.15)
Hence, using (3.13) and (3.15), Equality (3.14) writes as follows for j = 1, . . . , d -1:

ν ∂ 2 dd u j (0, f (0)) + ∂ 2 dj u d (0, f (0)) = ν (d -1)H(0) ∂ d u j (0, f (0)) + ν d-1 i=1 ∂ 2 ij f (0)∂ d u i (0, f (0)) + ∂ j p(0, f (0)).
Therefore, Equality (3.12) writes in the following form:

ν∇ |D(u)(0, f (0))| 2 • n(0, f (0)) = 4 d-1 j=1 ν(d -1)H(0) (∂ d u j (0, f (0))) 2 + ν d-1 i=1 ∂ 2 ij f (0)∂ d u i (0, f (0))∂ d u j (0, f (0)) + ∂ j p(0, f (0))∂ d u j (0, f (0)) .
Hence, since |∂ d u| 2 = Λ (see (3.3)), Lemma 3.2 follows.

Conclusion.

Finally, gathering (3.7), (3.10) and (3.11), we obtain

D 2 J(ω * ) • V • V ≥ cΛ V • n 2 H 1/2 (∂ω * ) -2 ∂ω * ν(d -1)HΛ + ν t b F b + ∇ τ p • b (V • n) 2 .
Hence, using the regularity of ∂ω * and of the solution u, the quantity

ν(d -1)HΛ + ν t b F b + ∇ τ p • b
belongs to L ∞ (∂ω * ) and the second term behaves like V • n 2 L 2 (∂ω * ) . Hence, a natural assumption is that the shape Hessian is coercive in the H 1/2 (∂ω * ) sense. This is the case if

ν(d -1)HΛ + ν t b F b + ∇ τ p • b < 0 on ∂ω * . (3.16)
Note that this condition couples geometrical effect with the solution itself. In practice, the condition (3.16) cannot easily be tested theoretically since it couples the curvature of the object with derivatives of the flow. It might be tested numerically, however it will requires a curved mesh for the surface and the fluid domain and high order elements to catch the desired effects of curvature and the derivatives of the couple (u, p). Such a precise computation requires specific numerical attention.

Statement of the result

We have proved the following theorem:

Theorem 3.3. Let us assume that Estimate (3.16) is satisfied. Then, there exists a constant c > 0 such that, for all V ∈ U with div V = 0,

j ′′ (0) = D 2 J(ω * ) • V • V ≥ c V • n 2 H 1/2 (∂ω * ) .
4 Stability of the drag minimization problem

The results

We recall that ω * ∈ O δ is assumed to have a C 4,α boundary (α ∈ (0, 1)). Moreover, in this section, the shape Hessian is assumed to be coercive at ω * in the following sense: for any V in the tangent space defined by the constraints that is with

∂ω * V • n = 0, it holds j ′′ (0) = D 2 J(ω * ) • V • V ≥ c V • n 2 H 1/2 (∂ω * ) . (4.1)
The main result of this section states the stability of the drag minimization problem:

Theorem 4.1. If ω * is a critical shape for J where (4.1) holds, there exists η > 0 such that, for all Θ ∈ U with Θ -

I C 2,1 (R d ) < η, |Θ(ω * )| = |ω * | and Θ = I, J(Θ(ω * )) > J(ω * ).
Consider the admissible vector field Y Θ ∈ U defined in subsection 2.2. By assumption

j ′′ (0) ≥ C Y Θ • n 2 H 1/2 (∂ω * ) .
By the order two Taylor expansion

J(ω) = j(1) = j(0) + 1 0 (1 -x)j ′′ (x)dx = J(ω * ) + 1 0 (1 -x)j ′′ (x)dx,
this theorem is a direct consequence of the following theorem (and of the assumption of the H 1/2 coercivity of the shape Hessian):

Theorem 4.2. There exists η 0 > 0 and a function w : (0, η 0 ) → R with lim rց0 w(r) = 0 (which depends only on Ω, ω * and the data) such that, for all η ∈ (0, η 0 ) and all Θ ∈ U with

Θ -I C 2,1 (R d ) < η and |Θ(ω * )| = |ω * |, there exists a divergence free vector field Y Θ ∈ U constructed in (2.
2) whose the flow Φ t defines a path (ω t := Φ t (ω * )) t∈[0,1] between ω * and Θ(ω * ), such that, for all t ∈ [0, 1], the following estimate holds:

j ′′ (t) -j ′′ (0) ≤ w(η) Y Θ • n 2 H 1/2 (∂ω * ) .
Indeed, there is a non negative η such that w(η) ≤ C/2, so that

J(ω) ≥ J(ω * ) + C 4 Y Θ • n 2 H 1/2 (∂ω * ) > J(ω * ). (4.2)
In the following section, we prove Theorem 4.2 that extends to the Stokes case [9, Theorem 3] (or [10, Theorem 2.1]). We follow the same strategy that the one used in these references.

Let Θ ∈ U with Θ -I C 2,1 (R d ) small enough. For t ∈ [0, 1], we consider the solution

(u t , p t ) ∈ C 2,α (Ω\ω t ) × C 1,α (Ω\ω t ) ∩ L 2
0 (Ω\ω t ) of the perturbed Stokes problem (3.1). We define v t := u t • Φ t and q t := p t • Φ t .

Hence, using the change of variables x = Φ t y in Problem (3.1), we prove that

       -div (∇v t A(t)) + div (q t B(t)) = 0 in Ω\ω * ∇v t : B(t) = 0 in Ω\ω * v t = g in ∂Ω v t = 0 in ∂ω * , (4.3)
with

Jac t := det (DΦ t ) ∈ C 1,1 (R d ) A(t) := Jac t ν (DΦ t ) -1 t DΦ t -1 ∈ C 1,1 (R d , M d,d ) B(t) := Jac t t DΦ t -1 ∈ C 1,1 (R d , M d,d ).

Preliminary results: bounds on the state function

Let us first focus on the proof of the following preliminary result:

u t • Φ t -u C 2 (Ω\ω * ) ≤ w Θ -I C 2,1 (R d ) sup t∈[0,1] p t • Φ t -p C 1 (Ω\ω * ) ≤ w Θ -I C 2,1 (R d ) .
Notice that we can assume the existence of a constant C such that:

w(η) ≥ Cη. (4.4)
since if w(η) < Cη, we can redefine w as w(η) = Cη (which tends to 0 with η).

First, let us note that, using a differentiability result similar than the one given in [3, Lemma 3.2], there exist η 2 > 0 and a constant c > 0, depending only on the data, such that, for all Θ ∈ U with Θ -

I C 2,1 (R d ) < η 2 , v t -u C 2,α (Ω\ω * ) ≤ c Θ -I C 2,1 (R d ) and q t -p C 1,α (Ω\ω * ) ≤ c Θ -I C 2,1 (R d ) .
Hence, we deduce immediately the following lemma: Lemma 4.4. There exists η 2 > 0 such that, for all Θ ∈ U with Θ -I C 2,1 (R d ) < η 2 , there exists a constant C depending only on the data and on

Θ C 2,1 (R d ) such that, for all t ∈ [0, 1], v t C 2,α (Ω\ω * ) ≤ C and q t C 1,α (Ω\ω * ) ≤ C.
w(η) := sup Θ∈U, Θ-I C 2,1 (R d ) <η t∈[0,1] v t -u C 2 (Ω\ω * ) + q t -p C 1 (Ω\ω * ) .
Lemma 4.4 guarantees that this quantity is well defined. Since the estimate of Proposition 4.3 is obvious, it remains to prove that lim rց0 w(r) = 0.

Let us proceed by contradiction assuming that lim rց0 w(r) = 0. Then, there exists a real a > 0 and some sequences

(t n ) n∈N ⊂ [0, 1], (Θ n ) n∈N ⊂ U with Θ n -I C 2,1 (R d ) < η 2 for all n ∈ N, (ω n := Φ Θn,tn (ω * )) n∈N , (v n := u n • Φ Θn,tn ) n∈N and (q n := p n • Φ Θn,tn ) n∈N such that for all n ∈ N Θ n -I C 2,1 (R d ) ≤ 1 n v n -u C 2 (Ω\ω * ) + q n -p C 1 (Ω\ω * ) ≥ a > 0. (4.5) Using Lemma 4.4, the sequence (v n ) n is bounded in C 2,α (Ω\ω * ) and (q n ) n is bounded in C 1,α (Ω\ω * ). Hence, since C 2,α (Ω\ω * ) (respectively C 1,α (Ω\ω * )) is compactly imbedding in C 2 (Ω\ω * ) (respectively in C 1 (Ω\ω * ))
, there exist two sub-sequences converging respectively in C 2 (Ω\ω * ) to some u lim and in C 1 (Ω\ω * ) to some p lim . Moreover, there exists a sub-sequence (t n ) ⊂ [0, 1] which converges to some t lim . Then, notice that

-div (∇v n A(t n , Θ n )) + div (q n B(t n , Θ n )) + ν∆u lim -∇p lim = -div (∇v n A(t n , Θ n )) + div (q n B(t n , Θ n )) + ν∆v n -∇q n -ν∆ (v n -u lim ) + ∇(q n -p lim ), in Ω\ω * and ∇v n : B(t n , Θ n ) -div u lim = ∇v n : B(t n , Θ n ) -div v n + div (v n -u lim ) in Ω\ω * .
Let us focus on the second equality. We proceed in the same way for the first one. Since v n → u lim in C 2 (Ω\ω * ), we obtain, passing to the limit in the previous equality and using Problem (4.3) solved by the couple (v n , q n ), that

div u lim = 0 in Ω\ω * .
Thus, using the boundary conditions, we obtain that the couple (u lim , p lim ) solves

       -ν∆u lim + ∇p lim = 0 in Ω\ω * div u lim = 0 in Ω\ω * u lim = g on ∂Ω u lim = 0 on ∂ω * .
The uniqueness of such a solution implies (u lim , p lim ) = (u, p) which contradicts the second inequality of (4.5).

Proof of Theorem 4.2

Again, we fix Θ ∈ U a diffeomorphism with Θ -

I C 2,1 (R d ) < η 2 such that |Θ(ω * )| = |ω * |.
As previously, we consider the vector field Y Θ , the associated flow Φ t and the domains (ω t := Φ t (ω * )) t∈[0,1] and we set m Y Θ := Y Θ • n. The constant η 2 is given by Lemma 4.4. We denote by C any constant depending only on

Θ C 2,1 (R d ) .
We have proved (see (3.6)) that

D 2 J(ω t ) • Y Θ • Y Θ = - 1 2 ∂ωt ν ∇(|D(u t )| 2 ) • Y Θ (Y Θ • n t ) =:I 1 (t) + 2 ∂ωt u ′ t • (σ(u ′ t , p ′ t )n t ) =:I 2 (t)
.

This decomposition is similar as the one in J 1 and J 2 for the Hessian at the critical shape.

It is natural since the two depends on the deformation parameter Y Θ through two distinct norms. Hence, we will prove the following two lemmas:

Lemma 4.5. We have

|I 1 (t) -I 1 (0)| ≤ C w( Θ -I C 2,1 (R d ) ) m Y Θ 2 L 2 (∂ω * ) .
Lemma 4.6. We have

|I 2 (t) -I 2 (0)| ≤ C w( Θ -I C 2,1 (R d ) ) m Y Θ 2 H 1/2 (∂ω * ) .
In both lemmas, w denotes the modulus of continuity introduced in Proposition 4.3. Using the previous characterization of D 2 J(ω t )•Y Θ •Y Θ and these two lemmas, Theorem 4.2 is proved. Hence, it now suffices to prove these lemmas. To simplify the notations, we set

η := Θ -I C 2,1 (R d ) ,
and, in the sequel, we denote all the constants depending only on the data and on the norm Θ C 2,1 (R d ) by C.

Proof of Lemma 4.5

We know that

∇u t • Φ t = t [DΦ t ] -1 ∇v t and D(u t ) • Φ t = t [DΦ t ] -1 D(v t ). Hence, replacing u t by |D(u t )| 2 , ∇ |D(u t )| 2 • Φ t = t [DΦ t ] -1 ∇(|D(u t )| 2 • Φ t ) = t [DΦ t ] -1 ∇ t [DΦ t ] -1 D(v t ) 2 .
Thus, using the change of variables x = Φ t y,

I 1 (t) = ∂ω * ν t DΦ -1 t ∇ t DΦ -1 t D(v t ) 2 • (Y Θ • Φ t ) [(Y Θ • Φ t ) • (n t • Φ t )] J(t) = ∂ω * ν |m Y Θ • Φ t | 2 t DΦ -1 t ∇ t DΦ -1 t D(v t ) 2 • n (n • (n t • Φ t )) J(t), (4.6)
where J(t) := det DΦ t / ( t DΦ -1 t )n is the surface jacobian. In order to study I 1 (t) -I 1 (0), let us define

a 0 (t) := m Y • Φ t , a 1 (t) := t DΦ -1 t ∇ t DΦ -1 t D(v t ) 2 • n, a 2 (t) := n • (n t • Φ t ) , a 3 (t) := J(t).
Using the estimate concerning n t • Φ t given in Lemma 2.1, we have

a 2 (t) L ∞ (∂ω * ) = n • (n t • Φ t ) L ∞ (∂ω * ) ≤ C a 2 (t) -a 2 (0) L ∞ (∂ω * ) = n • (n t • Φ t -n) L ∞ (∂ω * ) ≤ C η. (4.7)
Moreover, using the estimate concerning J(t) given in Lemma 2.1, we have

a 3 (t) L ∞ (∂ω * ) = J(t) L ∞ (∂ω * ) ≤ C a 3 (t) -a 3 (0) L ∞ (∂ω * ) = J(t) -1 L ∞ (∂ω * ) ≤ C η. (4.8)
Let us now prove that, for all t ∈ [0, 1],

a 1 (t) L ∞ (∂ω * ) ≤ C and a 1 (t) -a 1 (0) L ∞ (∂ω * ) ≤ C w(η), (4.9)
where w is the modulus of continuity given in Proposition 4.3. Using the fact that v t is bounded in C 2,α (Ω\ω * ) (see Lemma 4.4), we have immediately

a 1 (t) L ∞ (∂ω * ) = t DΦ -1 t ∇ t DΦ -1 t D(v t ) 2 • n L ∞ (∂ω * ) ≤ C. Let us define ξ(t) := ∇ t DΦ -1 t D(v t ) 2 .
Inserting ξ(t), we obtain

a 1 (t) -a 1 (0) L ∞ (∂ω * ) = t DΦ -1 t -I ξ(t) • n + ξ(t) -∇ |D(u)| 2 • n L ∞ (∂ω * ) ≤ C ξ(t) L ∞ (∂ω * ) t DΦ -1 t -I L ∞ (∂ω * ) + C ξ(t)-∇ |D(u)| 2 L ∞ (∂ω * ) ≤ C t DΦ -1 t -I L ∞ (∂ω * ) + C ξ(t) -∇ |D(u)| 2 L ∞ (∂ω * ) , since ξ(t) L ∞ (∂ω * ) ≤ C using (2.
3). By Estimate (2.3), we have t DΦ -1 t -I L ∞ (∂ω * ) ≤ Cη and using the estimate on the norm of v t given by Proposition 4.3,

ξ(t) -∇ |D(u)| 2 L ∞ (∂ω * ) ≤ w(η).
Hence, we obtain the second estimate of (4.9). Finally, Lemma 2.1 ensures that

a 0 (t) L 2 (∂ω * ) = m Y Θ • Φ t L 2 (∂ω * ) ≤ C m Y Θ L 2 (∂ω * ) a 0 (t) -a 0 (0) L 2 (∂ω * ) ≤ C η m Y Θ L 2 (∂ω * ) . (4.10)
According to the expression (4.6) of I 1 and the definition of a i (i = 0, . . . , 3),

I 1 (t)-I 1 (0) = ∂ω * ν a 2 0 (t) -a 2 0 (0) 3 i=1 a i (t)+ ∂ω * ν a 2 0 (t) 3 i=1 a i (t) - 3 i=1 a i (0) . (4.11)
The last integral of (4.11) is bounded by

C m Y Θ 2 L 2 (∂ω * ) 3 i=1 a i (t) - 3 i=1 a i (0) L ∞ (∂ω * )
.

Using the estimates on the norms of a i (i = 1, . . . , 3) given by (4.7), (4.8) and (4.9) and the fact that w(η) ≥ C η by Assumption (4.4), we have

3 i=1 a i (t) - 3 i=1 a i (0) L ∞ (∂ω * ) ≤ C 3 i=1 a i (t) -a i (0) L ∞ (∂ω * ) ≤ C(w(η) + η) ≤ C w(η).
Moreover, the first integral of (4.11) is bounded by

3 i=1 a i (t) L ∞ (∂ω * ) a 0 (t) -a 0 (0) L 2 (∂ω * ) a 0 (t) + a 0 (0) L 2 (∂ω * )
and using the estimates on the norms of a i (i = 0, . . . , 3) given by (4.7), (4.8), (4.9) and (4.10) and the remark that w(η) ≥ C η, it is bounded by

C w(η) m Y Θ 2 L 2 (∂ω * ) ,
which concludes the proof.

Proof of Lemma 4.6

Since

u ′ t = -∂ nt u t (Y Θ • n t ) on ∂ω t , u ′ t = -m Y Θ ∂ nt u t (n t • n) on ∂ω t . (4.12)
Let us define the following Steklov-Poincaré operator corresponding to D N (introduced in Lemma 3.1) on ∂ω t :

C t : H 1/2 (∂ω t ) -→ H -1/2 (∂ω t ) z t -→ C t (z t ) := σ(Z t , Π t )n t ,
where

(Z t , Π t ) ∈ H 1 (Ω\ω t ) × L 2 (Ω\ω t ) solves        -div (σ(Z t , Π t )) = 0 in Ω\ω t div Z t = 0 in Ω\ω t Z t = 0 on ∂Ω Z t = z t on ∂ω t . (4.13)
We obviously have

1 2 ν D(Z t ) 2 L 2 (Ω\ωt) = C t (z t ) , z t ∂ωt . (4.14)
Let us define z t := -m Y Θ ∂ nt u t (n t • n). Thus, according to the boundary condition of u ′ t (see (4.12)), σ(u ′ t , p ′ t )n t = C t (z t ) and, using the definition of I 2 , we obtain

I 2 (t) = 2 ∂ωt ∂ nt u t C t (z t )(Y Θ • n t ).
Hence, using Estimate (4.14) on the norm of D(•), we get

I 2 (t) = -2 ∂ωt z t C t (z t ) = -ν D(Z t ) 2 L 2 (Ω\ωt) ,
where (Z t , Π t ) ∈ H 1 (Ω\ω t ) × L 2 (Ω\ω t ) solves (4.13).

We want now to estimate

I 2 (t) -I 2 (0) = Ω\ω * ν |D(Z 0 )| 2 - Ω\ωt ν |D(Z t )| 2 = Ω\ω * ν |D(Z 0 )| 2 - Ω\ω * ν t DΦ -1 t D( Z t ) 2 Jac t ,
where Z t := Z t • Φ t . Since div Y Θ = 0, we get Jac t = 1. We assume for a while that the following lemma is proved:

Lemma 4.7. For all t ∈ [0, 1],

z t H 1/2 (∂ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) , z t -z 0 H 1/2 (∂ω * ) ≤ Cw(η) m Y Θ H 1/2 (∂ω * ) , Z t H 1 (Ω\ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) and Z t -Z 0 H 1 (Ω\ω * ) ≤ Cw(η) m Y Θ H 1/2 (∂ω * ) .
Using this lemma and the fact that t DΦ -1

t L ∞ (R d ) ≤ C and t DΦ -1 t -I L ∞ (R d ) ≤ Cη (see Estimate (2.3)), we have |I 2 (t) -I 2 (0)| = Ω\ω * ν |D(Z 0 )| 2 - Ω\ω * ν t DΦ -1 t D( Z t ) 2 ≤ C D(Z 0 ) + t DΦ -1 t D( Z t ) L 2 (Ω\ω * ) D(Z 0 ) -t DΦ -1 t D( Z t ) L 2 (Ω\ω * ) ≤ C D(Z 0 ) L 2 (Ω\ω * ) + t DΦ -1 t D( Z t ) 2 (Ω\ω * ) t DΦ -1 t D( Z t ) -t DΦ -1 t D(Z 0 ) L 2 (Ω\ω * ) + t DΦ -1 t D(Z 0 ) -D(Z 0 ) L 2 (Ω\ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) t DΦ -1 t D( Z t ) -t DΦ -1 t D(Z 0 ) L 2 (Ω\ω * ) + t DΦ -1 t D(Z 0 ) -D(Z 0 ) L 2 (Ω\ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) D( Z t ) -D(Z 0 ) L 2 (Ω\ω * ) + Cη D(Z 0 ) L 2 (Ω\ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) w(η) m Y Θ H 1/2 (∂ω * ) + η m Y Θ H 1/2 (∂ω * ) .
The proof is completed. It remains to prove Lemma 4.7.

Proof of Lemma 4.7. We recall that

z t := -m Y Θ ∂ nt u t (n t • n) and z t := z t • Φ t .
We use the following product lemma:

∀v ∈ H 1/2 (∂ω * ), w ∈ C 1 (∂ω * ), v w H 1/2 (∂ω * ) ≤ C v H 1/2 (∂ω * ) w C 1 (∂ω * )
for some constant C depending only on Ω\ω to get

z t H 1/2 (∂ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) ∂ nt u t (n t • n) C 1 (∂ω * ) .
Then, using the estimate on the norm of n t given by Lemma 2.1 and the estimate on the norm of v t given by Proposition 4.3, we obtain the first estimate of the lemma:

z t H 1/2 (∂ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) . (4.15)
Defining n t := n t • Φ t , we have

z t -z 0 = z t -m Y Θ ∂ n u = m Y Θ (∂ n u -∂ nt u t • Φ t ( n t • n)) on ∂ω * .
However,

∂ n u -∂ nt u t • Φ t ( n t • n) C 1 (∂ω * ) ≤ ∂ n u (1 -( n t • n)) C 1 (∂ω * ) + ( n t • n) (∂ n u -∂ nt u t • Φ t ) C 1 (∂ω * ) .
Moreover, using the estimate on the norm of n t given by Lemma 2.1, we have

∂ n u (1 -( n t • n)) C 1 (∂ω * ) ≤ Cη, since ∂ n u C 1 (∂ω * ) ≤ C. Now, let us bound ( n t • n) (∂ n u -∂ nt u t • Φ t ) C 1 (
( n t • n) (∂ n u -∂ nt u t • Φ t ) C 1 (∂ω * ) = ( n t • n) ∂ n u -t DΦ -1 t ∇v t n t C 1 (∂ω * ) ≤ C ∇u(n t -n) C 1 (∂ω * ) + ∇u -t DΦ -1 t ∇v t n t C 1 (∂ω * ) ≤ Cw(η),
since we can assume η ≤ Cw(η) (see (4.4)). the second inequality of the lemma is proved:

z t -z 0 H 1/2 (∂ω * ) ≤ Cw(η) m Y Θ H 1/2 (∂ω * ) . (4.16)
Since Z 0 solves the Stokes system with z 0 as interior boundary condition, the following energy estimate holds using the classical energy estimate concerning the Stokes system (see for example the books [START_REF] Galdi | An introduction to the mathematical theory of the Navier-Stokes equations[END_REF][START_REF] Fabre | Prolongement unique des solutions de l'equation de Stokes[END_REF]) and the first estimate (4.15)

Z 0 H 1 (Ω\ω * ) + Π 0 L 2 (Ω\ω * ) ≤ C z 0 H 1/2 (∂ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) . (4.17)
With this inequality, it remains to prove the inequality concerning Z t -Z 0

H 1 (Ω\ω * )
. Indeed, we will then have

Z t H 1 (Ω\ω * ) ≤ Cw(η) m Y Θ H 1/2 (∂ω * ) + Z 0 H 1 (Ω\ω * ) ≤ C m Y Θ H 1/2 (∂ω * ) .
The couple (Z t , Π t ) ∈ H 1 (Ω\ω t ) × L 2 (Ω\ω t ) solves (4.13). Thus, defining Z t := Z t • Φ t and Π t := Π • Φ t , we prove that

         -div ∇ Z t A(t) + div Π t B(t) = 0 in Ω\ω * ∇ Z t : B(t) = 0 in Ω\ω * Z t = 0 on ∂Ω Z t = z t on ∂ω * .
Using the problem solved by (Z 0 , Π 0 ), this problem can be rewritten as

           ν∆( Z t -Z 0 ) -∇( Π t -Π 0 ) = L 1 (t) Z t , Π t -S Z t , Π t in Ω\ω * -div Z t -Z 0 = L 2 (t)( Z t ) -div Z t in Ω\ω * Z t -Z 0 = 0 on ∂Ω Z t -Z 0 = z t -z 0 on ∂ω * , with L 1 (t) Z t , Π t := -div ∇ Z t A(t) + div Π t B(t) L 2 (t)( Z t ) := ∇ Z t : B(t) S Z t , Π t := -ν∆( Z t ) + ∇( Π t ).
Thus, using the energy estimate,

Z t -Z 0 H 1 (Ω\ω * ) + Π t -Π 0 L 2 (Ω\ω * ) ≤ C z t -z 0 H 1/2 (∂ω * ) (L 1 (t) -S) Z t , Π t H -1 (Ω\ω * ) + (L 2 (t) -div ) Z t ) L 2 (Ω\ω * )
. (4.18)

We have

(L 1 (t) -S) Z t , Π t = (L 1 (t) -S) Z t -Z 0 , Π t -Π 0 + (L 1 (t) -S) Z 0 , Π 0 (L 2 (t) -div ) Z t ) = (L 2 (t) -div ) Z t -Z 0 + 2 (t) -div ) Z 0 ) .
Let us assume for a while where Θ is defined in subsection 2.2. We have shown the precised estimate: Let us make a final comment: this result can be understood as a quantitative estimate of the deviation from minimality for sets close to ω * in a strong sense (here the C 2,α norm). Similar results have been introduced recently. The present work corresponds to the first step in extending to the drag functional ideas and result on isoperimetric problems, see [START_REF] Acerbi | Minimality via second variation for a nonlocal isoperimetric problem[END_REF]. In a second step, using additional regularity properties provided by the perimeter, Acerbi et al. manage to extend this result on smooth domains to less regular sets.

(L 1 (t) -S) (w, q) H -1 (Ω\ω * ) ≤ Cη w H 1 (Ω\ω * ) + q L 2 (Ω\ω * ) (L 2 (t) -div ) (w) L 2 (

Conclusion

This paper is, in some sense, a generalization to a vectorial of some previous work (see [START_REF] Dambrine | About stability of equilibrium shapes[END_REF][START_REF] Dambrine | On variations of the shape Hessian and sufficient conditions for the stability of critical shapes[END_REF]) where the stability of critical shapes is investigated for the scalar case. Indeed, we here focused on the minimization problem of the energy dissipated by a newtonian and incompressible fluid driven by the Stokes law. Then, we gave a sufficient condition for the shape Hessian to be coercive in a weaker norm than the differentiability norm. However, we showed that this coercivity is sufficient to prove that such a critical shape is a local minimizer. Finally, we obtained a lower bound of the variations of the drag in terms of the symmetric difference of the domain.

This paper may be extended to the case where the fluid motion is assumed to be governed by the (stationary) Navier-Stokes equations. The most difficult part in order to obtain analogous results seems to be the coercivity of the shape Hessian. Moreover, the present work can be adapted to other classical boundary conditions. Finally, it should be possible to adapt these results to unbounded domain. However, due to the approach using derivatives, removing the regularity assumptions on the boundary and the deformations seems out of reach without new ideas.

Lemma 3 . 2 .

 32 It holds ν∇ |D(u)| 2 • n = 4 ν(d -1)HΛ + ν t b F b + ∇ τ p • b on ∂ω * ,(3.11)where• b is the projection of ∂ n u on the tangential space of ∂ω * ;

Proposition 4 . 3 .

 43 There exists a function w : [0, 1] → R with lim rց0 w(r) = 0 such that sup t∈[0,1]

Proof of Proposition 4 . 3 .

 43 We use similar arguments than the ones used in [10, Proof of Proposition 4.1]. The main ingredient is the compact embedding of C 2,α (Ω\ω * ) into C 2 (Ω\ω * ) and of C 1,α (Ω\ω * ) into C 1 (Ω\ω * ). Let η 2 > 0 given by Lemma 4.4. Define, for η ∈ (0, η 2 ),

Theorem 5 . 2 .

 52 If ω * is a critical shape for J where (4.1) holds, there exists η > 0 and C > 0 depending only on Ω and ω * such that, for all Θ ∈ U with Θ -I C 2,1 (R d ) < η, |Θ(ω * )| = |ω * | and Θ = I, J(Θ(ω * )) > J(ω * ) + C L d (ω * ∆Θ(ω * )) 2 .

  Study of J 1 . We first introduce the Steklov-Poincaré operator of Displacement-to-Traction also called Dirichlet-to-Neumann operator and prove that this operator is positive. Lemma 3.1. We consider Ω and ω two Lipschitz open sets of R d (d = 2, 3) such that ω ⊂⊂ Ω and Ω\ω is connected. For h ∈ H 1/2 (∂ω), let us define the following Steklov-Poincaré operator:

  ∂ω * ) . Using again Lemma 2.1, Estimate (2.3) and Proposition 4.3, we obtain

  Ω\ω * ) ≤ Cη w H 1 (Ω\ω * ) .Then, using this inequality in Estimate (4.18) and Estimates (4.16) and (4.17) on the norms of z tz 0 and Z 0 , we obtainZ L 2 (Ω\ω * ) w H 1 (Ω\ω * ) + q L 2 (Ω\ω * ) B(t) L 2 (Ω\ω * ) + w H 1 (Ω\ω * ) + q L 2 (Ω\ω * ) ≤ C Φ H 1 0 (Ω\ω * ) η w H 1 (Ω\ω * ) + q L 2 (Ω\ω * )using Lemma 2.1 and Estimate (2.3). We proceed in the same way to prove that, for all ξ ∈ L 2 (Ω\ω * ), that our previous stability estimates (4.2) provides the existence of a non negative real η and of a constant C such that for any diffeomorphism Θ ∈ U with Θ -I C 2,1 (R d ) < η and |Θ(Ω)| = |Ω|, it holds J(ω) ≥ J(ω * ) + C 4 Y Θ • n 2 H 1/2 (∂ω * ) ≥ J(ω * ) + C L d (ω * ∆Θ(ω * )) 2 ,

	(4.19)
	Then,
	Ω\ω

(L 1 (t) -S) Z t , Π t H -1 (Ω\ω * ) + (L 2 (t) -div ) Z t L 2 (Ω\ω * ) ≤ Cη Z t -Z 0 H 1 (Ω\ω * ) + Π t -Π 0 L 2 (Ω\ω * ) + Cη Z 0 H 1 (Ω\ω * ) . t -Z 0 H 1 (Ω\ω * ) (1 -Cη) ≤ Cw(η) V H 1/2 (∂ω * ) ,

which concludes the proof. It remains to prove Estimate

(4.19)

. Let us consider Φ ∈ H 1 0 (Ω\ω * ). We have

Ω\ω * (L 1 (t) -S) (w, q) • Φ = Ω\ω * (-div (∇wA(t)) + div (qB(t)) + ν∆w -∇q) • Φ = Ω\ω * ((∇wA(t)) + qB(t) + ν∇w -qI) : ∇Φ ≤ C Φ H 1 0 (Ω\ω * ) A(t) * (L 2 (t) -div ) (w) ξ ≤ C ξ L 2 (Ω\ω * ) η w H 1 (Ω\ω * ) ,

which concludes the proof of Estimate

(4.19)

. so

Precised stability estimate

The main result of this section is a purely geometrical result valid for small smooth perturbations of a domain.

Proposition 5.1. In the assumptions of Section 2, there is a constant C that depends only on Ω such that

Proof. We first consider the case of a smooth positive deformation function h and prove (5.1) in this first case. Since h is positive, Ω is a subset of Ω and hence

To this function, is associated the vector field V by the construction (2.1) of subsection 2.2 (in general or by (2.2) if a volume constraint has to satisfied). We then apply the classical Lemma about variations of the volume

to get along the flow of the field V :

where J t is the surface jacobian. We estimate the middle term

thanks to Lemma 2.1 and conclude. The same proof works for smooth negative deformation field, then by density to any signed Lipschitz deformation. For a general deformation, expressed with a smooth normal deformation h, we split it into h = h + -h -its positive and negative part and define V + and V -by the construction (2.1). Then it suffices to notice that Ω per = Φ 1 (V -)[Φ 1 (V + )[Ω]]: the idea is to first apply the field V + to match the inflating areas then the field V -to contract the domain. Now, in the context of our problem, Cauchy Schwarz inequality provides


Fleur Rennesson

Frédéric Longuetaud

Mothe

A Krähenbühl

B Kerautret

I Debled-Rennesson

F Mothe

Knot Detection in

published or not. The documents may come

Introduction

Detect and identify automatically dig- ital objects in 3D volumetric images is a challenge, especially in medical imaging applications. Various methods are proposed to extract and measure internal characteristics of the human, like brain or blood vessel. Currently, the most effective techniques use deformable models [START_REF] Mcinerney | Deformable models in medical image analysis: A survey[END_REF].

In this work, we focus on a similar problem with 3D digital images of wood. The images are obtained by medical X-Ray Computer Tomography (CT) scanners and are illustrated on Fig. 1. On Fig. 1(a), a raw image is projected in an intensity interval which highlight four knots starting from the pith 3 . The related volume obtained by thresholding is rendered on Fig. 1(b).

A knot as in Fig. 1 and Fig. 2 is the first part of a branch, included within the tree stem due to the radial growth of the stem. The frequency and size of the knots are the first depreciation factors considered by wood suppliers for estimating the price of timber. Knottiness is also one of the main criteria considered in the visual grading of lumber. Knowledge of knot geometry and location would be valuable at sawmills for optimising cutting decisions or improving the grading of logs. X-ray CT has been recognised as being the most promising method to non-destructively analyse the internal structure of wood [START_REF] Schad | Nondestructive methods for detecting defects in softwood logs[END_REF]. CT scanners designed expressly for the wood industry are now available and some of the largest sawmills are now equipped with them. Moreover, such detailed data about knottiness are needed in the field of forest research to study for example the relationship between tree growth conditions and wood quality.

Several methods were proposed for measuring knottiness on the basis of CT images (see [3] for a review). But no method is entirely automatic and well adapted to all species of wood and to wood containing wet areas such as sapwood4 We proposed in [3] an entirely automated method to identify knots from X-ray CT images of a piece of wood (log or beam) and obtain data on knot geometry without any human intervention. However, the first step of the algorithm, image segmentation, was not studied in details. A simple thresholding operation was sufficient enough for processing the air-dried wood that were used in this study. For an application in industrial conditions, there is a need to improve the segmentation algorithm, especially for wet wood containing areas of sapwood (see Fig. 2(c)).

The main contribution of this paper is to automatically and quickly detect all knots of a given log based on the analysis of X-ray CT images. Various softwood species are studied, various moisture contents were considered leading to have large sapwood areas visible in the CT images (see Fig. 2(c)), and our method identifies all knots in all cases.

The proposed knot detection is based on a very precise detection of the localisation of each knot. To do this, we use intensity variation between two consecutive slices and study a specific intensity histogram. We proceed in two steps. The first one is to select slice intervals where intensity variations are important around the pith. In the second step, we keep a set of angular sectors also centered on the pith in each interval from the same criterion of intensity variation.

In the next section we list all tools needed by the proposed method. Section 3 introduces the definition of the cumulative z-motion histogram which allows to identify the meaning slice intervals and to deduce the angular sector of each knot. The last section shows results and presents perspectives of this work.

2 Suitable tools for X-Ray CT images of logs

Pith detection

Pith detection is an essential first step to wood knot detection. In particular, in the proposed method it will allow to focus on a restricted area around the pith. In this area, we are sure that intensity variations observed in images indicate the presence of knots. Indeed, the knots are connected to the pith where they have their origin. They are oriented radially in a log, from the pith to the bark, and generally upward. Some knots are totally included into the log and covered by the radial growth rings5 but in general, thay are still visible at the outside of the log in the form of branches.

To detect pith, we use the algorithm proposed by Fleur Longuetaud [START_REF] Longuetaud | Automatic detection of pith on ct images of spruce logs[END_REF] and improved by Boukadida et al. [START_REF] Boukadida | Pithextract: a robust algorithm for pith detection in computer tomography images of wood -application to 125 logs from 17 tree species[END_REF]. The idea is to detect the pith successively in each slice. Its location in slice n allows to steer and speed up pith computing on slice n + 1. The first step of the algorithm uses the roundness of growth rings. They are detected in each slice I k with a Sobel filter. This one computes G x and G y images, respectively the gradients of I k along x and y.

A new image I ′ k is computed with a thresholding of G norm (see Fig. 3). I ′ k allows to define a set of lines. For that, we are computing Θ, the direction of gradient in each pixel of I ′ k not null. After, we are building the map of the intersection locations for all lines taken two by two. The intersection locations are computed with an Hough transform [START_REF] Duda | Use of the hough transformation to detect lines and curves in pictures[END_REF]. The pith position is then defined like the location of the largest number of intersecting lines.

In a second step, the pith position is corrected in slices where pith neighbourhood is noisy. It is usually the slices with knots. These ones disturb the roudness of growth rings and this noise hugely affects the gradient direction and therefore Hough transform. In these slices, the pith is recomputed by linear interpolation of the previous and the next slice.

Intervals detection in an histogram

The method that we propose detects intervals centered on maximums. It is a problem that can be very difficult to solve according to noise level on histogram. The detection that we are now presenting is suitable for the histograms we have to deal with.

Let H be an

N -class histogram with H = (H k ) k∈[1,N ] .
Smoothing The first step consists to smooth H with an averaging filter of radius r f . The aim is to improve the computing of the discrete derivative during the following step. We get H l histogram where the k class is :

H l k = 1 2r f + 1 k+r f i=k-r f H i with r f : radius of histogram averaging mask.

Maximums computing

We defined M as the set of the index of H l maximums from a neighbourhood of r M radius.

M = {k|H l k > t M and ∀i ∈ [-r M , -1] ∪ [1, r M ], H l k > H l k+i } with t M : minimum value for one maximum.
r M : width of neighbourhood on which the maximums are defined.

Value of a maximum must be greater than t M . This allows to avoid the small maximums coming from the noise. And a maximum must have a value greater than the ones of its r M 2 left and right neighbours.

Intervals computing Let V be the set of searched intervals. We determine V from the M set of maximum indexes. An interval of classes is computed from each maximum, starting from each side of this maximum.

V = { [k -i, k + j] | k ∈ M and ∃ t i ∈ [k -i, k] such that    ∀t ∈ [t i , k], H l t > t V ∀t ∈]k -i, t i [, H l t ≤ t V and (H l t ) ′ > 0 (H l k-i ) ′ ≤ 0 and ∃ t j ∈ [k + 1, k + j] such that    ∀t ∈]k, t j ], H l t > t V ∀t ∈]t j , k + j[, H l t ≤ t V and (H l t ) ′ < 0 (H l k+j ) ′ ≥ 0 and j -i > w V } with (H l t ) ′ = H l t -H l t-1
and t V : derivative threshold. w V : minimum width of an interval.

Derivatives are tested on each side of maximums only, from the nearest neighbour with a value lower than t V . An interval is kept if its magnitude is larger than w V .

Z-motion histogram and knot detection

With detection methods of pith and intervals, presented in previous section, we can determine knot areas. They take the form of angular sectors of cylinder and our goal is to get one area by knot. Let us specify that we detect these areas in image sequences from medical X-Ray scanner. Each image corresponds to a circular slice of log (see Fig. 2). Pixel values corresponds to density on Hounsfield scale. This density is linearly correlated to the wood density. We referred either to wood density or to images intensity.

On X-Ray images of log, pixel intensity is approximately between -3000 and 3000. These bounds changed few from an image to an other. On Hounsfield scale, wood density is usually greater than -900 and always lesser than 530. All image values out of this interval will be considered as outlier and ignored.

Our method consists in applying twice the process of interval detection. The two processed histograms will be constructed from the same notion : the zmotion.

The z-motion We decided to use the z-motion by observing successively an image sequence. By scrolling an image sequence of one log, egg-shaped shapes are observed in motion. These movements occur from pith to bark, orthogonally to the growth rings, when the images are scrolled from the bottom to the top of the tree. These are the knot movements that have their origin at the pith location and have grown radially and upwards.

The z-motion Z of a slice I k is defined like the intensity variation between the two successive slices I k-1 and I k :

Z k = |I k -I k-1 |
With help of this notion, we now present the two steps of our method. Before this two steps, it is necessary to detect pith (see Sec. 2.1) on all slices. It allows to define a C circle centered on pith in each slice of log. Then for the first step, the histogram is constructed by cumulating the z-motion of each slice on the entire circle C. We then obtain the H slices histogram on which we apply the process of interval detection (see Sec. 2.2). During the second step, we subdivided the C circle in a set of angular sectors with same size. For each interval of H slices , we combine the z-motion separately on each angular sector.

Hereafter are more details about the construction of these two histograms.

The histogram of cumulative z-motion

The slice histogram Let I be a set of images of a log slice. Let C be a circle of radius r C . We defined the H slices (k) value of H slices histogram for the I k slice as follows :

H slices (k) = (i,j)∈C Z k (i, j)
where -C is centered on the pith position in I k slice.

-∀(i, j), I k (i, j)

∈ [i min , i max ] et I k-1 (i, j) ∈ [i min , i max ] -∀(i, j), Z k (i, j) ∈ [z min , z max ]
with r C : radius of the cumulative z-motion circle. [i min , i max ]: intensity interval. [z min , z max ]: interval of intensity variations.

The z-motion is only computed inside the C circle centered on pith. According to experts, the minimum wood density can be fixed to -900. It is our value for i min . For i max , we fixed its value to 530, a biological constant for the maximum wood density. Variations lesser than z min due to growth rings, sapwood, etc. are considered as noise. Variations in [z min , z max ] interval are potential movements of knots. Variations greater than z max are considered as outliers.

In this first step, each histogram value corresponds to the z-motion sum in a slice. Intervals are therefore slice intervals where the cumulative z-motion is important (see Fig. 5).

One of the most important parameters is the r C radius of C circle. It is determine on each slice. We start from pith in ten directions as long as pixel values are greater than -900. In fact, values of X-ray scanned images fall drastically below -900 when we leave the log. Then we compute the average radius r and r C is fixed to 0.8 • r.

The z-motion interval [z min , z max ] is fixed to [200,500] to just cumulate important z-motion. After histogram construction, it is smoothed with an averaging mask where r f = 2. Each histogram value becomes candidate to the maximum set if its value is greater than t M = 5% of the maximum value of H slices .

To compute intervals, we fix the t V threshold to the same value as t M . All the first neighbors of a maximum greater than t V belong to the interval of this maximum. From the first neighbor lower than t V , following neighbors are added to the initial interval until the discrete derivative goes through 0. Finally we just keep slice intervals with a width greater than w V = 10. We obtain a set of slice intervals V H slices with all knots.

The histogram of angular sectors This second step is applied on each slice interval obtained at previous step. For a given slice interval, we cumulate zmotion separately in each angular sector. We obtain the histogram of cumulative intensity by angular sector (see Fig. 6). Angular sectors are a partition of C circle in N sectors with a same angle. We choose N = 360.

Let [a, b] ∈ V slices be an interval of slice indexes. Let C be a circle with a r C radius split in 360 sectors with same angle. We defined the value H sectors (α) of H sectors histogram restricted to angle of α index as follows :

H sectors (α) = k∈[a,b] (i,j)∈Cα Z k (i, j) with -C centered on pith coordinates in I k -C α the C angular sector of α index ∈ [0, 359] -∀k, ∀(i, j), I k (i, j) > i min et I k-1 (i, j) > i min -∀k, ∀(i, j), Z k (i, j) ∈ [z min , z max ]
For an interval V q of V slices , we obtain an histogram H sector of z-motion distribution according to angular sector. Like in the first step, we use algorithm of interval detection to obtain a set of angular intervals which contains one knot. We keep the same parameters as at the first step, t M and t V excepted. In fact, the majority of parameters was fixed in studying scanned picture properties and in this second step, images are the same.

However, t M and t V correspond to histogram properties. On histograms of angular sectors, peaks are usually the same size, unlike the slice histograms. In fact, the knots of a slice interval have about about same volume. We fix rationally t M and t V at 50% of the H sectors maximum value without risking to forget a maximum and therefore a knot sector.

Results

The presented method is able to quickly isolate each log knot in angular sectors. The idea is to select an efficient segmentation method without introducing a long and complex processing step. It is important for sawmills that hope to scan logs at 2 meters per second. For biologists, the most important criterion is the precision. They want to obtain accurate measurements such as length, inclination, maximal diameter and volume of knots. The biologists and sawmills are both very interested by the exact number of knots in a log.

To fulfil this dual purpose, we have chosen to use the extraction of connected components in each angular sector based on a simple threshold. The connected components of one slice are first computed and are then merged with the ones of the previous slice.

The Fig. 7 presents several results of the connected component extraction. The main connected component of one slice is represented in pink to the left images and the corresponding 3D reconstruction is given in green on the right image. On logs without sapwood (see Fig. 7(a)), we obtain all knots with an high precision. In contrast, when logs contain sapwood (see Fig. 7(b)), the connected component extraction produces a component where the knot is connected to sapwood and growth rings. This happens when the density of sapwood and growth rings is close to the knot. Thresholding before extraction does not allow to remove them whithout distort the knot. The first parameter that our algorithm allow to compute, even with sapwood, is the number of knots. It is the number of detected angular sectors. We can also determine the inclination of all knots by extracting angle corresponding to the maximum of histogram for each angular sector (see Fig. 6). The length is determined in logs without sapwood. When there are sapwood and when a knot is connected to it, experts estimate that node leaves the stem and length can be estimate by the radius of log. The diameter can be exactly computed at any position in the 3D reconstruction of knots without sapwood. In the others cases, we can compute the maximum diameter because it is computed from the middle of knot and sapwood is not present at this place. Without sapwood, the volume can also be determined from the 3D reconstruction of knots. The sapwood induce an overestimation due to the set of parts connected to the knot.

To compare our method with other usual approaches of deformable model, we experiment one of them able to deal with automated topology changes [START_REF] Lachaud | Deformable model with adaptive mesh and automated topology changes[END_REF]. For this purpose, the initial model was created by using a cylindrical mesh defined from the outer shape contour (see Fig. 8(a)). Note that we use the euclidean metric model implementation. The model is efficient to extract knots in logs without sapwood (See Fig. 8(a-c)). However it is not able to isolate the knot parts from sapwood areas (See Fig. 8(d)). If the segmentation results are comparable for the logs without sapwood, the real advantage of our method is the computation time. We use histogram computation and extraction of connected components. The knot detection time is about the second with our approach whereas deformable model takes really more time, more than one hour. The proposed method is also able to detect and analyze the knots for the case of sapwood while the deformable fails to isolate such a part in this case. The source code of the proposed method is available online [START_REF] Krähenbühl | Tkdetection[END_REF].

Conclusion

The method proposed in this paper automatically detects and counts all the knots of a given piece of wood based on the analysis of X-Ray CT images. It allows to easily deduce several measures such as inclination, length, maximal diameter and volume of knots. All its steps are based on elementary operations such as histogram computation and therefore the method is very fast. In future works, we will improve our method with objective to free ourselves of sapwood.

Fig. 1 .

 1 Fig. 1. Example of wood image and 3D volumic rendering.

  Fig. 2. X-Ray CT images of stem : (a) 3D volumetric rendering , (b) dry and (c) with sapwood.

Fig. 3 .

 3 Fig. 3. Pith detection in two slices, with (below) and without (above) knot.

Fig. 4 .

 4 Fig. 4. Angular sesctors intervals on a given slice interval. The green sections represents the slice and the angular sector where the z-motion cumulation is maximum. We can see in initial log section the considered circle C and the localisation of pith.

Fig. 5 .

 5 Fig. 5. Histogram after detection of slice intervals in a spruce log.

Fig. 6 .

 6 Fig. 6. Histogram and circular projection after detection of sector interval.

Fig. 7 .

 7 Fig. 7. Detected angular sectors and connex components with 3D visualization.

Fig. 8 .

 8 Fig.8. Experiments of deformable models[START_REF] Lachaud | Deformable model with adaptive mesh and automated topology changes[END_REF] on wood with and without sapwood.

pith : central part of stem.

sapwood : wood area located at the periphery of the stem.

growth rings: concentric circles in stem. They are centered on pith.

Acknowledgements

The authors would like to thank Jacques-Olivier Lachaud and Benjamin Taton for kindly providing the source code of the deformable model [7].