\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2015 (2015), No. 186, pp. 1--13.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2015 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2015/186\hfil Asymptotic formulas] {Asymptotic formulas for the identification of small inhomogeneities in a fluid medium} \author[M. Abdelwahed, N. Chorfi, M. Hassine \hfil EJDE-2015/186\hfilneg] {Mohamed Abdelwahed, Nejmeddine Chorfi, Maatoug Hassine} \address{Mohamed Abdelwahed \newline Department of Mathematics, College of Sciences, King Saud University, Riyadh, Saudi Arabia} \email{mabdelwahed@ksu.edu.sa} \address{Nejmeddine Chorfi \newline Department of Mathematics, College of Sciences, King Saud University, Riyadh, Saudi Arabia} \email{nejmeddine.chorfi@yahoo.com} \address{Maatoug Hassine \newline D\'epartement de Math\'ematiques, Facult\'e des Sciences de Monastir, Monastir, Tunisia} \email{maatoug\_hassine@yahoo.fr} \thanks{Submitted June 1, 2015. Published July 10, 2015.} \subjclass[2010]{35R30, 35J25, 49Q12, 65R99} \keywords{Stokes system; inhomogeneities; viscous incompressible fluid flow; \hfill\break\indent reciprocity gap functional; asymptotic formula; sensitivity analysis} \begin{abstract} We consider a viscous incompressible fluid flow governed by the Stokes system. We assume that a finite number of small inhomogeneities (particles) are immersed in the fluid. The reciprocity gap functional is introduced to describe the boundary data. An asymptotic formula for the reciprocity gap functional is derived. The obtained formulas can form the basis for very effective computational identification algorithms, aimed at determining information about inhomogeneities from boundary measurements. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \newtheorem{corollary}[theorem]{Corollary} \allowdisplaybreaks \section{Introduction} The problem of determining interior information about a medium form boundary field measurements has received considerable attention in the applied mathematical, as well as in the engineering literature (see \cite{AK1,AMV,AVV,BHJM,KV,SU}). Examples of the later type are found in connection with the identification of cracks \cite{AV,ABV,AB,BBJ,NK}. Significant mathematical results on the determination of one or more small conductivity imperfections inside a conductor of known background conductivity have been established in \cite{BFV,CMV,FV}. The reconstruction of electromagnetic imperfections of small diameter form boundary measurement has been analyzed in \cite{VV}. A rigorous derivation of the leading order boundary perturbation resulting from the presence of a finite number of interior particles of small diameter for full Maxwell equations is provided in \cite{AVV}. A boundary integral formula for the reconstruction of imperfections of small diameter in an elastic medium is derived in \cite{AA,AKNT}. This work concerns the fluid mechanics area. Our aim is to design an efficient method to determine the location and size of a finite number of inhomogeneities of small volume immersed in a fluid medium using boundary measurements. The proposed method is based on a sensitivity analysis of the \emph{reciprocity gap functional} \cite{AB,BBJ} with respect to the presence of a small inhomogeneity. An asymptotic formula is derived giving the relation between the known boundary data (via the reciprocity gap functional) and the unknown inhomogeneities properties; location, size and shape. To present the leading term of the reciprocity gap functional variation we introduce the concept of \emph{Viscous Moment Tensor}. The concept is defined in away analogous to the polarization tensors in electro-magnetic \cite{AMV} and the elastic moment tensors in elasticity \cite{AK1}. The obtained asymptotic formula will serve as very useful tools for developing very effective algorithms for reconstruction of small inhomogeneities from boundary measurements. Such algorithms can be used in various applications like fiber-reinforced polymers \cite{CK,FPZ}, colloid \cite{G,P,ZP} and casting or injection filling \cite{CSO,DGB,S} where the design of the mixing of liquid metallic should be optimized. This article is organized as follows. In the next section we present the governing equations. The considered fluid is viscous and incompressible. The Stokes system is used to describe the fluid motion. In Section \ref{reciprocity}, we introduce the reciprocity gap functional and we establish a preliminary estimate. The main result is presented in Section \ref{results}. We introduce the Viscous Moment Tensor and we derive the asymptotic formulas. The case of single inhomogeneity is discussed in Section \ref{s-particle}. The case of multiples inhomogeneities is considered in Section \ref{m-particles}. Section \ref{proof} is devoted to the proof of the main result. The proof is based on some preliminary Lemmas. We complete this article with some concluding remarks in the last section. \section{Governing equations} \label{equations} Consider a viscous incompressible fluid occupying a bounded and smooth domain of $\mathbb{R}^{d}$, $d=2,3$ with a smooth and connected boundary $\Gamma =\partial \Omega$. We assume that a finite number of immiscible liquid inhomogeneities (particles) $\mathcal{F}^i$, $i=1,\dots ,m$ of small volume $\omega_\varepsilon^i \subset \mathbb{R}^{d}$ are immersed in the fluid. %see figure \ref{inhomg1}. For simplicity, we assume that the inhomogeneities are well separated and have constant physical properties. With each inhomogeneity $\mathcal{F}^i$ we associate its density $\rho_i$, its kinematic viscosity $\nu_i$ and its geometry form $\omega_\varepsilon^i = z_i +\varepsilon \omega^i$, where $\varepsilon$ is the common order of magnitude of the diameter of the inhomogeneities and $\omega^i \subset \mathbb{R}^d$ is a bounded, smooth domain containing the origin. The domains $\omega^i$ determine the relative size and shape of the inhomogeneities. The total collection of inhomogeneities thus takes the form $\omega_\varepsilon = \cup_{i=1}^{m} (z_i+\varepsilon \omega^i)$. The points $z_i \in \Omega$, $i=1,\dots ,m$ determine the location of the inhomogeneities. We shall assume that they satisfy \begin{equation}\label{a0p} | z_i-z_j| \ge d_0 >0, \quad \forall j\neq i \text{ and } \operatorname{dist}(z_i,\partial \Omega)\ge d_0 >0,\; i=1,\dots ,m. \end{equation} We also assume that the parameter $\varepsilon$ is sufficiently small so that the inhomogeneities are disjoint and their distance to $\mathbb{R}^d \backslash {\overline{\Omega}}$ is larger than $d_0/2$. Let $\rho >0$ and $\nu >0$ denote the density and the kinematic viscosity of the background fluid. We assume that \begin{equation}\label{a0} 0 < c_0 \leq \rho=\rho (x)\leq C_0 <\infty.\quad 0 < c_1 \leq \nu=\nu (x) \leq C_1 <\infty \quad \forall x \in \Omega, \end{equation} for some fixed constants $c_0$, $C_0$, $c_1$ and $C_1$. For simplicity, we assume that $\rho$ and $\nu $ are $\mathcal{C}^{\infty}(\overline{\Omega})$, but this latter assumption could be considerably weakened. Using this notation, we introduce the density and the viscosity \[ \rho_\varepsilon = \begin{cases} \rho & \text{if } x \in \Omega_\varepsilon = \Omega \backslash {\overline{\omega_\varepsilon}}\\ \rho_i & \text{if } x \in \omega^i_\varepsilon=z_i+\varepsilon \omega^i, \end{cases} \quad \nu_\varepsilon= \begin{cases} \nu & \text{if } x\in \Omega_\varepsilon = \Omega \backslash {\overline{\omega_\varepsilon}}\\ \nu_i & \text{if } x\in \omega^i_\varepsilon=z_i+\varepsilon \omega^i. \end{cases} \] In this work, we assume that both the continuous phase (the background fluid) and the dispersed phase (the inhomogeneities $\mathcal{F}^i$, $i=1,\dots ,m$) are immiscible viscous incompressible fluids governed by the Stokes equations. Just the gravitational force is considered. We assume that $\Gamma$ is partitioned into two parts $\Gamma_d$ and $\Gamma_n$ such that $\Gamma=\Gamma_d \cup \Gamma_n $ and $\Gamma_d \cap \Gamma_n = \emptyset$. In the presence of the inhomogeneities, the velocity vector $u_\varepsilon$ and the pressure $p_\varepsilon$ satisfy \begin{equation}\label{a1} \begin{gathered} - \nabla \cdot [2 \nu_\varepsilon \mathcal{D}(u_\varepsilon)] +\nabla p_\varepsilon = \rho_\varepsilon \mathcal{G} \quad \text{in } \Omega\\ \nabla \cdot u_\varepsilon = 0 \quad\text{in } \Omega\\ u_\varepsilon = u_d \quad \text{on } \Gamma_d\\ \sigma (u_{\varepsilon},p_\varepsilon) \mathbf{n} = g \quad \text{on } \Gamma_n, \end{gathered} \end{equation} where $\mathcal{G}$ is the constant gravity acceleration, $\mathcal{D}(u_\varepsilon) = \frac{1}{2} \big (\nabla u_\varepsilon + \nabla u_\varepsilon^{T} \big )$ is the rate of deformation tensor for the flow, $\sigma (u_{\varepsilon},p_\varepsilon) = 2 \nu_\varepsilon \mathcal{D}(u_\varepsilon) - p_\varepsilon I$ is the stress tensor, $I$ is the $d\times d$ identity matrix, $u_d$ is a given velocity on $\Gamma_d$ and $g$ is a given traction force exerted on the free surface $\Gamma_n$. As physical interpretation, the quantity $(\sigma (u_{\varepsilon},p_\varepsilon)\mathbf{n})(x)$ is the force at a point $x \in \partial \Omega$ which acts on the fluid in $\Omega$. Here, $\mathbf{n}$ denotes the unit normal to the boundaries $\partial \Omega$ and $\partial \omega_\varepsilon$ which is outer with respect to $\Omega$ and $\omega_\varepsilon$. In the absence of any inhomogeneities, the velocity $u_0$ and the pressure $p_0$ satisfy \begin{equation}\label{a3} \begin{gathered} - \nabla \cdot [2\nu \mathcal{D}(u_0)] +\nabla p_0 = \rho \mathcal{G} \quad\text{in } \Omega\\ \nabla \cdot u_0 = 0 \quad\text{in } \Omega\\ u_{0} = u_d \quad \text{on } \Gamma_d\\ \sigma (u_{0},p_0) \mathbf{n} = g \quad \text{on } \Gamma_n. \end{gathered} \end{equation} Alternatively, \eqref{a1} may be formulated as \begin{equation}\label{a2} \begin{gathered} - \nabla\cdot [2\nu \mathcal{D}(u_\varepsilon)] +\nabla p_\varepsilon = \rho \mathcal{G} \quad \text{in } \Omega\backslash {\overline{\omega_\varepsilon}}\\ \nabla\cdot u_\varepsilon = 0 \quad\text{in } \Omega\backslash {\overline{\omega_\varepsilon}}\\ - \nabla\cdot [2\nu_i \mathcal{D}(u_\varepsilon)] +\nabla p_\varepsilon = \rho_i \mathcal{G} \quad \text{in } z_i+\varepsilon \omega^i,\; i=1,\dots ,m\\ \nabla\cdot u_\varepsilon = 0 \quad\text{in } z_i+\varepsilon \omega^i,\; i=1,\dots ,m. \end{gathered} \end{equation} The above equations are to be solved subject to appropriate boundary conditions. At the exterior boundary $\Gamma$ we consider the boundary conditions described in \eqref{a1}. At the inhomogeneity surface $\partial \omega_\varepsilon^i$, kinematic and stress boundary conditions are imposed. The kinematic boundary condition, which stipulates the continuity of velocities at the interface, is \[ {u_\varepsilon ^+}_{|\partial \omega_\varepsilon^i} (\text{exterior fluid}) = {u_\varepsilon ^-} _{|\partial \omega_\varepsilon^i} (\text{interior fluid}) \text{ on } \partial \omega_\varepsilon^i,\; i=1,\dots ,m . \] The stress boundary condition requires that mechanical equilibrium be satisfied at the interface. The stress exerted on the interface are the hydrodynamics forces resulting from the interior and exterior fluids. Neglecting surface tension effects (at low Reynold number), the stress boundary condition at the interface is therefore \[ \big\{ 2\nu (x) \mathcal{D}(u_\varepsilon) - p_\varepsilon I \big\}^+ \mathbf{n} = \big\{ 2 \nu_i \mathcal{D}(u_\varepsilon) - p_\varepsilon I \big\} ^- \mathbf{n} \quad \text{on } \partial \omega_\varepsilon^i,\; i=1,\dots ,m \,. \] \section{Reciprocity gap functional}{\label{reciprocity}} We introduce the reciprocity gap functional, one can see \cite{AB} for Laplace equation or \cite{AS} for Stokes system. It is based on the boundary data. For the Stokes equations, we have $\mathcal{R}_\varepsilon :H^1(\Omega)^d\times L^2(\Omega) \to \mathbb{R}$ by \[ \mathcal{R}_\varepsilon(w,\xi)= \int_{\partial \Omega} \sigma(w,\xi)\mathbf{n} u_\varepsilon \,ds - \int_{\partial \Omega}\sigma(u_\varepsilon,p_\varepsilon)\mathbf{n} w\,ds, \] where $(u_\varepsilon,p_\varepsilon)$ is the solution to \eqref{a2}. If we consider the restriction of the reciprocity gap function to the space \[ \mathcal{V} (\Omega) =\Big\{ (w,\xi) \in H^1(\Omega)^d\times L^2(\Omega);\, - \nabla\cdot [2 \nu \mathcal{D}(w)] +\nabla \xi =0,\, \nabla\cdot w = 0 \text{ in }\Omega \Big\} \] we obtain the following preliminary estimate. \begin{proposition} \label{Prop-1} For all $(w,\xi) \in \mathcal{V} (\Omega)$, we have \[ \mathcal{R}_\varepsilon(w,\xi)=\mathcal{R}_0(w,\xi) + \int_{\omega_\varepsilon} 2 [\nu-\nu_\varepsilon] \mathcal{D}(u_\varepsilon):\mathcal{D}(w)\,dx - \int_{\omega_\varepsilon} (\rho-\rho_\varepsilon)\mathcal{G} w \,dx. \] \end{proposition} \begin{proof} Using Green's formula, \eqref{a1} and \eqref{a3} we obtain \begin{gather*} \int_{\Omega} 2\nu_\varepsilon \mathcal{D}(u_\varepsilon):\mathcal{D}(w) \,dx = \int_{\Omega} \rho_\varepsilon \mathcal{G}\, w \,dx + \int_{\Gamma} \sigma(u_\varepsilon,p_\varepsilon) \mathbf{n} \,w \,ds\quad \forall (w,\xi) \in \mathcal{V} (\Omega),\\ \int_{\Omega} 2\nu \mathcal{D}(u_0):\mathcal{D}(w) \,dx = \int_{\Omega} \rho \mathcal{G}\, w \,dx + \int_{\Gamma} \sigma(u_0,p_0)\mathbf{n} w \,ds \quad \forall (w,\xi) \in \mathcal{V} (\Omega). \end{gather*} From the fact that $- \nabla\cdot [2 \nu \mathcal{D}(w)] +\nabla \xi =0$, and $ \nabla\cdot w = 0$ in $\Omega$, one gets \begin{gather*} \int_{\partial \Omega} \sigma(w,\xi)\mathbf{n} u_\varepsilon \,ds =\int_{\Omega} 2\nu \mathcal{D}(w): \mathcal{D}(u_\varepsilon)\,dx \\ \int_{\partial \Omega} \sigma(w,\xi)\mathbf{n} u_0 \,ds = \int_{\Omega} 2\nu \mathcal{D}(w) :\mathcal{D}(u_0)\,dx. \end{gather*} Combining the previous equalities, it follows that \[ \mathcal{R}_\varepsilon(w,\xi)=\mathcal{R}_0(w,\xi) + \int_{\omega_\varepsilon} 2 [\nu-\nu_\varepsilon] \mathcal{D}(w):\mathcal{D}(u_\varepsilon) \,dx - \int_{\omega_\varepsilon} (\rho-\rho_\varepsilon)\mathcal{G} \,w \,dx, \] for all $(w,\xi) \in \mathcal{V} (\Omega)$. \end{proof} Let $(U(.,z),\,P(.,z)) \in (\mathbb{R}^d\times\mathbb{R}^d )\times\mathbb{R}^d$ denote the fundamental solution to the Stokes equations corresponding to a Dirac mass at the point $z$ and to coefficient $\nu$. That is, for all $1\le j\le d$, $(U^j(.,z),\,P^j(.,z))$ is a solution to \begin{equation}\label{a6} \begin{gathered} - \nabla_x \cdot [2\nu (x) \mathcal{D}_x(U^j)(x,z)] +\nabla_x P^j(x,z) = \delta_z e_j \quad\text{in } \Omega,\\ \nabla_x \cdot U^j(x,z) = 0 \quad \text{in } \Omega, \end{gathered} \end{equation} where $U^j$ denotes the $j^{th}$ column of $U$. \begin{remark} \label{rmk3.1} \rm In the case where $\nu$ is constant (see \cite{DL}), the fundamental solution ($U,\,P$) is given by \begin{gather*} U(x,z) = \frac{1}{4\pi \nu} \big( -\log (r) I + e_r e_r^T \big), \quad P(x,z) = \frac {x}{2 \pi r^2} \quad \text{if } d=2,\\ U(x,z) = \frac{1}{8\pi \nu r} \big( I + e_r e_r^{T} \big), \quad P(x,z) = \frac {x}{4 \pi r^3}, \quad \text{if } d=3, \end{gather*} with $r=\|x-z\|$, $ e_r=x/r$ and $e_r^{T}$ is the transposed vector of $e_r$. \end{remark} A \emph{Stokeslet} of strength $\mathbf{b} \in \mathbb{R}^d$ located at the point $z\in \mathbb{R}^d$ is a function of $x $ formed by the pair $(U(x,z)\mathbf{b},P(x,z)\cdot \mathbf{b})$. Let $\mathcal{S} (\Omega)$ be the following set, obtained by restriction to $\Omega$ of Stokeslets located at points in the exterior of $\Omega$ \[ \mathcal{S} (\Omega) = \{ (U(x,z)\mathbf{b},P(x,z).\mathbf{b}),\; \mathbf{b}\in \mathbb{R}^d,\; z\in \mathbb{R}^d\backslash \overline{\Omega}\}. \] It is clear that $\mathcal{S} (\Omega)$ is a subset of $\mathcal{V} (\Omega)$. For each fixed $1\leq j \leq d$, we denote by $\mathcal{R}^{j}_\varepsilon$ the following reciprocity gap functional associated to Stokeslets of strength $\mathbf{e}_j$, \[ \mathcal{R}^{j}_\varepsilon(z) = \int_{\partial \Omega} \sigma \big( U^j(x,z),P^j(x,z)\big)\mathbf{n} u_\varepsilon \,ds - \int_{\partial \Omega}\sigma(u_\varepsilon,p_\varepsilon)\mathbf{n} U^j(x,z)\,ds, \] for all $z\in \mathbb{R}^d\backslash \overline{\Omega}$. By Proposition \ref{Prop-1}, we deduce that the unknown parameters $m$, $z_k$ and $\omega^k$, $k=1,\dots ,m$ must satisfy the system. \begin{corollary}\label{corol61} For each $1\leq j \leq d$, the reciprocity gap function $\mathcal{R}^{j}_\varepsilon$ satisfies the expansion \begin{equation}{\label{syst-0}} \begin{aligned} \mathcal{R}^{j}_\varepsilon(z)-\mathcal{R}^{j}_0(z) &= \sum_{k=1}^m \int_{z_k+\varepsilon\omega^k} 2[\nu-\nu_\varepsilon] \mathcal{D}_x(U^j(x,z)) :\mathcal{D}(u_\varepsilon)\,dx \\ &\quad - \sum_{k=1}^m \int_{z_k+\varepsilon\omega^k} (\rho- \rho_\varepsilon)\mathcal{G} U^j(x,z) \,dx,\quad \forall z\in \mathbb{R}^d\backslash \overline{\Omega}. \end{aligned} \end{equation} \end{corollary} \section{Main results} \label{results} In this section we derive an asymptotic formula for the reciprocity gap function $\mathcal{R}^j_\varepsilon$. The obtained results will serve as very useful tools for the numerical reconstruction of the ``location'' and ``size'' of the inhomogeneities. We shall initially consider the case in which $\Omega$ contains a single inhomogeneity $\omega_\varepsilon=\varepsilon \omega$, that is centred at the origin. The case where $\Omega$ contains multiple inhomogeneities is presented in section \ref{m-particles}. \subsection{Single inhomogeneity} \label{s-particle} First, we introduce the concept of the Viscous Moment Tensor. \begin{definition} \label{def4.1} \rm The viscous moment tensor associated to the domain $\omega$ and viscosity ratio $ \nu (0)/\nu_1$ is given by \[ \mathcal{M}_{pq}^{kl}=(\frac{\nu(0)}{\nu_1}-1) \int_{\partial \omega }y_pe_q \Big ( E^{k,l}\mathbf{n}+ (1-\frac{\nu_1}{\nu(0)}) [\frac{\nu(0)}{\nu_1}\mathcal{D}_y(v^{k,l})- q^{k,l} I ]^+ \mathbf{n}\Big )\,ds(y), \] for $1\le k,l,p,q\le d$, where $(e_q)_{q=1}^{d}$ is the canonical basis of $\mathbb{R}^d$, $y_p$ denotes the $p^{th}$ component of $y$, and ($v^{k,l},q^{k,l}$), denotes the solution to \begin{equation} \label{b2p} \begin{gathered} - \nabla_y \cdot[\frac{\nu (0)}{\nu_1} \mathcal{D}_y(v^{k,l})] +\nabla_y q^{k,l} = 0,\quad \nabla_y \cdot v^{k,l} = 0, \quad \text{in } \mathbb{R}^d\backslash {\overline{\omega}} \\ - \nabla_y \cdot [\mathcal{D}_y(v^{k,l})] +\nabla_y q^{k,l} = 0,\quad \nabla_y\cdot v^{k,l} = 0 , \quad \text{in } \omega \\ v^{k,l} \text{ is continuous across } \partial \omega, \\ \Big( \frac{\nu (0)}{\nu_1} \mathcal{D}_y(v^{k,l})-q^{k,l} I\Big)^+ \mathbf{n} - \Big( \mathcal{D}_y(v^{k,l})-q^{k,l} I\Big)^- \mathbf{n}= -E^{k,l}\mathbf{n} \quad \text{ on } \partial \omega, \\ \lim_{|y|\to + \infty} v^{k,l}(y)=0, \end{gathered} \end{equation} where $E^{kl} \in \mathbb{R}^d\times \mathbb{R}^d$, $1\le k,l\le d$ is the symmetric matrix defined by \[ E^{kl}_{pq}= \frac{1}{2}(\delta_{pk} \delta_{ql} + \delta_{pl} \delta_{qk}), \quad 1\le p,q\le d. \] Here $\delta_{pq}$ denotes the Kronecker symbol. \end{definition} Since $v^{k,l}$ is continuous across $\partial \omega$, the solution ($v^{k,l},\,q^{k,l}$) can be represented with the help of a single layer potential (see e.g.\cite{DL}), namely there exists $\eta^{k,l} \in H^{-1/2} (\partial \omega)^d$ such that \[ v^{k,l}(y) = \int_{\partial \omega}U(y-x)\eta^{k,l}(x) \,ds(x),\quad q^{k,l}(y) = \int_{\partial \omega}P(y-x).\eta^{k,l}(x) \,ds(x). \] In general the functions $v^{k,l}$ and $q^{k,l}$ cannot be computed explicitly. One exception in the case when $\omega$ is a ball. The main result of this case is given by the following theorem. We derive an asymptotic formula connecting the inhomogeneity properties and the reciprocity gap functional variation. The proof is relegated to section \ref{proof}. \begin{theorem} \label{asympt-1} For all $z\in \mathbb{R}^d\backslash \overline{\Omega}$, we have \begin{equation} \label{asymp-th1} \begin{aligned} \mathcal{R}^{j}_\varepsilon(z)-\mathcal{R}^{j}_0(z) &= \varepsilon^d \Big \{ 2\nu(0) \,\mathcal{D}_x(U^j)(0,z):\mathcal{M}\mathcal{D}_x(u_0)(0) \\ &\quad - [\rho(0) - \rho_1] \,|\omega|\mathcal{G} U^j(0,z)\Big \} + o(\varepsilon^{d}),\quad j=1,\dots ,d. \end{aligned} \end{equation} \end{theorem} Note that the viscous moment tensor $\mathcal{M}$ depends on the viscosity ratio, the size and the shape of the inhomogeneities. The notion of polarization matrix has been introduced by Schiffer and Szeg\"o \cite{SC}, and since it has been extensively studied (see e.g.\cite{AK1} and the references therein). In particular, one can prove here that $\mathcal{M}$ is positive definite and symmetric in the following sense \[ \mathcal{M}^{pq}_{kl}=\mathcal{M}^{qp}_{kl},\quad \mathcal{M}^{pq}_{kl}= \mathcal{M}^{pq}_{lk}, \quad \mathcal{M}^{pq}_{kl}=\mathcal{M}^{kl}_{pq},\quad \forall p,q,k,l \in \{ 1,\dots ,d\}. \] For similar results and proofs one can consult \cite{CMV} for the conductivity problem and \cite{AK1} for the elasticity system. \begin{remark} \label{rmk4.1} \rm Using Green's formula and the jump relation on $\partial \omega$ one can derive the following expressions of $\mathcal{M}$: \begin{align*} \mathcal{M}_{pq}^{kl} &= (\frac{\nu(0)}{\nu_1}-1) \Big \{ \frac{1}{2} |\omega|(\delta_{pk} \delta_{ql} + \delta_{pl} \delta_{qk}) \\ &\quad + (1-\frac{\nu_1}{\nu(0)}) \int_{\partial \omega }y_p e_q [\frac{\nu(0)}{\nu_1}\mathcal{D}_y(v^{k,l})- q^{k,l} I ]^+ \mathbf{n}ds(y)\Big \} \\ &= (1-\frac{\nu_1}{\nu(0)}) \Big \{ \frac{1}{2} |\omega|(\delta_{pk} \delta_{ql} + \delta_{pl} \delta_{qk}) \\ &\quad + (\frac{\nu(0)}{\nu_1}-1) \int_{\partial \omega }y_p\,e_q [\mathcal{D}_y(v^{k,l})- q^{k,l} I ]^- \mathbf{n}ds(y)\Big \} \\ &= (\frac{\nu(0)}{\nu_1}-1) \int_{\partial \omega }y_p e_q \Big ([\mathcal{D}_y(v^{k,l})- q^{k,l} I ]^- \mathbf{n} \\ &\quad -\frac{\nu_1}{\nu(0)} [\frac{\nu(0)}{\nu_1}\mathcal{D}_y(v^{k,l})- q^{k,l} I ]^+ \mathbf{n}\Big )\,ds(y). \end{align*} \end{remark} \subsection{Multiple inhomogeneities} \label{m-particles} In the case of more than one inhomogeneity, say $\omega_\varepsilon = \cup_{i=1}^m \{ z_i+\varepsilon \omega^i\}$, the previous theorem may very simply be changed to proceed inductively one inhomogeneity at time. The result is described by the following Theorem. \begin{theorem} \label{asympt-2} For all $z\in \mathbb{R}^d\backslash \overline{\Omega}$, we have \begin{equation} \label{asymp-th2} \begin{aligned} \mathcal{R}^{j}_\varepsilon(z)-\mathcal{R}^{j}_0(z) &=\varepsilon^d \sum_{i=1}^{m}\Big \{ 2\nu(z_i) \mathcal{D}_x(U^j)(z_i,z):\mathcal{M}_i\mathcal{D}_x(u_0)(z_i)\\ &\quad - [\rho(z_i) - \rho_i] |\omega|\mathcal{G} U^j(z_i,z)\Big \} + o(\varepsilon^{d}),\quad j=1,\dots ,d, \end{aligned} \end{equation} where $\mathcal{M}_i$ is the viscous moment tensor corresponding to the domain $\omega_i$ and viscosity ratio $\nu(z_i)/\nu_i$. \end{theorem} For each $z\in \mathbb{R}^d\backslash \overline{\Omega}$, the variation $\mathcal{R}^{j}_\varepsilon(z)-\mathcal{R}^{j}_0(z)$ can be entirely estimated from numerical computation of the pair ($ U^j(x,z),\sigma \Big( U^j(x,z),P^j(x,z)\Big)\mathbf{n}$) and boundary measurement of the velocity $u_\varepsilon$ and the stress tensor $\sigma(u_\varepsilon,p_\varepsilon)$. Neglecting the smaller order term, the equation \eqref{asymp-th2} leads to a nonlinear system satisfied by the unknown parameters $m$, $z_k$ and $\omega^k$, $k=1,dots ,m$. \section{Proof of the main result} \label{proof} To prove the main result, we introduced in section \ref{prelim} some preliminary lemmas. The proof is presented in section \ref{proof-th}. Whenever no confusion is possible we shall use the simpler notation $U^j(x) =U^j(x,z)$ and $P^j(x) =P^j(x,z)$. \subsection{Preliminary estimates} \label{prelim} \begin{lemma} \label{prelim-1} There exists a positive constant $C$, independent of $\varepsilon$, such that for all $j=1,\dots,d$ \[ \big| \int_{\omega_\varepsilon} [\nu-\nu(0)]\mathcal{D}_x(U^j):\mathcal{D}_x(u_\varepsilon)\, dx \big| \leq C \varepsilon^{d+1}. \] \end{lemma} \begin{proof} Expanding $\nu(x)=\nu(0)+x\cdot \nabla_x \nu(\eta_x)$, $\eta_x\in \Omega$, and using the change of variable $x=\varepsilon y$, it follows that \[ \int_{\omega_\varepsilon} [\nu-\nu(0)]\mathcal{D}_x(U^j):\mathcal{D}_x(u_\varepsilon){ dx} = \varepsilon^{d+1} \int_{\omega} [y.\nabla_x \nu(\eta_x)]\mathcal{D}_x(U^j)(\varepsilon y ) :\mathcal{D}_x(u_\varepsilon)(\varepsilon y )\,dy . \] From the fact that $x_i\to\nabla_x \nu(x_i)$ is uniformly bounded on $\Omega$ we deduce \[ \big|\int_{\omega_\varepsilon} [\nu-\nu(0)]\mathcal{D}_x(U^j):\mathcal{D}_x(u_\varepsilon)\, dx\big| \le C \varepsilon^{d+1} \|\mathcal{D}_x (U^j)(\varepsilon y)\|_{L^2(\omega)} \|\mathcal{D}_x (u_\varepsilon)(\varepsilon y)\|_{L^2(\omega)} . \] Using Green's formula and equations \eqref{a1} and \eqref{a3}, we have \begin{align*} &\int_{\Omega} 2\nu_\varepsilon |\mathcal{D}(u_\varepsilon-u_0)|^2{ dx}\\ &=\int_{\omega_\varepsilon} 2(\nu-\nu_\varepsilon)\mathcal{D}(u_0):\mathcal{D}(u_\varepsilon-u_0){ dx}- \int_{\omega_\varepsilon} (\rho-\rho_\varepsilon) \mathcal{G} \,(u_\varepsilon-u_0){ dx}. \end{align*} From the fact that $\nu$ and $\rho$ are uniformly bounded on $\Omega$ (due to hypotheses \eqref{a0}), and $u_0$ and $\mathcal{D}(u_0)$ are uniformly bounded on $\omega_\varepsilon$ (due to elliptic regularity), one can obtain \[ \int_{\omega_\varepsilon} |\mathcal{D}_x(u_\varepsilon-u_0)|^2\, dx \le C\,\varepsilon^{d}. \] Changing variable, we have \[ \int_{\omega} |\mathcal{D}_x(u_\varepsilon-u_0)(\varepsilon y)|^2\, dy \le C, \] and therefore, \[ \int_{\omega} |\mathcal{D}_x (u_\varepsilon)(\varepsilon y)|^2\, dy \le C. \] Since $\operatorname{dist}(z,\omega_\varepsilon) \ge d_0>0$, it follows that \[ |\int_{\omega_\varepsilon} [\nu-\nu(0)]\mathcal{D}_x(U^j):\mathcal{D}_x(u_\varepsilon){ dx}| \le C\,\varepsilon^{d+1} \|\mathcal{D}_x U^j(\varepsilon y,z)\|_{L^2(\omega)} \le C\,\varepsilon^{d+1}. \] \end{proof} \begin{lemma} \label{prelim-2} We have the following asymptotic expansion \begin{align*} &\int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(U^j)(x,z):\mathcal{D}_x(u_\varepsilon)(x)\,dx\\ & = \varepsilon^{d} \Big\{\int_{\partial \omega }2\nu(0)\mathcal{D}_x(u_0 )(0)\mathbf{n} \mathcal{D}_x(U^j)(0,z)y \, ds(y) \\ &\quad +\int_{\partial \omega }[2\nu(0) \mathcal{D}_y(v )- q I ]^+(y) \mathbf{n} \mathcal{D}_x(U^j)(0,z)y \, ds(y) \Big \} + O \big(\varepsilon^{d+1/2}\big), \quad j=1,\dots,d. \end{align*} \end{lemma} \begin{proof} Let ($v,\,q$) denote the solution to \begin{equation}\label{b2} \begin{gathered} - \nabla_y\cdot [2\nu (0) \mathcal{D}_y(v)] +\nabla_y q = 0,\quad \nabla_y\cdot v = 0, \quad \text{in } \mathbb{R}^d\backslash {\overline{\omega}} \\ - \nabla_y\cdot [2\nu_1 \mathcal{D}_y(v)] +\nabla_y q = 0,\quad \nabla_y\cdot v = 0 , \quad \text{in } \omega \\ v \text{ is continuous across } \partial \omega, \\ \begin{aligned} &\Big( 2\nu (0) \mathcal{D}_y(v)-q I\Big)^+ \mathbf{n} - \Big( 2\nu_1 \mathcal{D}_y(v)-q I\Big)^- \mathbf{n}\\ &= -2 [\nu (0)-\nu_1]\mathcal{D}_x(u_0)(0)\mathbf{n} \quad \text{on } \partial \omega, \end{aligned} \\ \lim_{|y|\to + \infty} v(y)=0. \end{gathered} \end{equation} The existence of ($v,q$) can be established in a manner similar to that of ($v^{k,l},q^{k,l}$). Setting \[ w_\varepsilon (x) = u_\varepsilon(x)-u_0(x)-\varepsilon\, v(x/\varepsilon),\quad s_\varepsilon(x) = p_\varepsilon(x)-p_0(x)- \,q(x/\varepsilon), \] we have \begin{equation} \label{a25} \begin{aligned} &\int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(U^j):\mathcal{D}_x(u_\varepsilon)\,dx \\ &= \int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(w_\varepsilon): \mathcal{D}_x(U^j)\, dx + \int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(u_0): \mathcal{D}_x(U^j) \,dx \\ &\quad +\varepsilon \int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(v)(x/\varepsilon): \mathcal{D}_x(U^j) \, dx . \end{aligned} \end{equation} Now we shall estimate each term on the right hand side of \eqref{a25} separately. Using the change of variable $x=\varepsilon y $, the first term can be written as \[ \int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(w_\varepsilon): \mathcal{D}_x(U^j)\,dx = \varepsilon^{d-1}\int_{\omega} 2 \nu_1 \mathcal{D}_y(w_\varepsilon)(\varepsilon y ): \mathcal{D}_x(U^j)(\varepsilon y )\,dy \, . \] With the help of the Green's formula and the fact that ($w_\varepsilon,s_\varepsilon$) is solution of \begin{gather*} - \nabla_x\cdot [2 \nu_1 \mathcal{D}_x(w_\varepsilon)] +\nabla_x s_\varepsilon = - \nabla_x\cdot [2 [\nu-\nu_1] \mathcal{D}_x(u_0) ] - (\rho -\rho_1 )\mathcal{G} \quad \text{in } \omega_\varepsilon,\\ \nabla_x\cdot w_\varepsilon = 0 \quad \text{in } \omega_\varepsilon, \end{gather*} one can derive that \[ \| \mathcal{D}_y(w_\varepsilon)(\varepsilon y ) \|_{L^2 (\omega)}=O(\varepsilon^{3/2}). \] Then, using the previous estimate and the fact that $\mathcal{D}_x(U^j)(\varepsilon y ) = \mathcal{D}_x(U^j)(\varepsilon y,z ) $ is uniformly bounded on $\omega$, we obtain \begin{equation} \label{a26} \begin{gathered} \big| \int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(w_\varepsilon): \mathcal{D}_x(U^j)\,dx \big| \leq \varepsilon^{d-1} \| \mathcal{D}_y(w_\varepsilon)(\varepsilon y ) \|_{L^2 (\omega)} \|\mathcal{D}_x(U^j)(\varepsilon y ) \|_{L^2 (\omega)}\\ \le C \varepsilon^{d+1/2} . \end{gathered} \end{equation} Expanding the functions $\mathcal{D}_x(U^j)$ and $\mathcal{D}_x(u_0)$ in a Taylor series about the origin, we see that the second term in \eqref{a25} may be written as \begin{equation} \label{a27} \begin{aligned} \int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(u_0): \mathcal{D}_x(U^j) \,dx &= \varepsilon^d \int_{\omega} 2 \nu_1 \mathcal{D}_x(u_0)(\varepsilon y ) : \mathcal{D}_x(U^j)(\varepsilon y,z ) { dy }\\ &= 2 \nu_1 \,|\omega|\,\varepsilon^d \, \mathcal{D}_x(u_0)(0):\mathcal{D}_x(U^j)(0,z) + O\big( \varepsilon^{d+1}\big). \end{aligned} \end{equation} To estimate the third term, we again use the change of variable $x=\varepsilon y $, Taylor's theorem and the Green's formula \begin{align*} &\varepsilon \int_{\omega_\varepsilon }2\nu_1 \mathcal{D}_x(v )(x/\varepsilon) :\mathcal{D}_x(U^j)\,dx\\ &= \int_{\omega_\varepsilon }2\nu_1 \mathcal{D}_y(v ) :\mathcal{D}_x(U^j){ dx} \\ &= \int_{\omega_\varepsilon }2\nu_1 \mathcal{D}_y(v ) :\mathcal{D}_x(U^j)(0)\,dx + \int_{\omega_\varepsilon }2\nu_1 \mathcal{D}_y(v ) :[\mathcal{D}_x(U^j)(x)-\mathcal{D}_x(U^j)(0)]\,dx\\ &=\varepsilon^{d}\int_{\partial \omega}[2\nu_1 \mathcal{D}_y(v )- q I ]^- \mathbf{n} \mathcal{D}_x(U^j)(0)y\, ds(y) + O\big( \varepsilon^{d+1}\big). \end{align*} Using the jump relation on $\partial \omega$ we derive \begin{equation} \label{a28} \begin{aligned} &\int_{\partial \omega_\varepsilon}[2\nu_1 \mathcal{D}_y(v )- q I ]^- (x/\varepsilon) \mathbf{n} U^j(x,z)\, ds(x) \\ &= \varepsilon^{d}\int_{\partial \omega}[2\nu(0) \mathcal{D}_y(v )- q I ]^+ \mathbf{n}\mathcal{D}_x(U^j)(0)y\, ds(y) \\ &\quad + \varepsilon^{d} \int_{\partial \omega} 2(\nu(0)-\nu_1) \mathcal{D}_x(u_0)(0) \mathbf{n} \mathcal{D}_x(U^j)(0)y\, ds(y) + O\big( \varepsilon^{d+1}\big)\\ & = \varepsilon^{d}\int_{\partial \omega}\Big([2\nu(0) \mathcal{D}_y(v )- q I ]^+ \mathbf{n} + 2\nu(0)\mathcal{D}_x(u_0)(0)\mathbf{n}\Big)\mathcal{D}_x(U^j)(0)y\, ds(y)\\ &\quad - 2 \nu_1 |\omega|\varepsilon^d \mathcal{D}_x(u_0)(0): \mathcal{D}_x(U^j)(0,z) + O\big( \varepsilon^{d+1}\big). \end{aligned} \end{equation} Substituting \eqref{a26}, \eqref{a27} and \eqref{a28} in \eqref{a25}, we obtain \begin{align*} &\int_{\omega_\varepsilon} 2 \nu_1 \mathcal{D}_x(U^j)(x,z):\mathcal{D}_x(u_\varepsilon)(x)\,dx \\ & = \varepsilon^{d} \Big\{\int_{\partial \omega }2\nu(0)\mathcal{D}_x(u_0 )(0)\mathbf{n} \mathcal{D}_x(U^j)(0,z)y \,\,ds(y) \\ &\quad -3cm +\int_{\partial \omega }[2\nu(0) \mathcal{D}_y(v )- q I ]^+(y)\mathbf{n} \mathcal{D}_x(U^j)(0,z)y \,ds(y) \Big \} + O \big(\varepsilon^{d+1/2}\big), \end{align*} for $j=1,\dots,d$. \end{proof} \begin{lemma} \label{prelim-3} We have the estimate \[ \int_{\omega_\varepsilon} (\rho -\rho_1) \mathcal{G} U^j{ dx} = \varepsilon^{d} (\rho(0) -\rho_1) |\omega|\mathcal{G} U^j(0,z)+ \mathcal{O}(\varepsilon^{d+1} ). \] \end{lemma} \begin{proof} Expanding $\rho(x)=\rho(0)+x\cdot \nabla_x \rho(\theta_x)$, $\theta_x\in \Omega$, and using the change of variable $x=\varepsilon y$, we have \begin{align*} &\int_{\omega_\varepsilon} (\rho -\rho_1) \mathcal{G} U^j\,dx \\ &= \int_{\omega_\varepsilon} (\rho(0) -\rho_1) \mathcal{G} U^j\,dx + \int_{\omega_\varepsilon} [x.\nabla_x \rho(\theta_x)]\mathcal{G} U^j\,dx \\ &= \varepsilon^{d} \int_{\omega}(\rho(0) -\rho_1)\mathcal{G} U^j(\varepsilon y,z)\, dy + \varepsilon^{d+1} \int_{\omega} [y.\nabla_x \rho(\theta_x)]\mathcal{G} U^j(\varepsilon y,z)\,dy . \end{align*} Due to Taylor's theorem and the fact that $x_i\to\nabla_x \rho(x_i)$ is uniformly bounded on $\Omega$ we derive \begin{equation}\label{p1} \int_{\omega_\varepsilon} (\rho -\rho_1) \mathcal{G} U^j\,dx = \varepsilon^{d} (\rho(0) -\rho_1) |\omega|\mathcal{G} U^j(0,z)+ \mathcal{O}(\varepsilon^{d+1} ). \end{equation} \end{proof} \subsection{Proof of Theorem \ref{asympt-1}} \label{proof-th} By Proposition \ref{Prop-1}, we have \[ \mathcal{R}^j_\varepsilon(z)-\mathcal{R}^j_0(z) =\int_{\omega_\varepsilon} 2[\nu-\nu_1]\mathcal{D}(U^j):\mathcal{D}(u_\varepsilon)\,dx - \int_{\omega_\varepsilon} (\rho -\rho_1) \mathcal{G} U^j\,dx. \] From Lemmas \ref{prelim-1} and \ref{prelim-2}, it follows that \begin{equation} \label{p2} \begin{aligned} &\int_{\omega_\varepsilon} 2[\nu-\nu_1]\mathcal{D}(U^j):\mathcal{D}(u_\varepsilon)\,dx \\ &=\int_{\omega_\varepsilon} 2[\nu-\nu(0)]\mathcal{D}(U^j):\mathcal{D}(u_\varepsilon)\,dx + \int_{\omega_\varepsilon} 2[\nu(0)-\nu_1]\mathcal{D}(U^j):\mathcal{D}(u_\varepsilon)\,dx \\ &=\quad \frac{\nu(0)-\nu_1}{\nu_1}\varepsilon^{d} \Big\{\int_{\partial \omega }2\nu(0)\mathcal{D}_x(u_0 )(0)\mathbf{n}\mathcal{D}_x(U^j)(0,z)y \, ds(y) \\ &\quad +\int_{\partial \omega }[2\nu(0) \mathcal{D}_y(v )- q I ]^+\ mathbf{n}\mathcal{D}_x(U^j)(0,z)y \, ds(y) \Big \} + O \big(\varepsilon^{d+1/2}\big), \end{aligned} \end{equation} for $j=1,\dots,d$. The above equation can be rewritten as \begin{equation} \label{p3} \begin{aligned} &\int_{\omega_\varepsilon} 2[\nu-\nu_1]\mathcal{D}(U^j):\mathcal{D}(u_\varepsilon)\,dx \\ & = \frac{\nu(0)-\nu_1}{\nu_1}\,\varepsilon^{d} \mathcal{D}_x(U^j)(0,z) : \Big\{2\nu(0) \int_{\partial \omega }y \otimes \mathcal{D}_x(u_0 )(0)\mathbf{n}\, ds(y) \\ &\quad + \quad \int_{\partial \omega }y \otimes [2\nu(0) \mathcal{D}_y(v ) - q I ]^+\mathbf{n}\, ds(y) \Big \} + O \big(\varepsilon^{d+1/2}\big),\quad j=1,\dots,d. \end{aligned} \end{equation} Using the definitions of $E^{k,l}$, ($v,\,q$) and ($v^{k,l},q^{k,l}$), it is easy to show that \begin{gather*} \mathcal{D}_x(u_0)(0)=\sum_{1\le k,l\le d} [\mathcal{D}_x(u_0)(0)]_{kl}\,E^{k,l},\\ v(y)=\frac{\nu(0)-\nu_1}{\nu_1} \sum_{1\le k,l\le d} [\mathcal{D}_x(u_0)(0)]_{kl}\,v^{k,l},\\ q(y)=2(\nu(0)-\nu_1) \sum_{1\le k,l\le d} [\mathcal{D}_x(u_0)(0)]_{kl}\,q^{k,l}, \end{gather*} where \[ [\mathcal{D}_x(u_0)(0)]_{kl} = \frac{1}{2} \Big( \frac{\partial u^k_0}{\partial x_l}(0) + \frac{\partial u^l_0}{\partial x_k}(0) \Big). \] The integral term in \eqref{p3} may be decomposed as \begin{equation} \label{p4} \begin{aligned} &2\nu(0)\int_{\partial \omega }y \otimes \mathcal{D}_x(u_0 )(0)\mathbf{n}{ds(y)} +\int_{\partial \omega }y \otimes [2\nu(0) \mathcal{D}_y(v ) - q I ]^+(y) \mathbf{n}{ds(y)}\\ &= 2\nu(0)\sum_{1\le k,l\le d} [\mathcal{D}_x(u_0)(0)]_{kl} \Big\{\int_{\partial \omega }y \otimes E^{k,l}\mathbf{n}\,ds(y) \\ &\quad +\frac{\nu(0)-\nu_1}{\nu(0)} \int_{\partial \omega }y \otimes [\frac{\nu(0)}{\nu_1} \mathcal{D}_y(v^{k,l})- q^{k,l} I ]^+(y) \mathbf{n}\,ds(y) \Big\} \end{aligned} \end{equation} Finally, inserting \eqref{p4} in \eqref{p3} and using Lemma \ref{prelim-3} we conclude that \begin{align*} \mathcal{R}^j_\varepsilon(z)-\mathcal{R}^j_0(z) &=\varepsilon^d \Big \{ 2\nu(0) \mathcal{D}_x(U^j)(0,z):\mathcal{M}\mathcal{D}_x(u_0)(0) \\ &\quad- (\rho(0) - \rho_1) |\omega|\mathcal{G} U^j(0,z)\Big \} + O(\varepsilon^{d+1/2}),\quad j=1,\dots,d. \end{align*} which implies the desired asymptotic formula. \subsection*{Conclusion} %\label{conclusion} The asymptotic formulas derived in Section \ref{results} can serve as very useful tools for the numerical reconstruction of the ``location'' and ``size'' of the inhomogeneities. If for instance ``the normal component of the stress tensor, $\sigma(u_\varepsilon,p_\varepsilon)$ is prescribed on $\Gamma_n$ and measured on $\Gamma_d$'' and ``the velocity field $u_\varepsilon$ is prescribed on $\Gamma_d$ and measured on $\Gamma_n$'', then the function \[ \mathcal{R}^j_\varepsilon(z)-\mathcal{R}^j_0(z)= \int_{\partial \Omega} \sigma(U^j,P^j) \mathbf{n} (u_\varepsilon-u_0) \,ds - \int_{\partial \Omega}\sigma(u_\varepsilon-u_0,p_\varepsilon-p_0)\mathbf{n} U^j\,ds, \] for $j=1,\dots,d$, may be considered as a measured datum on $\partial \Omega$. The parameters $\rho$ and $\nu$ are assumed to be known, and we may easily compute ($u_0,p_0$). From the asymptotic formula in Theorem \ref{asympt-2} it now follows that, up to terms of smaller order, we are in possession of the values of the quantity \[ \sum_{i=1}^{m}\big\{ 2\nu(z_i) \,\mathcal{D}_x(U^j)(z_i,z):\mathcal{M}_i\mathcal{D}_x(u_0)(z_i) - [\rho(z_i) - \rho_i] |\omega|\mathcal{G} U^j(z_i,z)\Big \}, \] for $z\in \mathbb{R}^d\backslash \overline{\Omega}$ and $j=1,\dots,d$. A first task of the identification process, is then to determine (as well as possible) the number $m$ of ``poles'' (centers of inhomogeneities), and their locations $z_k,\, 1\leq k\leq m$. A second task would be to determine other information about the inhomogeneities, such as their sizes. A detailed account of this work and some numerical results illustrating the identification method will be the subject of a forthcoming article. \subsection*{Acknowledgements} The authors would like to extend their sincere appreciation to the deanship of Scientific Research at King Saud University for funding this research group No (RG - 1435-026). \begin{thebibliography}{00} \bibitem{AV} G. Alessandrini, A. Diaz Valenzuela; \emph{Unique determination of multiple cracks by two measurements}, SIAM. J. Control Optim. 34 (3), 1996, 913-921. \bibitem{ABV} G. Alessandrini, E. Beretta, S. Vessela; \emph{Determining linear cracks by boundary measurements: Lipschitz stability}, SIAM J. Math. Anal. 27, 1996, 361-375. \bibitem{AA} C. Alves, and H. Ammari; \emph{Boundary integral formulae for the reconstruction of imperfections of small diameter in an elastic medium}, SIAM, J. Appl. Math. 62, 2001, 94-106. \bibitem{AS} C. Alves, A. L. Silvestre; \emph{On the determination of point-forces on a Stokes system}, Mathematics and computers in Simulation 66 (2004), 385-397. \bibitem{AK1} H. Ammari, H. Kang; \emph{Reconstruction of Small Inhomogeneities from Boundary Measurements}, Lecture Notes in Mathematics, Volume 1846, Springer-Verlag, Berlin 2004. \bibitem{AKNT} H. Ammari, H. Kang, G. Nakamura, K. Tanuma; \emph{Complete asymptotic expansions of solutions of the system of elastostatics in the presence of an inclusion of small diameter and detection of an inclusion}, J. Elasticity 67, 2002, 97-129. \bibitem{AMV} H. Ammari, S. Moskow, M. Vogelius; \emph{Boundary integral formulas for reconstruction of electromagnetic imperfections of small diameter}, ESAIM, Cont. Opt. Cal. Variat. vol. 9, 2004, 49-66. \bibitem{AVV} H. Ammari, M. Vogelius, D. Volkov; \emph{Asymptotic formulas for perturbations in the electromagnetic fields due to presence of inhomogeneities of small diameter II. The full Maxwell's equations}, J. Math. Pures Appl. 80, 2001, 769-814. \bibitem{AB} S. Andrieux, A. Ben Abda; \emph{Identification of planar cracks by complete overdetermined data: inversion formulae}, Inverse Problems, 12, 1996, 553-563. \bibitem{BBJ} A. Ben Abda, H. Ben Ameur, M. Jaoua; \textit{Identification of 2D cracks by elastic boundary measurements}, Inverse Problems 15, 1999, 67-77. \bibitem{BHJM} A. Ben Abda, M. Hassine, M. Jaoua, M. Masmoudi; \emph{Topological sensitivity analysis for the location of small cavities in Stokes flow}, to appear in SIAM J. Cont. Optim. \bibitem{BFV} E. Beretta, E. Francini, M. Vogelius; \emph{Asymptotic formulas for for steady state voltage potentials in the presence of thin inhomogeneities. A rigorous error analysis}, J. Math. Pures Appl. 82, 2003, 1277-1301. \bibitem{CMV} D. J. Cedio-Fengya, S. Moskow, M. Vogelius; \emph{Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction}, Inverse Problems. vol. 14, 1998, 553-595. \bibitem{CK} S. T. Chung, T. H. Kwon; \emph{Numerical simulation of fiber orientation in injection moulding of short-fiber-reinforced thermoplastics,} Polym. Eng. Sci. 7 (35), 1995, 604-618. \bibitem{CSO} R. Codina, U. Schaffer, E. O\~nate; \emph{Mould filling simulation using finite elements}, Int. J. Numer. Meth. Fluid Flow 4, 1994, 291-310. \bibitem{DL} R. Dautray, J. Lions; \emph{Analyse math\'emathique et calcul num\'erique pour les sciences et les techniques}. MASSON, collection CEA, 1987. \bibitem{DGB} G. Dhatt, D. M. Gao, A. Ben Cheikh; \emph{A finite element simulation of metal flow in moulds}, Int. J. Numer. Meth. Eng. 30, 1990, 821-831. \bibitem{FPZ} X. Fan, N. Phan-Thien, R. Zheng; \emph{A direct simulation of fiber suspensions}, J. Non-Newtonian Fluid Mech. 74, 1998), 113-135. \bibitem{FV} A. Friedman, M. Vogelius; \emph{Identification of small inhomogeneties of extreme conductivity by boundary measurements: a theorem on continuous dependence}, Arch. Rat. Mech. Anal., 105, 1989, 299-326. \bibitem{G} T. Gotz; \emph{Simulating particles in Stokes flow}, J. Comput. Appl. Math. 175 (2005), 415-427. \bibitem{KV} R. Kohn, M. Vogelius; \emph{Determining conductivity by boundary measurements}, Comm. Pure Appl. Math., 37 (1984), pp. 289-298, II. Interior results, Comm. Pure Appl. Math., 38 (1985), pp. 643-667. \bibitem{NK} N. Nishimura, S. Kobayashi; \emph{A boundary integral equation method for an inverse problem related to crack detection}, Int. J. Num. Methods Eng., 32(1991), pp. 1371-1387. \bibitem{P} C. Pozrikidis; \emph{Dynamic simulation of the flow of suspensions of two-dimensional particles with arbitrary shapes}, Eng. Anal. Boundary Elements 25 (1) (2001), pp.19-30. \bibitem{SC} M. Schiffer, G. Szeg\"o; \emph{Virtual mass and polarization}, {Amer. Math. Soc.} 67(1949), 130-205. \bibitem{S} K. M. Shyue; \emph{A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state}, J. Comput. Phys. 156 (1999), pp. 43-88. \bibitem{SU} J. Sylvester, G. Uhlmann; \emph{A global uniqueness theorem for inverse boundary value problem}, Annals of Math., 125(1987), pp. 153-169. \bibitem{VV} M. Vogelius, D. Volkov; \emph{Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities}, Math. Model. Numer. Anal. 34 (2000), 723-748. \bibitem{ZP} H. Zhou, C. Pozrikidis; \emph{Adaptive singularity method for Stokes flow past particles}, J. Comput. Phys. 117 (1995), pp. 79-89. \end{thebibliography} \end{document}