On the detection of small moving disks in a fluid

. We are interested in determining the positions and the velocities of moving rigid solids in a bounded cavity ﬁlled with a perfect ﬂuid. We assume that the solids are small disks and that they move slowly. Using an integral formulation, we ﬁrst derive the asymptotic expansion of the DtN map of the problem as the diameters of the disks tend to zero. Then, combining a suitable choice of exponential type data and the DORT method (which is usually used in inverse scattering for the detection of point-like scatterers), we propose a reconstruction method for the unknown positions and velocities.

We also introduce the unit normal n to ∂F directed towards the exterior of the fluid and τ the unit tangent vector to ∂F such that τ = n ⊥ (throughout this paper, we set x ⊥ := (−x 2 , x 1 ) for all x = (x 1 , x 2 ) ∈ R 2 ). We assume that every solid moves with a velocity V m ∈ R 2 . At every given time, we can consider the Eulerian velocity field U(x), x ∈ F, of the fluid. We assume that the flow is irrotational, We shall address the inverse problem of determining the positions and the velocities of the solids by using actuators and sensors located on the outer boundary Γ. We assume that we can prescribe the normal velocity of the fluid on Γ, leading to the modification of (1.2b) into where F is given. Regarding the available output, we assume that we can measure the tangential velocity U · τ on Γ. According to conditions (1.1a) and (1.1b), we classically introduce the stream function ψ : F → R such that U = −∇ ⊥ ψ in F . Then, (1.1a) and ( With these settings, the measurement reads U · τ = −∂ n ψ. The detection problem under consideration is then to recover the centers, the diameters, and the velocities of the solids from the DtN map Λ : f ∈ H 1/2 (Γ) −→ ∂ n ψ ∈ H −1/2 (Γ).
Remark 1. Instead of using the stream function, one can equivalently formulate the problem in terms of the potential function and obtain Neumann boundary conditions instead of Dirichlet in (1.4).
Although such geometric inverse problems of detecting moving solids in a fluid appear in many applications, the associated literature is quite limited, as most contributions deal with the case of motionless solids (i.e., obstacles). In particular, the case of a fluid described by the Stokes equations has been investigated using optimization methods by Caubet et al. in [12,5,11,10] and more recently by Bourgeois and Dardé in [6] using the quasi-reversibility method combined to a level set method. The case of small obstacles has been studied in Caubet and Dambrine [9]. Regarding moving obstacles, Conca et al. show in [15,14] that the position and the velocity for a single disk moving in a perfect fluid can be recovered from one measurement of the velocity on part of the boundary. Linear stability estimates are also provided. In Conca, Malik, and Munnier [16], the authors consider a moving rigid solid immersed in a potential fluid and provide examples of detectable (ellipses for instance) and undetectable shapes. Conca, Schwindt, and Takahashi obtained in [17] an identifiability result in the case of a rigid solid immersed in a viscous fluid. In this work, we restrict the analysis to the case where the solids are small disks (see Figure 1) of typical size ε, so that rotation plays no role. For such configurations, we show an identifiability result and we provide a reconstruction method. More precisely, we denote each closed disk D ε m , m = 1, . . . , M, and we denote by γ ε m its boundary. We assume that D ε m is centered at r m and is of radius εR m , where the parameter ε is meant to tend to 0. Denoting by F ε = Ω \ ∪ M m=1 D ε m the domain occupied by the fluid and letting ψ ε be the corresponding solution of (1.4), our inverse problem can be formulated as follows: Knowing the DtN map is it possible, and if so how, to recover the number M of disks, their positions r m , m = 1, . . . , M, their (rescaled) radii R m , and their velocities V m ?
We answer this question in several steps. In section 2, we first derive the asymptotic expansion of the DtN operator Λ ε as ε → 0 + (Theorem 2.1) using a boundary integral formulation of the forward problem. In section 3, we combine this expansion with the DORT method 1 to recover the number of disks and their positions, provided they are distant enough. Initially introduced by Fink and Prada [18], this method has been justified mathematically and used in the framework of wave systems for the detection of distant point-like scatterers in acoustics [21,26,7] and electromagnetics [4]. Once the positions have been determined, the velocities and rescaled radii can be easily recovered using suitable data. Finally, we collect in section 4 some numerical examples to illustrate the efficiency of the proposed reconstruction method.

Asymptotic expansion for small disks.
Although the literature dealing with small inhomogeneities or small inclusions is quite rich (see, for instance, the review paper [2] and the book [3] by Ammari and Kang), the needed asymptotics is not, as far as we know, directly available for our problem (1.4), nor for the equivalent Neumann conjugate problem satisfied by the velocity potential. Quoting only the tightly related works, let us mention Il in [22] who studied an elliptic boundary value problem set in a three dimensional (3D) domain containing a small hole, Maz'ya, Nazarov, and Plamenevskij who studied the Laplace problem with only one small hole (see [24, p. 59] for the Dirichlet case and [24, p. 291] for the Neumann case) and Friedman and Vogelius [19,Lemma 3.3] who obtained the first order asymptotics for several infinitely conducting small inhomogeneities (i.e., constant Dirichlet condition on the boundaries of the small holes and Neumann condition on the exterior boundary).
The rest of this section is devoted to the proof of this result. The first ingredient is the following reciprocity identity, which can be easily obtained using Green's formula: where we have set This formula shows in particular that the asymptotics of the bilinear form associated with Λ ε can be obtained from the asymptotics of ψ ε (which is given) and ∂ n ψ ε on γ ε . In order to obtain these asymptotics, we use a superposition result, usually referred to as Kirchhoff principle in the context of fluid dynamics. Thus, we have and v ε m and w ε m solve for every m = 1, . . . , M the following boundary value problems: where the constants a ε m , b ε m,k , and c ε m,k are such that where these integrals should be understood in the sense of duality We detail in section 2.1 the derivation of the asymptotics of u ε , based on a boundary integral formulation of problem (2.4)-(2.5). Then, as the proofs are quite similar, we only sketch in section 2.2 the proofs for v ε m and w ε m . These three asymptotic expansions are used to obtain the asymptotics of , and the result follows then from (2.2). For the sake of clarity and for reader's convenience, we preferred to give a constructive proof of these expansions.

Asymptotic expansion for the function u ε .
The next result collects well-known properties of the single layer potential that we use in what follows (see, e.g., McLean [25], Rjasanow and Steinbach [27], or Steinbach [28,Chapter 6]).
log |x| denote the Green function of the operator −Δ in R 2 . Let γ denote the smooth boundary of a bounded (possibly multiply connected) domain of R 2 . We introduce the single layer potential associated to a given density p: Then, the following assertions hold true: The trace of S γ p on γ is given by the boundary integral operator 3. The normal derivative of S γ p on γ is given by the so-called jump condition in which n denotes the exterior unit normal to γ (pointing from the interior bounded domain inside γ to the unbounded exterior) and the + and − signs refer, respectively, to the normal derivatives coming from the exterior or the interior of γ. 4. Assume that γ has only one connected component. Then, there exists a unique density ψ eq ∈ H −1/2 (γ), called the equilibrium density of γ, such that S γ ψ eq is constant on γ and satisfies the normalization condition γ ψ eq (y) dσ y = 1. We seek u ε as a single layer potential: where p m ε and q ε are densities in H −1/2 (Γ) and H −1/2 (γ ε m ) to be determined.
Taking into account Remark 2, system (2.4)-(2.5) satisfied by u ε is equivalent to the following system of (M + 1) coupled integral equations on γ ε 1 , . . . , γ ε M and Γ: Defining the curve γ = r + ε −1 (γ ε − r ) and setting r m = r − r m , a simple rescaling leads us to rewrite (2.11a)-(2.11b) in the equivalent form: The circulation-free condition (2.11c), once the moving disks are rescaled, becomes For every ε > 0, we introduce the following a priori decomposition of the densities: Likewise, we write the constants a ε as where the terms appearing in the right-hand sides of (2.14)-(2.15) have to be determined. Classically, this can be achieved by plugging the above a priori expansions in (2.12) and identifying the terms of same order in the expansion. More precisely, (2.12) read while the circulation free condition (2.13) implies the following conditions: For every X, Y ∈ R 2 , X = 0, we have Using the above Green's function asymptotics in (2.16) and identifying the terms of order 0 in ε yields In particular, recalling the notation U f introduced in (2.1), we have S Γ q 0 = U f in Ω. Choosing p ,0 = 0 and a ,0 = U f (r ) for every = 1, . . . , M, the first equation in (2.18) is then fulfilled. Now, we identify the terms of order 1 in ε of system (2.12), and get From the last equation, we have q 1 = 0. Taking the scalar product of the first equation of (2.19) with the equilibrium density ψ eq of γ (see Proposition 2.2), we get that ∇U f (r ) · γ (x − r )ψ eq (x) dσ x = a ,1 . Since the moment of order 1 of the equilibrium density of a circle coincides with its center (this follows in particular from relation (B.4) in Lemma B.1 of Appendix B, with V = 0 and c = (2πR ) −1 ), we have a ,1 = 0. The first equation of (2.19) then yields S γ p ,1 (x) = ∇U f (r ) · (r − x), and according to (B.4) in Lemma B.1, we thus have p , . Let us now consider the terms of order 2 in (2.12b) (this will be enough to derive the expected asymptotic expansion of Λ ε ). We easily obtain that Applying once again Lemma B.1, with V = −∇U f (r ) and c = ∇U f (r ) · r , it follows easily from (B.4) that for every m ∈ {1, . . . , M}, Consequently, (2.20) shows that q 2 solves the boundary integral equation Let us focus now on the remainders P ε ∈ H −1/2 (γ ), Q ε ∈ H −1/2 (Γ), and A ε ∈ R appearing in (2.14) and (2.15), and prove that they are bounded. This will provide a justification of the formal a priori expansion of p ε , q ε and a ε . To do so, we need the following general result, proved in Appendix A.
Proof. On Γ, we have Multiplying by Q ε , we get (recall that Q ε 2 Lemma 2.3 implies the existence of a constant C > 0 independent of ε > 0 such that On γ , we have Multiplying by the equilibrium density ψ eq of γ , integrating on γ , and applying Lemma 2.3, we get, since γ P ε (y) dσ y = 0, where the constant C > 0 does not depend on ε. Multiplying now (2.24) by P ε , integrating on γ , and using once again Lemma 2.3, we get Combining the estimates (2.23) and (2.26), we immediately obtain that which clearly implies that P ε Summing up, we have proved the following result. Proposition 2.5. With the above notation, we have for every = 1, . . . , M, εp ε (r + ε(y − r )) = εp ,1 (y) + ε 2 P ε (y), y∈ γ , (2.27a) q ε (y) = q 0 (y) + ε 2 q 2 (y) + ε 3 Q ε (y), y ∈ Γ, (2.27b) , and A ε are uniformly bounded in ε. Now, we can compute the contribution [u ε , U g ] γ ε to the expansion of Λ ε (see (2.2) and (2.6)). On γ ε , we have on the one hand u ε = a ε and, on the other hand, p ε = [∂ n u ε ] = ∂ n u ε (the first equality follows from the jump relation (2.9) and the second one from the fact that the continuous harmonic extension of u ε inside D ε is the constant function a ε ). Using the asymptotic expansions of Proposition 2.5, simple computations lead to the following result.

Asymptotic expansion for the functions v ε m and w ε m .
Let us deal with the function v ε 1 , as the other cases (v ε m , m 2 and w ε m , m = 1, . . . , M) can be treated in a similar way. For the sake of simplicity, we drop the subscript referring to this disk's number in the proof and we merely denote this function by v ε . Similarly, we rename D ε the disk D ε 1 , and γ ε its boundary, while the other disks are renamed D ε m , m = 1, . . . , N := M − 1, with boundaries γ ε m . With this notation, v ε solves where the constants b ε and b ε m , m = 1, . . . , N, are such that We seek v ε as a single layer potential: , and q ε ∈ H −1/2 (Γ) are to be determined. Problem (2.28) is equivalent to the system of integral equations: with the constraints (see Remark 2) γ ε p ε (y) dσ y = γ ε m p ε m (y) dσ y = 0 for all m = 1, . . . , N. Rescaling the above equations, we obtain with obvious notation, with the constraints For all = 1, . . . , N, we introduce the formal asymptotic expansions of the quantities: Using these expressions in system (2.29), the order 0 (in ε) reads which, thanks to (2.29d), leads to q 0 = p 0 = p ,0 = 0 and b 0 = b ,0 = 0. Identifying the terms of order 1 in ε of system (2.29) leads to Multiplying the first equation by the ψ eq γ of γ and integrating over γ, we find that b 1 = − γ ψ eq γ (x)x 1 dσ x = −r · e 1 , since the first momentum of the equilibrium density of a circle is its center. From the other equations we deduce that q 1 = 0, p ,1 = 0, and b ,1 = 0. Finally, S γ p 1 (x) = (x − r) · e 1 for x ∈ γ, and according to Lemma B.1, we have p 1 (x) = (2/R)(x − r) · e 1 . Finally, the identification of the second order terms in ε in (2.29c) yields Identity (B.4) of Lemma B.1 shows that γ (r−y)p 1 (y) dσ y = −2πR 2 e 1 , and therefore, q 2 is the unique solution of the integral equation due to Calderón [8], we use a suitably chosen family of exponential type excitations f and test functions g. Next, we make use of the so-called DORT method to recover the unknown positions r m . Based on time-reversal techniques, this method has been successively used for the detection of well separated point like scatterers from far field measurements, for many wave type systems in acoustics [21,26,7] and electromagnetics [4]. The underlying idea of this method is that the eigenfunctions of some finite rank integral operator (the so-called time-reversal operator, that can be computed from the measurements) generate waves that selectively focus on each target. For the elliptic Laplace problem considered here, let us emphasize that time reversal is not performed experimentally, but only used numerically in the reconstruction algorithm through the use of suitably chosen excitations and test functions. Of course, one can also use other reconstruction methods which are classically used for the detection of point-like targets in scattering theory, like the MUSIC algorithm [13,23,20].
One can easily check that (see, for instance, [1, p. 360 where r min := min m =n |r m − r n | is the minimal distance between the solids. In other words, for every m = 1, . . . , M, ϕ m,0 , ϕ m,1 , and ϕ m,2 constitute approximate eigenfunctions of A as k → +∞ (for high frequency and distant disks), the corresponding eigenvalues being, respectively, 2πR 2 m (simple eigenvalue) and −πR 2 m (with multiplicity two). These approximate eigenfunctions can be used to recover the positions of the disks by constructing the corresponding Herglotz wave we see that Hence, u m, generates a wave that selectively focuses on the target m as k → ∞ and allows us to recover its location. In particular, for = 0, we have and thus, u m,0 reaches its maximum exactly at the point x = r m . These functions, which are well defined since P m (r m ) = 0, enjoy the properties

Numerical tests.
In this section, we present some numerical results to illustrate our reconstruction method. The solutions of the boundary value problems involved, in particular for the data generation, are computed using a MATLAB boundary integral equation solver. 2 As a typical configuration, we consider three small disks located in a (smoothened) rectangular domain Ω. The three small disks are of radii (m + 1)ε, where ε = 10 −3 and m = 1, 2, 3. They are located, respectively, at the points (0.7, −0.1), (−0.7, 0.3), and (0, −0.3) and their velocities are (0, 1), (1, 1), and (1, 0). As a preliminary step, we first generate numerically the data that will be used for the reconstruction. Given a wavenumber k and a uniform discretization of [0, 2π] with mesh size h = 1/(N + 1), we compute the matrix A = (A i,j ) 1 i,j N corresponding to the discretization of the kernel A(α, β). More precisely, we have A i,j = −1/k 2 (Λ ε f i,j , g i,j ), where f i,j , g i,j are the functions obtained in (3.3) for η = k(α i − β j )/2, where α i and β j belong to the chosen discretization of [0, 2π]. This requires solving the forward problem (1.4) for N 2 different right-hand sides f i,j , 1 i, j N . The data are generated using N = 40 angles and then perturbed them artificially by adding 10% of noise.
Our DORT based reconstruction procedure is applied to these noisy data. We compute the eigenvalues and eigenfunctions of the discretized integral operator A h associated to the kernel A. As expected from the theoretical analysis of section 3, it turns out that this integral operator has indeed three significant eigenvalues for a wide range of values of k.
We show on Figure 2 the dependence of the positive eigenvalues with respect to the wavelength λ = 2π/k of the oscillating source (without noise). According to this figure, we can recover the number of disks as soon as k is chosen in the shaded region on Figure 2, roughly corresponding to a constant number of significant eigenvalues. The wavenumber is now chosen to be k = 10, which yields accurate location identification. The centers of the disks can be obtained as the points where the Herglotz waves Finally, let us emphasize that although our reconstruction method is theoretically justified only for small disks (ε → 0), it turns out to still be efficient numerically for "extended" disks as long as the centers are concerned.