RECONSTRUCTION OF OBSTACLES IMMERSED IN AN INCOMPRESSIBLE FLUID

. We consider the reconstruction of obstacles inside a bounded do- main ﬁlled with an incompressible ﬂuid. Our method relies on special complex geometrical optics solutions for the stationary Stokes equation with a variable viscosity.


Introduction
Let Ω ⊂ R 3 be an open bounded domain with smooth boundary. Assume that D is a subset of Ω such that D ⊂ Ω and Ω \ D is connected. Let u = (u 1 , u 2 , u 3 ) be a vector-valued function satisfying where n is the unit outer normal to ∂Ω. On the other hand, the condition u| ∂D = 0 is the no-slip condition. If ν is a positive constant this is the classical stationary Stokes system describing the motion of a homogeneous, incompressible, viscous fluid in the domain Ω \ D. Here we will also treat the more general case ν(x) ∈ C 4 (Ω), 0 < δ ≤ ν in Ω for some constant δ > 0. A physical description of (1) with a variable 2000 Mathematics Subject Classification: Primary: 35R30, 35Q30; Secondary: 76D07. Key words and phrases: Stokes system, complex spherical waves, inverse problem. The first author was financially supported by the Deutsche Forschungsgemeinschaft DFG. The second author is partially supported by NSF. The third author was supported in part by the National Science Council of Taiwan. viscosity is given in [8]. The domain D is considered as an obstacle immersed in an incompressible fluid which fills Ω \ D. For a variable ν(x), the well-posedness of (1) is proved in [8]. More precisely, for any f ∈ H 1/2 (∂Ω) satisfying (2), there exist u ∈ H 1 (Ω) and p ∈ L 2 (Ω) (unique up to a constant) satisfying (1). Denote σ(u, p) = νS(∇u) − p. Then we have σ(u, p)n| ∂Ω ∈ H −1/2 (∂Ω), which is called the Cauchy forces. We define Λ D (f ) = {σ(u, p)n| ∂Ω : (u, p) is a solution of (1) with u| ∂Ω = f }.
Note that Λ D maps H 1/2 (∂Ω) to a "set" of H −1/2 (∂Ω). Nonetheless, any two elements of Λ D (f ) for a fixed f differ by a constant multiple of n. In this problem, we are interested in the inverse problem of determining D from the boundary data (f, Λ D (f )) ∈ H 1/2 (∂Ω) × H −1/2 (∂Ω).
Our study is inspired by Alvarez, Conca, Fritz, Kavian and Ortega's paper [2] in which they considered the object determination problem for the Stokes as well as for the Navier-Stokes equations (with constant ν). They were able to prove uniqueness of this inverse problem for the stationary and the non-stationary case. Part of the work is also motivated by a recent paper of the last two authors [20] where we considered the reconstruction problem of the inclusion for the isotropic elasticity using the so-called complex spherical waves, which was first defined in [9] for the Schrödinger equation.
On the other hand, for the 2D-Oseen equation Kress and Meyer gave in [12] a solution to the unique identification problem for an obstacle if the velocity of the fluid is measured on some arc outside the obstacle. They also study a Newton iteration scheme in order to approximate the solution. The object determination problem in the conductivity equation is much more well studied. We will not try to give a full account of these developments here. For detailed references, we refer to [9].
In this paper we focus on the reconstruction method for the proposed inverse problem. The method we will employ relies on the use of complex geometrical optics solutions. Such kind of solutions were first used in Sylvester and Uhlmann [17] for the inverse conductivity problem. Since then their method became the base for treating many inverse problems. For the elasticity equation such solutions were first used by Nakamura and Uhlmann in [14]. Recently, complex geometrical optics solutions with more general phase functions were constructed in [13] for the Schrödinger equation and in [7] for the Schrödinger equation with magnetic potential. The method used in [13] and [7] relies on Carleman type estimates, which is a more flexible tool in treating lower order perturbations. In this paper we apply the ideas of [7], [13], and [20] in order to construct complex geometrical optics solutions for (1). The real part of the phase function here is a radial function. With the help of these special solutions we are able to detect unknown objects inside the fluid with known background viscosity.
To construct complex geometrical optics solutions for the free system associated with (1) (see (11)), the divergence-free restriction causes serious difficulties. A key observation is that (11) can be transformed to a system with the Laplacian as the leading part. This observation is motivated by the isotropic elasticity system and was also used in [8]. Another difficulty caused by the divergence-free condition is the use of local measurements as one did for the conductivity equation [9] and the isotropic elasticity [20]. The use of local measurements is important in real applications. The divergence-free condition puts a global restriction on the Dirichlet data f of (1), i.e., the compatibility condition (2). This condition does not have a local character. Even though we are not able to localize the measurements exactly, we can still do this approximately in the sense that the measurements outside of a fixed region can be made as small as we wish. This fact is due to the decaying property of complex geometrical optics solutions.

