On the detection of several obstacles in 2D Stokes flow: topological sensitivity and combination with shape derivatives

We consider the inverse problem of detecting the location and the shape of several obstacles immersed in a ﬂuid ﬂowing in a larger bounded domain Ω from partial boundary measurements in the two dimensional case. The ﬂuid ﬂow is governed by the steady-state Stokes equations. We use a topological sensitivity analysis for the Kohn-Vogelius functional in order to ﬁnd the number and the qualitative location of the objects. Then we explore the numerical possibilities of this approach and also present a numerical method which combines the topological gradient algorithm with the classical geometric shape gradient algorithm; this blending method allows to ﬁnd the number of objects, their relative location and their approximate shape.


Introduction
We deal in this work with the inverse problem of determining the number, the position and the shape of relatively small objects inside a two dimensional fluid. We assume that the fluid motion is governed by the steady-state Stokes equations. In order to reconstruct the obstacles, we assume that a Cauchy pair is given on a part of the surface of the fluid, that is a Dirichlet boundary condition and the measurement of the Cauchy forces. Hence, the identifiability result of Alvarez et al. [5,Theorem 1.2] implies that this problem could be seen as the minimisation of a cost functional, which in our case will be a Kohn-Vogelius type cost functional.
The small size assumption on the objects leads us to perform asymptotic expansions on the involved functional. For this, we will use the notion of topological gradient which will determine a criteria in order to minimise the cost functional. The topological sensitivity analysis consists in studying the variation of a cost functional with respect to the modification of the topology of the domain, for example when we insert 'holes' (or objects) in the domain. It was introduced by Schumacher in [43] and Sokolowski et al. in [47] for the compliance minimisation in linear elasticity.
Topological sensitive analysis related to Stokes equations have been studied in the past by several authors, especially relevant are the works of Guillaume et al. [29], Maatoug [30], Amstutz [6] with steady-state Navier-Stokes equations and [7] with generalization for some non-linear systems and Sid Idris [44] which develops a detailed work in the two-dimensional case. In all of these works the focus is set to find topological asymptotic expansions for a general class of functionals where the system satisfies only Dirichlet boundary conditions. Closer works to our problem have been presented in the past by Ben Abda et al. [10] and by Caubet et al. [20]. In the first reference they consider a Neumann boundary condition on the small objects obtaining general results in two and three dimensional cases, with a complete development of the theory only on the three dimensional case. In the work of Caubet et al., they deal with the same problem as the one we consider here but only again on the three-dimensional case. In our two-dimensional case, due to the impossibility to have an asymptotic expansion of the solution of Stokes equations by means of an exterior problem (phenomena which is related to the Stokes paradox), we have to approximate it by means of a different problem. The deduction of this approximation is strongly influenced by the recent work of Bonnaillie-Noël et al. [11]. Indeed the same problem appears for the Laplace equation: it is based on the fact that the existence of a solution of the boundary value problem is not guaranteed except when u 0 (z) = 0. The classical analysis of elliptic equation in unbounded domain is made in the functional setting of weighted Sobolev spaces. It is known that (1.1) has a unique solution in a space containing the constants, hence this solution is the constant u 0 (z) which prohibits the condition at infinity if u 0 (z) = 0. Taking into account of this, we can define the asymptotic expansion for the Stokes system which is a crucial part in order to obtain the desired expansion for the functional involved. It is important to remark that (for a given real number u 0 (z)) several technical results which lead to the main result are different to the ones in the three-dimensional setting which involves additional difficulties to our problem. From the obtained theoretical results, we present some numerical simulations in order to confirm and deepen our theoretical results by testing the influence of some parameters in our algorithm of reconstruction such as the shape and the size of the obstacles. We also propose an algorithm which joins the topological optimization procedure with the classical shape optimization method using the previous computation of the shape gradient for the Kohn-Vogelius functional made by Caubet et al. in [21]. This blending method allows not only to obtain the number and qualitative location of the objects, moreover it allows to approximate the shape of this ones. Nevertheless, we precise that the geometrical shape optimization step will fail if the previous topological step doesn't give the total number of objects.
To conclude, we also mention the recent developments on topological sensitivity based iterative schemes made by Carpio et al. in [16,18,19]. We also refer to some works using the level set method by Lesselier et al. in [34,25,26]. Combinations of several shape optimization methods was also recently tested by several authors. Allaire et al. propose in [3] to couple the classical geometrical shape optimization through the level set method and the topological gradient in order to minimize the compliance. The same combination is made for another problem by He et al. in [32]. In [15], Burger et al. use also this combination for inverse problems. There, the topological gradient is incorporated as a source term in the transport Hamilton-Jacobi equation used in the level set method. Concerning the minimization of the compliance, Pantz et al. propose in [42] an algorithm using boundary variations, topological derivatives and homogenization methods (without a level set approach).
Motivations This obstacle inverse problem arises, for example, in mold filling during which small gas bubbles can be created and trapped inside the material (as it is mentioned in [10]) We can also mention the fact that the most common devices used to spot immersed bodies, such as submarines or banks of fish, are sonars, using acoustic waves: Active sonars emit acoustic waves (making themselves detectable), while passive sonars only listen (and can only detect targets that are noisy enough). To overcome those limitations, one want to design systems imitating the lateral line systems of fish, a sense organ they use to detect movement and vibration in the surrounding water (as emphasized in [23]).
Organization of the paper The paper is organized as follows. First, we introduce the adopted notations. Then, in Section 2, we present in details the considered problem and give the main idea used to study it: we introduce some perturbed domains and the considered Kohn-Vogelius functional. Section 3 is devoted to the statement of the main result: we give the topological asymptotic expansion of this functional. In Section 4, we prove the asymptotic expansion of the solution of the considered Stokes problems when we add small obstacles inside. Then, we use the resulting estimates to prove the main results by splitting the functional in Section 5. We make numerical attempts in Section 6. We explore the efficiency of this method and point out these limits. Finally we propose a new algorithm in Section 7 which combines the topological gradient techniques with shape gradient techniques in order to be able to find the number of objects, their approximate location and their approximate shape. Technical results needed to justify the expansions are postponed in appendices.

The problem setting
Introduction of the general notations. For a bounded Lipschitz open set Ω ⊂ R 2 , 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, we denote by W 1,p α (Ω) the weighted Sobolev spaces defined in Appendix B Definition B.1. For k ∈ N, we denote · k,Ω the norm · H k (Ω) and |·| 1,Ω the semi-norm of H 1 (Ω). We also use the notations · 1/2,∂Ω and · −1/2,∂Ω to define respectively the norms · H 1/2 (∂Ω) and · H −1/2 (∂Ω) . We represent the duality product between H −1/2 (∂Ω) and H 1/2 (∂Ω) using the notation · , · ∂Ω . Finally, n represents the external unit normal to ∂Ω and we define the space We here precise that the notation Ω means Ω p(x)dx which is the classical Lebesgue integral. Moreover, we use the notation ∂Ω p to denote the boundary integral ∂Ω p(x)ds(x), where ds represents the surface Lebesgue measure on the boundary. The aim is to simplify the notations when there is no confusion.

Framework.
Let Ω be a bounded Lipschitz open set of R 2 containing a Newtonian and incompressible fluid with coefficient of kinematic viscosity ν > 0. Let ω ⊂ R 2 a fixed bounded Lipschitz domain containing the origin. For z ∈ Ω and 0 < ε << 1, we denote ω z,ε := z + εω.
The aim of this work is to detect some unknown objects included in Ω. We assume that a finite number m * of obstacles ω * z,ε ⊂ Ω have to be detected. Moreover, we assume that they are well separated (that is: and have the geometry form where ε k is the diameter and ω * k ⊂ R 2 are bounded Lipschitz domains containing the origin. The points z * k ∈ Ω, 1 ≤ k ≤ m * , determine the location of the objects. Finally, we assume that, for all 1 ≤ k ≤ m * , ω * z k ,ε k is far from the boundary ∂Ω. In order to determine the location of the objects, we make a measurement g ∈ H −1/2 (O) on a part O of the exterior boundary ∂Ω with O ∂Ω. Then, we denote ω * ε := m * k=1 ω * z k ,ε k and consider the following overdetermined Stokes problem Here u represents the velocity of the fluid and p the pressure and σ(u, p) represents the stress tensor defined by σ(u, p) := ν ∇u + t ∇u − pI.
We assume here that there is no body force and consider the homogeneous Dirichlet boundary conditions on the obstacles, which is the so-called no-slip boundary conditions. Notice that, if div u = 0 in Ω, we have with D(u) := ∇u + t ∇u . Thus we consider the following geometric inverse problem: Find ω * ε ⊂⊂ Ω and a pair (u, p) which satisfy the overdetermined problem (2.2). (2.3) To study this inverse problem, we consider two forward problems: where ω ε := m k=1 ω z k ,ε k for a finite number m of objects located in z 1 , . . . , z m . These two forward problems are classically well-defined. We refer to [13,27] for the results of existence and uniqueness of (u ε D , p ε D ). Notice that the compatibility condition (2.1) associated with Problem (2.4) is satisfied. The existence and the uniqueness of (u ε M , p ε M ) is for example detailed in [20,Theorem A.1]. We underline the fact that p ε M does not need to be normalized to be unique due to the Neumann boundary conditions imposed on O.
One can remark that if ω ε coincides with the actual domain ω * ε , then u ε D = u ε M in Ω\ω ε . According to this observation, we propose a resolution of the inverse problem (2.3) of reconstructing ω * ε based on the minimization of the following Kohn-Vogelius functional We then define . We can notice that, integrating by parts the expression of F KV . This expression shows that the error can be expressed by an integral involving only the part of the boundary where we make the measurement and reveals the coupling of the solutions via this functional. Finally, we can notice that the Dirichlet error is weighted by the Neumann error, and vice versa.
Remark 2.1. In order to guarantee that the inverse problem of finding ω * ε and a pair (u, p) satisfying (2.2) has a solution, we have to assume the existence of such a ω * ε . This means that the measurement g is perfect, that is to say without error. Then, according to the identifiability result [4, Theorem 1.2] proved by Alvarez et al., the domain ω * ε is unique. Notice that in [4], ω * ε is assumed to have a C 1,1 boundary but we can only assume that it has a Lipschitz boundary in the Stokes case (see [9,Theorem 2.1]). Hence, if we find ω * ε such that J KV (Ω \ ω * ε ) = 0, then u ε D = u ε M in Ω \ ω * ε , i.e. u ε D satisfies (2.2) and thus ω ε = ω * ε is the real domain. In the following, for ε = 0, we will consider as a convention that ω 0 = ∅ (instead of ω 0 = m k=1 {z k }, which comes from the definition of ω ε ), and therefore: Ω 0 = Ω. Then, we will denote (u 0 (Ω) the respective solutions of the following systems:

The main result
From now on, we consider that we seek a single obstacle ω z,ε := z + εω, located at a point z ∈ Ω. Notice that in the case of several inclusions, we proceed by detecting the objects one by one. Thus, after detecting a first obstacle ω z 1 ,ε 1 , we work replacing the whole domain Ω by Ω\ω z 1 ,ω 1 (and then we have ∂ω z 1 ,ε 1 ⊂ ∂ (Ω\ω z 1 ,ω 1 ) \O) and the results presented below (in particular the topological derivative) are still valid for a new inclusion ω z,ε . Note that, the asymptotic expansion of the solution of elliptic boundary value problem in multiply perforated domains is studied in [12,36].

Introduction of the needed functional tools
We recall that the topological sensitivity analysis consists in the study of the variations of a design functional J with respect to the insertion of a small obstacle ω z,ε at the point z ∈ Ω (with no-slip boundary conditions). The aim is to obtain an asymptotic expansion of J of the form where ε > 0, ξ is a positive scalar function intended to tend to zero with ε and where Ω z,ε := Ω\ω z,ε , with ω z,ε := z + εω. We summarize the notations concerning the domains in Figure 1.
O Ω ∂Ω ω z,ε Ω z,ε The computation of the topological gradient δJ in this work is mainly based on the paper by Caubet and Dambrine [20] which deals with the presented problem in the threedimensional setting. The work of Bonnaillie-Noël and Dambrine [11], which deals with asymptotic expansions for Laplace equation in a domain with several obstacles, was the basis for the choice of the approximating problem in the two-dimensional setting. We also have been inspired strongly by the works of Sid Idris in [44] and [28,29] (written with Guillaume), where the authors study topological asymptotic expansions for Laplace and Stokes equations in two and three dimensions, which provides us several techniques specially useful for the technical proofs presented in the appendix. Finally, let us point out the works of Amstutz [6,7], where the author develops a topological asymptotic expansion for a cost functional in the context of a fluid governed by the stationary Navier-Stokes equations, which contribute to understand better the possibilities for the asymptotic expansion of the solutions for our considered systems. It is important to mention that in all these situations the problem involves only Dirichlet boundary conditions.
We recall the expression of the fundamental solution (E, P ) to the Stokes system in R 2 given by is the canonical basis of R 2 and δ is the Dirac distribution.

The result
The following theorem gives us the expression of the topological gradient of the Kohn-Vogelius functional J KV : For z ∈ Ω, the functional J KV admits the following topological asymptotic expansion where u 0 D ∈ H 1 (Ω) and u 0 M ∈ H 1 (Ω) solve respectively Problems (2.4) and (2.5) with ω ε = ∅ and o(f (ε)) is the set of functions g(ε) such that lim ε→0 in the general asymptotic expansion (3.1).
Remark 3.2. Notice that, contrary to the 3 dimensional case [20, Theorem 3.1] the topological gradient doesn't depend on the geometry of ω. The formula applies for all shapes in 2D. This phenomena is closely related to the Stokes paradox as been pointed in [1,2,44] and is coherent with the results obtained by several authors in similar problems, for example [6,8,10,28,29].
For simplicity in what follows we will work with an origin-centered inclusion, that means: ω z,ε = ω 0,ε =: ω ε also consider Ω ε := Ω 0,ε . The procedure for all z ∈ Ω is exactly the same just by taking account the change of variable y = z + εx, instead of y = εx that we will use.

Asymptotic expansion of the solution of the Stokes problem
In order to provide an asymptotic expansion of the Kohn-Vogelius functional J KV , we need first an asymptotic expansion of the solution of the Stokes problems (2.4) and (2.5).
Unlike the three-dimensional case, the two-dimensional problem cannot be approximated by an 'exterior problem', which in general in this case doesn't have a solution which vanishes at infinity. This kind of problem has been treated by Bonnaillie-Noël and Dambrine in [11] for the Laplace equation in the plane: we will follow this procedure in order to find a suitable approximation for the Stokes problem.
We recall that we here focus on the detection of a single obstacle (see the beginning of Section 3). This section is devoted to the proof of the following proposition: Proposition 4.1. The respective solutions u ε D ∈ H 1 (Ω z,ε ) and u ε M ∈ H 1 (Ω z,ε ) of Problems (2.4) and (2.5) admit the following asymptotic expansion (with the subscript = D and = N respectively): where (U , P ) ∈ H 1 (Ω) × L 2 0 (Ω) solves the following Stokes problem defined in the whole where E is the fundamental solution of the Stokes equations in R 2 given by (3.2). The We recall that we will detail the proof in the case z = 0 (see Remark 3.3).

Defining the approximation
As we mentioned above, the approximation should be done in a different setting compared to the three-dimensional case, following the same strategy as in [11]. This basically consists in building 'a correction term' to the solution given by E(x − z)u 0 which has a logarithmic term and then tends to infinity at infinity and is not of finite energy in R 2 \ ω. Therefore it has to be considered only in Ω. To this, we consider the pair (U , P ) ∈ H 1 (Ω) × L 2 0 (Ω) solution of Problem (4.1) and we combine these solutions with unknown scale parameters a(ε) and b(ε). Imposing the desired scales to the error function, we will be able to determine the scale factors a(ε) and b(ε) which define completely the approximation for u ε . Here, we will detail the Dirichlet case, the treatment of Neumann case is analog.
Consider the solution (U D , The idea is to combine this solution and the function C D to build a proper corrector. To build this, we search coefficients a(ε) and b(ε), such that the error r ε D defined by: Notice that the remainder r ε D satisfies: where p r ε D is defined in analogous way with pressure terms, that is For x ∈ ∂Ω, we have: Let us assume for a while that ω is a disk. Then, for x ∈ ∂ω ε , there exists X ∈ ∂B(0, 1) such that x = εX and we have We can expand the terms U D (εX) and u 0 D (εX) via Taylor developments: and thus, we get (noticing that O(ε) is contained in o(1)): where i = 1, 2. Therefore, we have the linear system in unknowns (a(ε), b(ε)): We easily get that b(ε) = −a(ε) which implies: This vectorial equality implies two possible choices for a(ε), recalling that where c 1 and c 2 are two positive constants. This leads that a(ε) can be expressed as a(ε) = 1 C−log ε for another positive constant denoted by C in the two possible cases, and then, we get the following scale: It is important to notice, as been pointed in [11,Remark 2.2], that this construction is performed in the case of a disk, where |x| = ε for x ∈ ∂ω ε . In the general case, ω is not a ball and then log |x| = log ε for all x ∈ ∂ω ε and one has to add correctors as performed by Maz'ya et al. in [37, Section 2.4, p. 60-64]. This correction of log ε is of order zero, is then negligible with respect to the logarithmic term. The linear system in (a(ε), b(ε)) remains unchanged and so h ε is still the same rational fraction. Hence, we approximate u ε D by: Analogously, we approximate u ε M by:

An explicit bound of r ε D and r ε M with respect to ε
The Dirichlet case Notice that, in this case, the rest r ε D satisfies: The key point to obtain a bound of r ε D is the following lemma. We postpone its technical proof in Appendix C.
There exists a constant c > 0 (independent of ε) such that: Using this lemma, we get: Notice that: We have used a Taylor development of u 0 D in the last equality and ζ x is some point in the line which joins 0 and εX. Now, recalling that ∇u 0 D is uniformly bounded and using the boundness of U D and the other terms by their definition, we get that: Therefore: which concludes the proof of Proposition 4.1 with = D.
The Neumann case In this case, the rest r ε M satisfies: where the pressure associated to C M is defined explicitly by the expression In order to be able to bound this rest, we use the following lemma. Its technical proof is postponed in Appendix D.
There exists a constant c > 0 (independent of ε) such that: Thanks to this lemma, we know that there exists a constant c > 0 independent of ε, such that: We have In fact, for and, choosing η such that η 1,Ω\B(0,1) = φ 1/2,O , we obtain that The same procedure for U M in Ω ε instead of Ω\B(0, 1) gives the bound for σ(U M , P M ).
Remark 4.4. Notice that we need to consider the set Ω \ B(0, 1) for σ(C M , Π M )n in order to obtain a bound independent of ε: we need to consider a set sufficiently away from zero, due to the definition of C M . For U M , we don't have this problem because it is defined in the whole Ω.
Now we need estimates for the functions U M and C M . Notice that, from the well posedness of the problem (4.1) with = N , we have: For C M we will need a bound for the term |C M | 1,Ω\B(0,1) , for this, first notice that |∇C M | = O(1/r) and let R big enough such that Ω ⊂ B(0, R), therefore: We finally get: Then, from (4.12), The other term of (4.11) is treated identically as in the Dirichlet case (see (4.6)) and therefore, we get: which concludes the proof of Proposition 4.1 with = N .

Proof of Theorem 3.1
We recall that we will detail the proof only for the case of an origin-centered inclusion, i.e. z = 0 (see Remark 3.3).

A preliminary lemma
First we need an estimate of the norm · 1/2,∂ωε of an uniformly bounded function. Here · 1/2,∂ωε has to be seen as the trace norm Lemma 5.1. Let ε ∈ (0, 1/2). If u ∈ H 1 (Ω) is such that its restriction to ω 1 (i.e. ω ε for ε = 1) is C 1 , then there exists a constant c > 0 independent of ε such that Proof. From Theorem A.1, there exists a constant c > 0 independent of ε such that Since u is uniformly bounded on ∂ω ε , we use the change of variables y = εx to prove that there exists a constant c > 0 independent of ε such that Moreover, using the changes of variables x = εX and y = εY , the fact that u(εX) = are some points in the lines which join 0 to εX and εY respectively due to a Taylor expansion), there exists c > 0 independent of ε such that Therefore, we get:

Splitting the variations of the objective
Now, we turn our attention to the Kohn-Vogelius functional given by We first recall the following decomposition:

Asymptotic expansion of A M
We follow here a similar strategy as the one used in the 3D case detailed for example in [20], in contrast to that work we rely on the Stokes fundamental solution properties and the definition of the approximation problem instead of single layer formulas present in the 3D case. We know using elliptic regularity that ∇u 0 M is uniformly bounded on ω ε . Thus We recall that: , C M is given by (4.1 bis) (with = N ) and Π M is given by (4.8). Then the following equality holds Let us first focus on the first term in the right-hand side of (5.3). Using the same argument as the one used in the deduction of (4.12), we get: Therefore, using the explicit upper bound of u 0 M 1/2,∂ωε given by Lemma 5.1, we have Then, using the explicit upper bound of r ε M 1,Ωε given by Proposition 4.1, we obtain For the other term we study each term separately. For this recall that: u 0 We get the last equality because ∇u 0 M is uniformly bounded and: Therefore: Gathering (5.2), (5.5) and (5.6), we obtain

Asymptotic expansion of A D
We recall that and that: where (U D , P D ) ∈ H 1 (Ω) × L 2 0 (Ω) solves (4.1), C D is given by (4.1 bis) (with = D and z = 0) and the pressure associated to C D is defined explicitly by the expression Proceeding as in the previous section 5.3, we prove that Moreover, using Green's formula, we have Now, let us study Proceeding as in the previous section 5.3 (see inequality (5.5)), we use an inequality similar to (5.4), the asymptotic expansion of u ε D given by Proposition 4.1 and Lemma 5.1 to obtain For the other term, we do similar computations as in A M to prove that In order to have a suitable pair (measure g, domain ω * ), we use a synthetic data: we fix a shape ω * (more precisely a finite number of obstacles ω * 1 , . . . , ω * m ), solve the Stokes problem (2.4) in Ω\ω * using another finite elements method (here a P 2-P 1 finite elements discretization) and extract the measurement g by computing σ(u, p)n on O.
In the practical simulations that we present, we add circular objects. In order to determine the radius of these spheres, we use a thresholding method. For an iteration k, it consists in determining the minimum argument P * of the topological gradient δJ KV in Ω\ k j=1 ω j and in defining the set P of the points P ∈ Ω\ k j=1 ω j such that δJ KV (P ) = δJ KV (P * ) + 0.25 * |δJ KV (P * )| .
Then we fix a minimum radius r min := 0.01 and we define the radius of the k th sphere by r k := max r min , min P ∈P (|x P − x P * | , |y P − y P * |) . (6.1) Notice that this method obviously depends on the mesh. We use the classical topological gradient algorithm (see for example [22,29,31,6]) that we recall here for reader's convenience: Algorithm 1. fix an initial shape ω 0 = ∅, a maximum number of iterations M and set i = 1 and k = 0, 2. solve Problems (2.4) and (2.5) in Ω\ k j=0 ω j , 3. compute the topological gradient δJ KV using Formula (5.9), i.e.
We add to this algorithm a stop test (in addition of the maximum number of iterations). In every iteration, we compute the functional J KV . This non-negative functional has to decrease at each iteration. Thus, we stop our implementation when it is not the case, i.e. when J KV Ω\ k+1 j=0 ω j > J KV Ω\ k j=0 ω j . Notice that with this algorithm, we add only one object at each iteration. This method can be slower than the one proposed by Carpio et al. in [17]: they can add several obstacles simultaneously adding points where the topological derivative is large and negative, selecting well calibrated thresholds. The same authors in [16] detailed this approach: they introduce a non-monotone scheme that allows to add and remove points, to create and destroy contours at each stage and even to make holes inside an object. However, in our case, adding only one object at each iteration seems to be more appropriate because otherwise objects can be added wrongly. Moreover, notice that step 5 comes to the assumption that the objects are well separated. Finally, since we assumed that the obstacles are far away from the exterior boundary, we have to take away the added objects on it. Then, if the minimum of the topological gradient is on the exterior boundary, we push the added inclusion inside with a depth 0.005 in the rectangular cases. In origin-centered circular domain we push the added inclusion inside in a quantity proportional to the point, i.e. if the detected point is (x * , y * ) we force it to be (η · x * , η · y * ) where η is usually 0.95 or 0.9.

First simulations
First we want to detect three circles ω * 1 , ω * 2 and ω * 3 centered respectively in (0.475, −0.235), (−0.475, −0.225) and (0.470, 0.150) (i.e. near from the exterior boundary) with shared radius r * = 0.013. The detection is quite efficient (see Figure 2). Indeed we detect three objects with shared radius r = 0.01 for ω * 2 and ω * 3 and r = 0.015 for ω * 1 , we summarized the results in Table 1. Here, we stop the algorithm because of the functional increases as we can see in Figure 3. Notice that some iterations are being made just to adjust the size of a detected object. We can also remark that the values of the cost functional are still relatively high and this refers to the fact that, up to our knowledge, there does not exist a theoretical result of convergence of this algorithm yet.
In this first simulation, the objects are very far away from each other. But what happens when the obstacles are close from each other? Figure 4 shows that the detection of close objects is efficient if the distance between the obstacles is big enough. Indeed, we want to detect three circles ω * 4 , ω * 5 and ω * 6 centered respectively in (−0.475, −0.225), (0.470, 0.100) and (0.470, 0.130) with shared radius r * = 0.01. We obtain just two circles with shared radius r = 0.01 as summarized in Table 2. However if we increase the distance of the near circles enough, considering now, for example, the circle ω * 6bis centered at (0.470, 0.205) we get an efficient detection of the three circles, as we summarize in Figure 5 and Table 3. The distance needed for an efficient 'differentiation' between the objects is relatively high: the required distance in this case is about 2r min . Now the question we asked is: can we detect other shapes than spheres? Thus, we want   The next example deals with a more complex geometry, we have to detect a circle ω * 10 and a non convex object ω * 11 composed by several circle arcs as a boundary. The algorithm is capable to detect both objects and increase the radio of the approximating ball for the non convex object in order to cover it properly. The results are adjoint in Figure 7. In conclusion of these first simulations, this method permits to give us the number of objects we have to determine and their qualitative location if they are separated enough. Moreover, it is efficient to detect different types of shapes, including objects with corners, or even non convex obstacles, in the sense that this topological algorithm is able not only to find the number and relative location of this objects, it is also able to determine their 'relative size' (with respect to its topological set diameter, for example).

Influence of the distance to the location of measurements
As been pointed out in [20] in the 3 dimensional case, the distance to the location of measurements is fundamental in order to get a good detection of the objects. In the following table 4, we notice that, when we move the object away from the boundary of measurements, we get a worse estimate of their location, and in a extreme case a completely wrong detection: more objects than the expected ones. This simple example shows that in our case we have the same problem as in 3-dimensional case: when we try to detect an object which is 'far away' from the boundary, the detection tends to locate it near the boundary (one of the coordinates is correctly estimated) but, as the distance increase, we get some problematic behavior, as we see in our example when the algorithm declares more objects than the real ones. This phenomenon of bad detection can be explained by the regularizing behavior of the Stokes equations (which is related to the behavior of the fundamental solution (3.2)). We emphasize this difficulty of detection pointing out that the functional J KV and its topological gradient are less sensitive to the addition of obstacles when they are far away from the exterior boundary.

Influence of the size of the objects
We now want to study how the size of an object (or several objects) modifies the quality of the detection given by our algorithm. In order to do that, we start by testing how is the detection of a single circle while we increase the radius. Notice that we consider the circle near to the boundary in order to get the best possible approximation as we have seen in the previous section. The following table 5 resumes this first test.  From this we can notice that, when the object is relatively small, the detection is quite efficient, but the quality is decreasing when the object becomes 'too big'. Notice that the main error is linked with the size of the approximation object, and not with their relative position. A more extreme example is putting a 'very big sized' object. In that case, which can be seen in Figure 8, we notice that the detection is completely wrong: we get an incorrect estimate of the number of objects. Interesting results we get when we introduce several objects with higher size, as we can see in Figure 9: the coordinate location is relatively good, but the size approximation tends to stuck in one of the objects. The algorithm choses one of the objects in order to Figure 9: Detection of several 'big sized' objects adjust its size on each iteration.
We can conclude that the detection of the objects depends strongly on the number of them and the size of them. If we are trying to detect a single object we get a reasonable tolerance on the size of the object in order to get a good estimates of their position and size, and only 'big objects' are badly detected. In the case of several objects, the algorithm tends to predict their relative location but only the size of one object is improved between iterations.

Simulations with noisy data
We now want to study how robust is our algorithm in presence of noisy data. For this, we decompose the measurement g = g 1 e 1 + g 2 e 2 (where (e 1 , e 2 ) is the canonical basis of R 2 ) and we consider the following noisy versions of g 1 and g 2 : where u 1 , u 2 are random variables given by an uniform distribution in [0, 1) and σ > 0 is a scaling parameter. Notice that this definition implies that the data g 1 and g 2 are contaminated by some relative error of amplitude σ in L 2 (O). Then, the noisy data will be: From this table we can observe that our algorithm is able to detect with precision the number and relative position of several small obstacles near the boundary O where the measurements are taken, when the boundary data g is contaminated with a moderated amount of noise. When the boundary data contains a higher level of noise, the relative position becomes less precise and finally the algorithm detects an incorrect number of obstacles and therefore the detection becomes completely wrong.

A blending method which combines the topological and geometrical shape optimization algorithms
The previous numerical simulations show that, using the topological gradient algorithm, one can detect the number of objects and their qualitative location but we do not have informations about the shapes of the objects. Hence it can provide initial shapes for an optimization method based on the boundary variation method for which we have to know the number of connected objects we have to reconstruct (see [21]). We present here a combination of these two approaches in order to find the number of objects, their locations and their shapes. As mentioned in the introduction, combinations of several shape optimization methods was recently tested by several authors. The most of them used the level set method (see [3,32,15]). We also mention the algorithm proposed by Pantz et al. in [42] which uses boundary variations, topological derivatives and homogenization methods. We here present an algorithm only based on the classical shape gradient and the topological gradient, without using the level set method or some homogenization methods.
We first recall some theoretical results concerning the computation of the shape derivative of the Kohn-Vogelius functional (see [21] for details). We precise that, in this part, in order to simplify the notation, we will not use the index ε : hence we will use ω = ω ε , u D = u ε D and u M = u ε M .

Shape derivative of the Kohn-Vogelius functional
Let d 0 > 0 fixed (small). We define O d 0 the set of all open subsets ω of Ω with a C 1,1 boundary such that d(x, ∂Ω) > d 0 for all x ∈ ω and such that Ω\ω is connected. The set O d 0 is referred as the set of admissible geometries. We also define Ω d 0 an open set with a C ∞ boundary such that To define the shape derivatives, we will use the velocity method introduced by Murat and Simon in [39]. To this end, we need to introduce the space of admissible deformations For details concerning the differentiation with respect to the domain, we refer to the papers of Simon [45,46] and the books of Henrot and Pierre [33] and of Sokołowski and Zolésio [48]. We consider a domain ω ∈ O d 0 . Then, we have the following proposition (see [21,Proposition 2]): Proposition 7.1 (First order shape derivative of the functional). For V ∈ U , the Kohn-Vogelius cost functional J KV is differentiable at ω in the direction V with where (w, q) is defined by Moreover, Proposition 4 in [21] explains the difficulties encountered to solve numerically this problem. Indeed, the gradient has not a uniform sensitivity with respect to the deformation direction. Hence, since the problem is severely ill-posed, we need some regularization methods to solve it numerically, for example by adding to the functional a penalization in terms of the perimeter (see [14] or [24]). Here we choose to make a parametric regularization using a parametric model of shape variations.

Framework for the numerical simulations
We follow the same strategy than in [21] that we recall for readers convenience. We restrict ourselves to star-shaped domains and use polar coordinates for parametrization: the boundary ∂ω of the object can be then parametrized by where x 0 , y 0 ∈ R and where r is a C 1,1 function, 2π-periodic and without double point. Taking into account of the ill-posedness of the problem, we approximate the polar radius r by its truncated Fourier series for the numerical simulations. Indeed this regularization by projection permits to remove high frequencies generated by cos(kθ) and sin(kθ) for k >> 1, for which the functional is degenerated. Then, the unknown shape is entirely defined by the coefficients (a i , b i ). Hence, for k = 1, . . . , N , the corresponding deformation directions are respectively, This equality is simply that

Algorithm
The first step is the use of the previous topological gradient algorithm described in Section 6.1. It permits to obtain the number of objects and their qualitative location which represents an initial shape ω 0 for a reconstruction using a boundary variation method. Then, the geometrical optimization method used for the numerical simulation is here the classical gradient algorithm with a line search (using the Wolfe conditions: see for example [40, eq. (3.6) page 34]): Algorithm 1. fix a number of iterations M and take the initial shape ω 0 (which can have several connected components) given by the previous topological algorithm, 2. solve problems (2.4) and (2.5) with ω ε = ω i , 3. extract ∇u D , ∇u N , p D and p M on ∂ω i and compute ∇J KV (Ω\ω i ) using formula (7.1), 4. use the Wolfe conditions to compute a satisfying step length α i , 5. move the coefficients associated to the shape: ω i+1 = ω i − α i ∇J KV (ω i ), 6. get back to the step 2. while i < M .
We precise that we here use the adaptive method described in [21,Section 4.3]. It consists in increasing gradually the number of parameters during the algorithm to a fixed final number of parameters. For example, if we want to work with nineteen parameters (which will be the case here), we begin by working with two parameters during five iterations, then with three parameters (we add the radius) during five more iterations, and then we add two search parameters every fifteen iterations. The algorithm is then the same than the one described above only replacing step 5. by where ω i (1 : m) represents the m first coefficients parametrizing the shape ω i (the same notation holds for ∇J KV (Ω \ ω i ) (1 : m)). The number m grows to the fixed final number of parameters following the procedure described previously.
To finish, we precise that we use the finite elements library Mélina (see [35]) to make this geometrical shape optimization part.

Numerical simulations
The framework is the following: we assume the kinematic viscosity ν is equal to 1, the exterior boundary is assumed to be the unit circle centered at the origin and we consider the exterior Dirichlet boundary condition where n = (n 1 , n 2 ) is the exterior unit normal. Notice that f is such that the compatibility condition (2.1) is satisfied. We assume that we make the measurement on the whole disk ∂Ω except the lower right quadrant. Here, we want to detect two squares ω * 12 and ω * 13 centered respectively at (−0.6 , 0.3) and (0.6 , 0.3) with a distance between the center and the vertices equal to 0.2. The first step, which is the topological approach, leads to two circles of radius 0.15 centered respectively at (−0.573 , 0.328) and (0.533 , 0.328) (see the 'initial shape' in Figure 10). Since the real objects are "big", we impose here r min = 0.15 in the topological algorithm (see (6.1)). This means that, practically, we assume that we know the characteristic size of the objects, i.e. if the objects are small or big. Then, the shape optimization algorithm leads to a good approximation of the shapes, at least for one of the obstacle (see Figure 10). We also underline the fact that, after the topological step, the cost of Figure 10: Detection of ω * 12 and ω * 13 with the combined approach (the initial shape is the one obtained after the "topological step") -Zoom on the improvement with the geometrical step for ω * 13 the functional is here about 1.26 and that, after the geometrical step, we obtain a cost about 2 · 10 −2 which qualitatively means that we improved the detection.
In conclusion, this blending method which combines the topological and the geometrical shape methods leads to good result in the identification of obstacle immersed in a fluid: we detect both the number of obstacles, their locations and their shapes.

Conclusion
Using a Kohn-Vogelius approach, we have detected the number of potential objects immersed in a two dimensional fluid and their qualitative location. To do this, we have computed the topological gradient of the considered Kohn-Vogelius functional using an asymptotic expansion of the solution of the Stokes equations in the whole domain when we add small obstacles inside: we adapted the usual 3D techniques to the two dimensional setting case, in which the classical asymptotic expansion of the solution of the Stokes equations are no longer valid. We obtain a formula valid for any geometry of small obstacles, which is a particular characteristic of the two dimensional setting of the problem. We have made some numerical attempts which have shown that 'not too big' obstacles close to the part of the boundary where we make the measurements can be detected. Once these restrictions are satisfied, the detection is quite efficient, even for objects with corners or non convex shapes. Finally, we have proposed and implemented an algorithm which combines the topological sensitive analysis approach with the classical shape derivative approach. This blending method led us to detect the number of objects using a topological step and, if this first step actually gives the total number of obstacles, a geometrical shape optimization step detects their approximate location and approximate shape from only the boundary measurements. This method gives interesting results in the simulations.

A A result concerning the space of traces
Here we recall a result used in the paper concerning the boundary values of functions, in particular when domains depend on a parameter (see [38,Chapter 4]): Let Ω and ω be two bounded simply connected domains of R N (N ≥ 2) of class C 0,1 . Let p ∈ (1, +∞), ε ∈ (0, 1/2) and ω ε := εω. Let us assume that ω ε ⊂ Ω and that there exists a constant c > 0 depending only of N , p, ω and Ω such that d(ω ε , ∂Ω) > cε. Then for p > N, Consider now the space . It is standard to check that The dual space of where p is such that Notice that these spaces are reflexive Banach spaces with respect to the norms:

B.2 The exterior Stokes problem in two dimensions
The following results are presented in [44], we present them here for reader's convenience. We first recall the following lemma concerning the Stokes problem in the whole space R 2 : Then every solution which is a tempered distribution should be a polynomial.
Proof. Applying Fourier transform to (B.1) we immediately notice that the support ofû andp is contained in {0}. Therefore, those distributions should be a finite sum of Dirac deltas, which implies that u and p are polynomials.
In such case, we have: Now, let us define: v := E * T , q := P * T , where (E, P ) is the fundamental solution of the Stokes system given by (3.2) and * denotes the convolution product. Then, Now notice that the pair (u − v, p − q) solves (B.1), then by the previous lemma this solution should be a polynomial. Then: where U 1 and P 1 are polynomials and t = D(u)n.
Using a Taylor development for u we get a logarithmical term, due to: where θ(y, x) = y − αx with α ∈ (0, 1), then: But log ∈ W 1,2 0 (ω c ), which implies: Also, due to U 1 ∈ W 1,2 0 (ω c ), we must have U 1 = λ, where λ is a constant. Therefore, we have: A similar reasoning gives p(y) = O(1/r), where r = y , and P 1 = 0. Therefore we have: The study of the function W in (B.3) will be useful for important results: we study a priori estimates for this function, in a similar way as Guillaume in [28].
Using these estimates we can bound the following quantities: From the previous inequalities we can estimate the size of the function W x ε in Ω ε , indeed, by change of variables (recall the equalities given by (B.6), (B.7)), we get that, for small ε: C Proof of Lemma 4.2 The following results are based in the ones presented in [44,Chapter 3]. We will use the notations introduced in section B.2.2. The statements will be presented for general domain Ω z,ε , but for sake of simplicity we will work in the proofs in the case z = 0, the general case comes from a linear change of coordinates. The proof of Lemma 4.2 is decomposed in the following three lemmas.
Combining the last two inequalities we get the desired result.
Proof of Lemma 4.2. If ϕ is constant on ∂ω ε and Φ = 0, the previous lemma gives the desired result. If Φ = 0 another previous lemma gives the desired estimate. So, let's focus on the case where ϕ is not constant. Let V the bounded solution of on ∂ω.
Using the previous lemmas we can bound the terms of this function, and λ can be bounded thanks to (B.5). Finally we have that W x ε satisfies the desired estimate by Lemma B.3 and we conclude by triangle inequality.