Preliminaries
In this section we will derive important integral inequalities that play a key role in our inverse problem. The derivation here is motivated by Ikehata's work [10]. It turns out that the integral inequalities do not depend on the pressure p. We shall take advantage of this property in the construction of complex geometrical optics solutions. For the obstacle D, in addition to the conditions described above, we assume ∂D ∈ C 2 . Let (u, p) be a solution of (1), then using Green's formula we obtain Now we let (u 0 , p 0 ) be a solution of (4) Let us define (ũ,p) = (u, p) for x ∈ Ω \ D and (ũ,p) = (0, 0) for x ∈ D. Since u = 0 on ∂D,ũ ∈ H 1 (Ω). So we can see that Next, subtracting (5) from (3) leads to On the other hand, we can deduce that Combining (6) and (7) yields The regularity theorem implies for some constant C > 0. In view of Korn's inequality, putting together (8) and (9) gives Now we would like to construct special solutions (u 0 , p 0 ) satisfying (11) div(νS(∇u 0 )) − ∇p 0 = 0 in Ω, div u 0 = 0 in Ω.

Construction of complex spherical waves
Our aim here is to construct special solutions u 0 , called complex spherical waves, from (13). The method here is completely analogue to the elasticity system treated in [20], which is based on the Carleman method introduced in [7] and [13]. We set the matrix operator P h = −h 2 P . More precisely, we have Later on we shall denote the matrix operator The construction here is simpler than the one given in [14] and [16] where the technique of intertwining operators were first introduced. Furthermore, we do not need to work with C ∞ coefficients here. Our first goal is to derive a Carleman estimate with semiclassical H −2 norm for P h . The conjugation of P h with e ϕ/h is given by We first consider the leading operator (hD + i∇ϕ) 2 and denote where A = (hD) 2 − (∇ϕ) 2 and B = ∇ϕ • hD + hD • ∇ϕ. The Weyl symbols of A and B are given as respectively. Let Ω 0 be an open bounded domain such that Ω ⊂ Ω 0 . Accordingly, we extend ν to Ω 0 by preserving its smoothness. We now let ϕ have nonvanishing gradient in Ω 0 and be a limiting Carleman weight in Ω 0 : i.e., ∇ 2 ϕ, ∇ϕ ⊗ ∇ϕ + ξ ⊗ ξ = 0 when ξ 2 = (∇ϕ) 2 and ∇ϕ · ξ = 0. Let us denote ϕ ε = ϕ + hϕ 2 /(2ε), where ε > 0 will be chosen later. Also, we denote by a ε and b ε the corresponding symbols as ϕ is replaced by ϕ ε . Then one can easily check that Arguing as in [13], we get where β(x, ξ) is linear in ξ. Therefore, at the operator level, we have (14) i where β w denotes the Weyl quantization of β.
With the help of (14), we can now estimate Here and below, we define the norm · and the inner product · | · in terms of L 2 (Ω). Integrating by parts, we conclude On the other hand, we observe that and the obvious estimate Using (14), (15), (16), and (17) gives Thus, taking h and ε (h ≪ ε) sufficiently small, we arrive at
The equation (26) is a system of Cauchy-Riemann type. Using the results in [4], [6], or [15], one can find an invertible 4 × 4 matrix G(x) ∈ C 2 (Ω) satisfying (26). So L can be chosen from columns of G. Then, R is required to satisfy Note that e −iψ/h P h L h 2 . Thus Theorem 3.1 implies that which leads to So if we write L = ℓ d and R = r s with ℓ, r ∈ C 3 , then (30) w = e −(ϕ+iψ)/h (ℓ + r) and g = e −(ϕ+iψ)/h (d + s) where r and s satisfy the estimate (29). Therefore, u 0 = ν −1/2 w + ν −1 ∇g − g∇ν −1 is the needed complex spherical wave, i.e., the real part of the phase function of u 0 is a radial function.

Remark 3.2.
Even though the four-vector ℓ d is nonzero in Ω, we cannot conclude that d never vanishes in Ω. However, for any point y ∈ Ω, it is easy to show that there exists a small ball B δ (y) of y with B δ (y) ⊂ Ω such that d does not vanish in B δ (y). We will use this fact in studying our inverse problem in the next section.

Reconstruction of the Obstacle
In this section we shall use complex spherical waves constructed in the previous section and the integral inequality (10) to determine some part of D by boundary measurements (u| ∂Ω , σ(u, p)n| ∂Ω ). To this end, we set u 0 h,t = e log t/h u 0 , with parameters t > 0 and h > 0. We choose the function u 0 to be of the form u 0 = ν −1/2 w + ν −1 ∇g −g∇ν −1 where w = e −(ϕ+iψ)/h (ℓ+r) and g = e −(ϕ+iψ)/h (d+s) are the functions given in (30). With u 0 h,t , we can define p 0 h,t = e log t/h p 0 with p 0 given by (12). Then (u 0 h,t , p 0 h,t ) satisfies (11) for all h, t and since ν is known, we can determine is the observable data. Therefore, we can determine by boundary measurements. We now use the behavior of E t (h) as h → 0 at different t's to reconstruct some part of ∂D. Our method shares the same spirit as Ikehata's enclosure method (see Ikehata's survey article [11]). Furthermore, since we probe the region by complex spherical waves, it is possible to recover some concave parts of inclusions. For other related results, we refer to [20] and references therein.
Then there exists h 0 > 0 such that for t > 0 and h ≤ h 0 the following assertions are true: provided the function d does not vanish in a neighborhood of y.
Proof. The key to obtain estimates of the terms E t (h) is to plug in the special solutions u 0 h,t and to look at the leading terms in view of small h. This leading terms come from the expression S(∇u 0 h,t ). In the case dist(D, x 0 ) > t these are the terms which behave like h −4 and are given by So from the second inequality of (10) we get for h ≤ h 0 and therefore the first assertion is proved.
In order to treat the second assertion we choose a small ball B ǫ ⋐ B t (x 0 )∩D and the special solution u 0 h,t such that the corresponding function d in the definition of u 0 h,t never vanishes in B ǫ . Then the boundary value of u h,t on ∂Ω is given by Since for h ≤ h 0 small enough and j = 1, 2, 3 the inequality and therefore assertion (b).
Analyzing the behavior of E t (h) when y ∈ ∂D ∩ ∂B t (x 0 ) we choose u 0 h,t , as before, such that d does not vanish in the ball B ǫ (y). Pick a small cone with vertex at y, say Γ, so that there exists an η > 0 satisfying Thus, we get from the first inequality of (10) that On the other hand, we can choose a cone Γ with vertex at Combining (31) and (34) yields the assertion (c).

Remark 4.2.
In the proofs of (b) and (c) we need to choose d which is nonvanishing in small subdomains of Ω. Since d depends only on the known background medium, they can be chosen to be nonvanishing near any point in Ω at our disposal. If fact, it suffices to take d which is nonvanishing near the probe front {|x − x 0 | = t}. Different choices of d will give rise to different Dirichlet data f h,t and therefore different measurements. In real applications, we believe that the concerns in (b) and (c) can be ignored.
Next we would like to discuss the possibility of localizing boundary measurements. It turns out that we cannot produce a similar result as in [9] and [20]. This is due to the compatibility condition for the Dirichlet data (2) or, equivalently, the divergence-free condition. Those restrictions are global ones. They cannot be localized. Even though we cannot exactly localize the measurements as we wish, we can still do it approximately. We now describe the method. Let where δ > 0 is sufficiently small. We are going to use the measurements f δ, With such constant c δ,h,t , the boundary value f δ,h,t satisfies the compatibility condition (2). We can see that the constant c δ,h,t tends to zero as h → 0. Indeed, we can estimate in Ω \ B t+δ/2 (x 0 ). We argue that the use of f δ,h,t is better than f h,t because f δ,h,t can be kept constant on most part of ∂Ω and this constant is as small as we want. Let us now define Proof. The main idea is to prove that the error caused by the remaining part of the measurement g δ,h,t := (1 − φ δ,t )f h,t | ∂Ω + c δ,h,t n is as small as any given polynomial order. Let (w δ,h,t , p δ,h,t ) be a solution of (4) with boundary value g δ,h,t . We now want to compare w δ,h,t with (1 − φ δ,t )u 0 h,t . Note that a pressure function associated with u 0 h,t is given by p 0 h,t = e log t/h (div(ν 1/2 w) + 2∆g). For simplicity, we denote v : We can check that and therefore F L 2 (Ω) ≤ Ca Also, we have The system (39) is a nonhomogeneous Stokes system, which is an elliptic system in the sense of Agmon, Douglis, and Nirenberg [1]. Here we are concerned with the regularity estimate of (39): In order to derive (40) for some 0 < a 2 < 1. Using the second inequality of (10) for (Λ D −Λ 0 )(ḡ δ,h,t ), g δ,h,t with u 0 being replaced by w δ,h,t , we get from (41) and the decaying property of u 0 h,t that (Λ D − Λ 0 )(ḡ δ,h,t ), g δ,h,t ≤ Ca 1/h 3 for some 0 < a 3 < 1.
Next we consider (b) and the first inequality of (c) in Theorem 4.1 for E δ,t (h). Choosing ζ = 1 in (42) we get that Thus, combining (b) of Theorem 4.1 with (44) implies that the same fact holds for E δ,t (h). Likewise, the first inequality of (c) in Theorem 4.1 holds true for E δ,t (h).
To end this section, we give an algorithmic description of our method.
Step 1 Pick a point x 0 near ch(Ω). Construct complex spherical waves u 0 h,t (with associated p 0 h,t ).
Step 4 If E δ,t (h) increases to ∞ as h → 0, then the front {|x − x 0 | = t} intersects the obstacle. Decrease t to make more accurate estimate of ∂D.