# Integral Representation of Hydraulic Permeability

Chuan Bi  Department of Mathematical Sciences, University of Delaware, Newark,
DE 19716, USA  M. Yvonne Ou  Department of Mathematical Sciences, University
of Delaware, Newark, DE 19716, USA  Shangyou Zhang  Department of Mathematical
Sciences, University of Delaware, Newark, DE 19716, USA

###### Abstract

In this paper, we show that the permeability of a porous material [40] and
that of a bubbly fluid [29] are limiting cases of the complexified version of
the two-fluid models posed in [29]. We assume the viscosity of the inclusion
fluid is $z\mu_{1}$ and the viscosity of the hosting fluid is
$\mu_{1}\in\mathbb{R}^{+}$, $z\in\mathbb{C}$. The proof is carried out by the
construction of solutions for large $|z|$ and small $|z|$ with an iteration
process similar to the one used in [16, 21] and the analytic continuation.
Moreover, we also show that for a fixed microstructure, the permeabilities of
these three cases share the same integral representation formula (IRF) (99)
with different values of contrast parameter $s:=1/(z-1)$, as long as $s$ is
outside the interval
$[-\frac{2E_{2}^{2}}{1+2E_{2}^{2}},-\frac{1}{1+2E_{1}^{2}}]$, where the
positive constants $E_{1}$ and $E_{2}$ are the extension constants that depend
only on the geometry of the periodic pore space of the material.

Version of:

Keywords: hydraulic permeability, Stokes equations, Composite materials,
Integral representation formula, Stieltjes class

Classification: 35Q35, 35Q70

##  1 Introduction

Darcy’s law, which was first proposed by H. Darcy in 1856 [19] based on
experimental observation of water flowing through beds of sand, describes the
relationship between the spontaneous flow discharge rate of steady state
through a porous medium, the viscosity of the fluid, and the pressure drop
over a distance. Later, theoretical/mathematical derivations of Darcy’s law
were presented in many works, e.g. M. Poreh et. al [36], S.P. Neuman[32], E.
Sanchez-Palencia[38], J.L. Lions [28], J.B. Keller[25] and J-L Auriault et al
[5], just to name a few.

In the setting of a periodic pore microstructure, as the period goes to zero,
the convergence to the Darcy’s law of the Stokes system with no-slip boundary
condition posed on the boundary of the pore space was proved by L. Tartar
using the energy method [40]. G. Allaire implemented the two-scale convergence
method introduced by G. Nguetseng[33] to derive the Darcy’s law and show the
convergence [3, 4]. Prior to the proof of Darcy’s law in the ’80s, H.
Brinkman[14] studied the viscous force exerted by a flowing flow on a dense
swarm of particles by adding a diffusion term to the Darcy’s law so as to take
into account the transitional flow between boundaries. H. Brinkman’s method
was further studied in [39, 37, 30]. In the case of a porous material where
the solid region is much smaller than the fluid part, T. Levy [26] and E.
Sanchez-Palencia[38] proposed the same form of Darcy’s law but with a
different representation of the permeability tensor $K$. Later on, G.
Allaire[2] showed the continuity of the transition between the two forms of
Darcy’s laws by considering various ratios between the size of the solid
inclusion and the size of the separation. Moreover, instead of considering the
porous materials as a periodic structured material, A. Beliaev [9] considered
the porous materials as a random and stochastically homogeneous material and
deduced the same Darcy’s law. G. Allaire [1] generalized the homogenization to
handle the more realistic micro-geometries of the porous medium where both the
solid part and the fluid part are connected. Furthermore, in terms of the the
fluid-solid interface conditions, a slip boundary condition is considered by
G. Allaire in [1, 17]. In the case of the fluid flow through a porous medium
subject to a time-harmonic pressure gradient, the permeability depends on the
frequency and is referred to as the dynamic permeability. The theory of
dynamic permeability is established [8, 22, 5, 12] and further developed by
M.-Y. Ou [35].

The goal of this paper is to study how the permeability tensor derived from
the homogenization approach for porous materials [40, 41] depend on the
microstruture of the pore space. Details of this will be presented in Section
1.1.

The main tool we use will be the integral representation formula (IRF) for
composite materials. Composite materials are materials made from more than one
constituent materials with different physical or chemical properties. The
effective properties of composites, such as elasticity, conductivity and
permeability are of great interest in different application fields.
Homogenization theory for composite materials has been extensively studied in
[10, 31, 34]. Mathematically, for a two-component composite material,
themicrostructural information is carried into the analytical formulation of
the effective properties of the composite. Bergman pioneered the study of
analyticity of the effective dielectric constant [11], and in terms of
integral representation of effective material properties, a rigorous basis of
integral representation of the effective conductivity is established by K.
Golden and G. Papanicolaou [21], the effective elastic constants by Y. Kantor
and D. Bergman [23], the effective diffusivity in convection-enhanced
diffusion was derived by M. Avellaneda and A. Majda [6, 7]. Further
enlargements of the domain of analyticity of the IRF of elasticity tensor to
the case where one phase is a void or a hard inclusion is studied by Bruno and
Leo [15, 16].

Unlike the problem setup for calculating effective material properties such as
effective conductivity, elasticity, and diffusivity, where the physical
property of interest is well defined both in the micro-scale and the
macro(homogenized)-scale, the permeability of porous material is by definition
an effective property and hence it makes sense only in the macro-scale. To
overcome this difficulty, we consider a porous material as the limit case of a
two-fluid mixture.

Specifically, we will start with the two-fluid mixture problem studied in
[29], where the effective property is called the _self-permeability_. We will
derive the IRF for the self-permeability and show that the permeability for a
porous material is equal to the limit of the self-permeability when the
viscosity of one phase becomes infinite. Similar to the hard/soft inclusion
case studied in [15, 16], we will extend the domain of analyticity of the IRF
to $\infty$ and to $0$ by an iterative process. As a result, the IRF derived
here is valid for porous materials with a solid skeleton as well as for fluid-
bubble mixtures. Hence it provides a theoretical connection between the
permeability defined in [40] and the self-permeability for the bubbly fluid
studied in [29] and any mixture in between these two limiting cases.

The paper is organized as follows. The permeability of a porous material is
defined in Section 1.1. Section 2 starts with the definition of the _self-
permeability_ $K$of a two-fluid mixture and the corresponding cell problem,
followed by an analysis of the cell problem and the construction of the
solution in the vicinity of the two limiting cases of $z=\infty$ and $z=0$. In
Section 3, the IRF of $K$ is obtained by applying the theory of matrix-valued
Stieltjes functions. In this section, the spectral representation of $K$ is
also derived. The relationships between the moments of the measure in the IRF
and the geometry of the pore space are derived by comparing these two
representations. Section 4 presents the numerical solutions of the cell
problem of a special pore structure, which validate the theoretical results
given in Section 3.

Einstein summation convention is applied unless stated otherwise.

###  1.1 Definition of Permeability from Homogenization

Following the convention of homogenization, the space coordinates for the cell
problem in the open unit cell $Q=(0,1)^{n}$ for $n=2,3$, are denoted by
$\mathbf{y}=(y_{1},y_{2},y_{3})$. Let $\Omega$ be a smooth bounded open set
and $Q$ an open unit cube made of two open sets $Q_{1}$, $Q_{2}$ and the
interface $\Gamma=\mbox{cl}(Q_{1})\cap\mbox{cl}(Q_{2})$ with cl($A$) being the
closure of a set $A$. Moreover, $\widetilde{Q_{i}}$ denotes the $Q$-periodic
extension of $Q_{i}$, $i=1,2$. Following [1], we assume that (1) $Q_{1}$ and
$Q_{2}$ have strictly positive measures in cl($Q$). (2) The set
$\widetilde{Q_{i}}$ is open with $C^{1}$ boundary and is locally located on
one side of its boundary, $i=1,2$, and $\widetilde{Q_{1}}$ is connected. (3)
$Q_{1}$ is connected with a Lipschitz boundary. In addition, we consider the
case of inclusion, i.e. $Q_{2}\cap\partial Q=\emptyset$.

Consider $\epsilon>0$ much smaller than the size of $\Omega$ and $\epsilon
Q$-periodically extend $\epsilon Q_{1}$ in the entire space.
$\Omega_{\epsilon}$ denotes the intersection of $\Omega$ and this $\epsilon
Q$-periodically extended structure. In [40], the permeability is derived from
the Stokes equation in $\Omega_{\epsilon}$, which reads: find
$\mathbf{u}^{\epsilon}\in H_{0}^{1}(\Omega_{\epsilon})^{n}$ and
$p^{\epsilon}\in L^{2}(\Omega_{\epsilon})/\mathbb{R}$ such that

$\left\\{\begin{split}-\mu\triangle\mathbf{u}^{\epsilon}+\nabla
p^{\epsilon}&=\mathbf{f}\quad\text{ in }\Omega_{\epsilon}.\\\
\text{div}\mathbf{u}^{\epsilon}&=0\quad\text{ in }\Omega_{\epsilon}\\\
\end{split}\right.$ (1)

where $\mathbf{f}\in L^{2}(\Omega)$ is independent of $\epsilon$ and the
viscosity $\mu$ is a constant ($\mu$ is set to 1 in [40, 41]). See Figure 1
for an example of the unit cube. Note that the superscript $\epsilon$ is used
to signify that the solutions $\mathbf{u}^{\epsilon}$ and $p^{\epsilon}$
depend on $\epsilon$.

Figure 1: A sample illustration of a periodic cell.

To be able to prove the convergence of $(\mathbf{u}^{\epsilon},p^{\epsilon})$
as $\epsilon\rightarrow 0$, it is necessary to extend these solutions from
$\Omega_{\epsilon}$ to $\Omega$ so they are defined in the same spatial
domain. In [40, 41], $\mathbf{u}^{\epsilon}$ was extended by zero and
$p^{\epsilon}$ by a properly defined extension operator with their extensions
denoted by $\hat{\mathbf{u}^{\epsilon}}$ and $\hat{p^{\epsilon}}$,
respectively. As $\epsilon\rightarrow 0$,

$\frac{\hat{\mathbf{u}^{\epsilon}}}{\epsilon^{2}}\rightharpoonup\mathbf{U}\mbox{
weakly in
}L^{2}(\Omega)^{n},\,\text{div}\,\mathbf{U}=0,\mathbf{U}\cdot\mathbf{n}=0\mbox{
on }\Gamma,\mbox{ and }\hat{p^{\epsilon}}\rightarrow p\mbox{ in
}L^{2}(\Omega)/\mathbb{R}$

and the limit functions satisfy the following Darcy’s law [40]

$\mathbf{U}=\frac{\mbox{\boldmath$K$}^{(D)}}{\mu}(\mathbf{f}-\nabla p)$ (2)

where the permeability tensor $\mbox{\boldmath$K$}^{(D)}$ is defined as

$K_{ij}^{(D)}=\int_{Q_{1}}\mathbf{u}_{D}^{j}\cdot\mathbf{e}_{i}d\mathbf{y},\,i,j=1,\cdots
n$ (3)

with $\mathbf{e}_{i}$ denoting the unit vector in the $i$-th direction and
$\mathbf{u}_{D}^{j}$ the unique solution of the following boundary value
problem

$\left\\{\begin{split}\mu\Delta_{\mathbf{y}}\mathbf{u}_{D}^{j}-\nabla_{\mathbf{y}}p^{j}=-\mathbf{e}_{j}\quad&\text{in
}Q_{1}\\\ \text{div}_{\mathbf{y}}\mathbf{u}_{D}^{j}=0\quad&\text{in }Q_{1}\\\
\mathbf{u}_{D}^{j}=\mathbf{0}\quad&\text{on }\Gamma\end{split}\right.$ (4)

in the space $\mathring{H}(Q_{1}):=\left\\{\mathbf{v}:\mathbf{v}\in
H^{1}(Q_{1})^{n}\biggr{|}\;\text{div}_{\mathbf{y}}\mathbf{v}=0,\mathbf{v}|_{\Gamma}=\mathbf{0},Q\text{-periodic}\right\\}$.
Note that the subscript $i$ of $\mathbf{u}$ and $p$ signifies the solutions
corresponding to the force term $\mathbf{e}_{i}$.

Since $\mu$ is set to 1 in [40, 41], the permeability $K$ presented there is
related to $\mbox{\boldmath$K$}^{(D)}$ by
$\mbox{\boldmath$K$}^{(D)}={\mbox{\boldmath$K$}}/{\mu}$. For future analysis,
we will derive here the quadratic form representation of the permeability. We
start by observing that for incompressible fluid, we have

$\triangle\mathbf{u}=\mbox{div}(\nabla\mathbf{u}+\nabla^{T}\mathbf{u})$

Therefore, (3) can be expressed as

$K_{ij}^{(D)}=\int_{Q_{1}}\mathbf{u}^{j}_{{D}}\cdot\left(\nabla_{\mathbf{y}}p^{i}-{\mu}\Delta_{\mathbf{y}}\mathbf{u}^{i}_{{D}}\right)d\mathbf{y}=\int_{Q_{1}}{\mu}\nabla_{\mathbf{y}}\mathbf{u}^{j}_{{D}}:\left(\nabla_{\mathbf{y}}\mathbf{u}^{i}_{{D}}+\nabla_{\mathbf{y}}^{T}\mathbf{u}^{i}_{{D}}\right)d\mathbf{y}$
(5)

after applying Divergence theorem, periodicity of $\mathbf{u}$, $p$ and no-
slip conditions on $\Gamma$. Here we have used the Frobenius inner product of
matrices $\mathbf{A}:\mathbf{B}=\sum_{i,j=1}A_{ij}B_{ij}$. In terms of the
usual notion of the symmetric part and the antisymmetric part of vector field
$\nabla\mathbf{u}$

$\displaystyle
e(\mathbf{u}):=\frac{1}{2}\left(\nabla\mathbf{u}+\nabla^{T}\mathbf{u}\right),\,\,\tilde{e}(\mathbf{u}):=\frac{1}{2}\left(\nabla\mathbf{u}-\nabla^{T}\mathbf{u}\right){,}$
(6)

the right-hand side of equation (5) becomes
$\int_{Q_{1}}2\mu(e(\mathbf{u}_{D}^{j})+\tilde{e}(\mathbf{u}_{D}^{j})):e(\mathbf{u}_{D}^{i})d\mathbf{y}=\int_{Q_{1}}2\mu
e(\mathbf{u}_{D}^{j}):e(\mathbf{u}_{D}^{i})d\mathbf{y}$ because the Frobenius
product of a symmetric matrix and an antisymmetric matrix must be 0. Therefore
we have the quadratic form of permeability tensor $\mbox{\boldmath$K$}^{(D)}$

$K_{ij}^{(D)}=\int_{Q_{1}}2\mu
e(\mathbf{u}_{D}^{j}):e(\mathbf{u}_{D}^{i})d\mathbf{y}.$ (7)

##  2 Approximation of flow in porous medium by a two-phase Stokes flow

In this section, we consider the system for porous materials (4) as one of the
limiting cases of the two-fluid problem described below, which is the same as
the one studied in [29] with the exception that the fluid viscosity here can
be complex-valued. It is easy to check that the homogenization process in [29]
stays valid after making small modifications to accommodate the complex valued
viscosity described below.

Let $\Omega$, $Q$ and $\epsilon$ be the same as in Section 1.1. $Q_{2}$ is
still the inclusion in the periodic cell. Consider the ${\epsilon Q}$-periodic
extension of $\epsilon Q_{1}$ (${\epsilon}Q_{2}$) and denote by
$\Omega_{1{\epsilon}}$ ($\Omega_{2\epsilon}$) its intersection with $\Omega$.
We note that $\Omega_{1{\epsilon}}$ (region of the hosting fluid) is the same
as $\Omega_{\epsilon}$ in the previous section. Suppose $\Omega_{1{\epsilon}}$
is occupied by fluid with viscosity $\mu_{1}>0$ and $\Omega_{2}$ by fluid with
viscosity $z\mu_{1}$ with $z\in{\mathbb{C}}$. The interface
$\tilde{\Gamma}=\partial\Omega_{1{\epsilon}}\cap\partial\Omega_{2{\epsilon}}$
is such that $\Omega_{1{\epsilon}}\cup{\tilde{\Gamma}\ \cup\
}{\Omega_{2{\epsilon}}}=\Omega$. For the ease of notation, we define the
stress tensor $\boldsymbol{\tau}(\mathbf{u},\mbox{\boldmath$\mu$})$ of a fluid
with viscosity $\mu$, velocity field $\mathbf{u}$ and pressure field $p$ as

$\boldsymbol{\tau}(\mathbf{u},p,\mu)=2{\mu}e(\mathbf{u})-p\mathbf{I},\quad\mbox{$\mathbf{I}$
is the identity matrix.}$ (8)

Let $\chi_{i}$ be the characteristic function of $\Omega_{\epsilon i}$,
$i=1,2$, consider the viscosity function

${\xi}^{{\epsilon}}(\mathbf{x};z)=(\chi_{2}(\mathbf{x})z\mu_{1}+\chi_{1}(\mathbf{x})\mu_{1}),\quad
z\in{\mathbb{C}}.$ (9)

The two-fluid problem is given by the following Stokes system

$\left\\{\begin{split}\text{div}\left(2{{\xi}}^{\epsilon}(\mathbf{x};z)e(\mathbf{u}^{\epsilon})\right)-\nabla
p^{\epsilon}&=-\mathbf{f}\quad\text{ in }\Omega\backslash\tilde{\Gamma}\\\
\text{div}\mathbf{u}^{\epsilon}&=0\quad\text{ in }\Omega\\\
\mathbf{u}^{\epsilon}&=\mathbf{0}\quad\text{ on }\partial\Omega\\\
\llbracket\mathbf{u}^{\epsilon}\rrbracket=0,\,\mathbf{u}^{\epsilon}\cdot\mathbf{n}&=0\quad\text{on
}\tilde{\Gamma}\\\
\llbracket\mbox{\boldmath$\pi$}\rrbracket\cdot\mathbf{n}=\left(\llbracket\mbox{\boldmath$\pi$}\cdot\mathbf{n}\rrbracket\cdot\mathbf{n}\right)\mathbf{n}&\equiv\llbracket\mbox{\boldmath$\pi$}\cdot\mathbf{n}\rrbracket-\mathbf{n}\times\mathbf{n}\times\llbracket\mbox{\boldmath$\pi$}\cdot\mathbf{n}\rrbracket\quad\text{on
}\tilde{\Gamma}\end{split}\right.$ (10)

where
$\mbox{\boldmath$\pi$}=\mbox{\boldmath$\tau$}(\mathbf{u}^{\epsilon},p^{\epsilon},\boldsymbol{\xi}^{\epsilon})$,
$\mathbf{f}$ is a square integrable momentum source independent of $\epsilon$,
$\llbracket\cdot\rrbracket$ the jump across the interface $\tilde{\Gamma}$,
and $\mathbf{n}$ is the outward unit normal of $\partial\Omega_{2\epsilon}$.
The second jump condition in (10) means the traction can only jump in the
normal direction. Also note that the superscript $\epsilon$ is used to signify
that the solutions $\mathbf{u}^{\epsilon}$ and $p^{\epsilon}$ depend on
$\epsilon$.

It is shown in [29] that as $\epsilon\to 0$, $\mathbf{u}^{\epsilon}$ and the
properly normalized $p^{\epsilon}$, which is denoted by $\hat{p}^{\epsilon}$,
converge as follows

$\frac{{\mathbf{u}^{\epsilon}}}{\epsilon^{2}}\to\mathbf{u}^{0}\quad\text{weakly
in }L^{2}(\Omega)^{n},\;\hat{p^{\epsilon}}\to P\quad\text{strongly in
}L^{2}(\Omega)/\mathbb{R}$

where $\mathbf{u}^{0}$ and $P$ satisfy the homogenized system:

$\left\\{\begin{split}\mathbf{u}^{0}&=-\mbox{\boldmath$K$}(\nabla
P-\mathbf{f})\quad\text{in }\Omega\\\ \text{div
}\mathbf{u}^{0}&=0\quad\text{in }\Omega\\\ \end{split}\right.$ (11)

where the components of $K$, which is referred to as the _self-permeability_
in [29], is defined as

$K_{ij}(z)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\int_{Q}{\mathbf{u}^{j}}\cdot\mathbf{e}_{i}d\mathbf{y},\,{i,j=1,\dots,n}$
(12)

with $\mathbf{u}^{{i}}$ being the unique solution to the cell problem posed in
the function space $H(Q)$, which is defined in (14),

$\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2{\mbox{\boldmath$\mu$}}(\mathbf{y};z)e(\mathbf{u}^{i})-p^{i}\mathbf{I}\right)+\mathbf{e}_{i}&=\mathbf{0}\quad\text{in
}Q_{1}\cup Q_{2}\\\
\llbracket\mbox{\boldmath$\pi$}\rrbracket\cdot\mathbf{n}&=\left(\llbracket\mbox{\boldmath$\pi$}\cdot\mathbf{n}\rrbracket\cdot\mathbf{n}\right)\mathbf{n}\text{
on }\Gamma\end{split}\right.$ (13)

where
$\mbox{\boldmath$\mu$}(\mathbf{y};z)=\mu_{1}\chi_{1}(\mathbf{y})+z\mu_{1}\chi_{2}(\mathbf{y})$
with $\chi_{m}$ being the characteristic functions of $Q_{m}$, $m=1,2,$ and
$\mbox{\boldmath$\pi$}=\mbox{\boldmath$\tau$}(\mathbf{u}^{k},p^{k},\mbox{\boldmath$\mu$})$,
cf. (8). Note that the superscript $i$ is used to signify that
$\mathbf{u}^{i}$ and $p^{i}$ are solutions to the cell problem (13) with the
force term $-\mathbf{e}_{i}$, $i=1,\dots,n$.

###  2.1 Function Spaces

Let $\mathcal{R}(Q_{2})$ denote the space of rigid body displacements in
$Q_{2}$, i.e. $\mathbf{u}=\mathbf{A}\mathbf{y}+\mathbf{b}$ with constant skew-
symmetric matrix $A$ and constant vector $\mathbf{b}$ in $Q_{2}$. We start
with the space of admissible functions for the velocity

$\displaystyle H(Q)$ $\displaystyle:=\left\\{\mathbf{v}:\mathbf{v}\in
H^{1}(Q_{1}\cup
Q_{2})^{n}\biggr{|}\;\text{div}_{\mathbf{y}}\mathbf{v}=0,\;\mathbf{v}\cdot\mathbf{n}=0\text{
in }H^{-\frac{1}{2}}(\Gamma),\right.$
$\displaystyle\qquad\left.{}\llbracket\mathbf{v}\rrbracket_{\Gamma}=\mathbf{0},\;(\mathbf{v},\mathbf{\eta})_{H^{1}(Q_{2})}=0,\forall\mathbf{\eta}\in\mathcal{R}(Q_{2}),\,\mathbf{v}\text{
is }Q\text{- periodic}\right\\}$ (14)

where $\mathbf{n}$ is the outward unit normal of $\partial Q_{2}$. $H(Q)$ is
endowed with the inner product

$(\mathbf{u},\mathbf{v})_{Q}=\int_{Q}2\mu_{1}e(\mathbf{u}):\overline{e(\mathbf{v})}d\mathbf{y}.$
(15)

The induced norm is denoted by
$\left\lVert\mathbf{u}\right\rVert_{Q}^{2}:=(\mathbf{u},\mathbf{u})_{Q}$. Note
that we have $H(Q)\cap\mathcal{R}(Q)=\\{\mathbf{0}\\}$ because $\mathbf{A}=0$
due to the $Q$-periodicity and $\mathbf{u}\cdot\mathbf{n}=0$ implies
$\mathbf{b}=\mathbf{0}$. We observe that if $\mathbf{u}\in H(Q)$ then
$\mathbf{u}\in H^{1}(Q)^{n}$ by the following argument. Obviously,
$\mathbf{u}\in L^{2}(Q)^{n}$. To prove $\frac{\partial u_{i}}{\partial
y_{j}}\in L^{2}(Q)$ for $i,j=1,2,3$, let $\phi$ be any $C^{\infty}$ test
function compactly supported in $Q$ and $h$ be the $i$-th component $u_{i}$
for any $i$. Then

$\int_{Q}h\nabla\phi
d\mathbf{y}=-\left(\int_{Q_{1}\cap\text{Supp}(\phi)}\phi\nabla
hd\mathbf{y}+\int_{Q_{2}\cap\text{Supp}(\phi)}\phi\nabla hd\mathbf{y}\right)$

here we used $\llbracket h\rrbracket=\mathbf{0}$. Now we can define a
candidate function $\mathbf{g}$ such that

$\begin{split}\mathbf{g}|_{Q_{i}}&:=\nabla h|_{Q_{i}},\,i=1,2\end{split}$ (16)

then clearly $\mathbf{g}\in L^{2}(Q)^{n}$ and
$\left<h,\nabla\phi\right>=-\left<g,\phi\right>$, where $\left<\cdot\right>$
denotes the usual $L^{2}$ inner product. Therefore ${h}\in H^{1}(Q)$ and hence
$u_{i}\in H^{1}(Q),\ i=1,\dots,n$.

Next, we show that $\|\cdot\|_{Q}$ is equivalent to the usual $H^{1}$ norm,
i.e., there exist constants $B_{1}$ and $B_{2}$ such that

$B_{1}\left\lVert\mathbf{u}\right\rVert_{H^{1}(Q)}\leq\left\lVert\mathbf{u}\right\rVert_{Q}\leq
B_{2}\left\lVert\mathbf{u}\right\rVert_{H^{1}(Q)}$ (17)

Because $H^{1}(Q)\cap\mathcal{R}(Q)=\\{0\\}$, by Theorem 2.5 in [34], there
exists a Korn’s constant $C_{1}$ such that

$C_{1}\left\lVert\mathbf{u}\right\rVert_{H^{1}(Q)}\leq\frac{1}{\sqrt{2\mu_{1}}}\left\lVert\mathbf{u}\right\rVert_{Q}$
(18)

where $C_{1}$ depends only on $Q$. Therefore, we can take
$B_{1}=\sqrt{2\mu_{1}}C_{1}$. To emphasize the dependence on $Q$, we will
write it as $B_{1}(Q)$. On the other hand, according to the orthogonal
decomposition that $\nabla\mathbf{u}=e(\mathbf{u})+\tilde{e}(\mathbf{u})$, see
(6),

$\left\lVert\mathbf{u}\right\rVert_{H^{1}(Q)}^{2}\geq\left\lVert\nabla\mathbf{u}\right\rVert_{L^{2}(Q)}^{2}=\left\lVert
e(\mathbf{u})\right\rVert_{L^{2}(Q)}^{2}+\left\lVert\tilde{e}(\mathbf{u})\right\rVert^{2}\geq\left\lVert
e(\mathbf{u})\right\rVert_{L^{2}(Q)}^{2}=\frac{1}{2\mu_{1}}\|\mathbf{u}\|_{Q}^{2}$

therefore $B_{2}=\sqrt{2\mu_{1}}$. The reason for introducing the $H(Q)$-norm
is that the self-permeability in (12) can be represented in terms of the inner
product. More specifically, using (12), (13) and the fact that
$\overline{\mathbf{e}_{i}}=\mathbf{e}_{i}$, by a calculation similar to (5)
and taking into account the interface condition $\mathbf{u}\cdot\mathbf{n}=0$
and the jump conditions in (13), (12) can be expressed in the following form

$\begin{split}K_{ij}(z)=\int_{Q}2{\mu}(\mathbf{y};\overline{z})\overline{e(\mathbf{u}^{{i}}(z))}:{e(\mathbf{u}^{{j}}(z))}d\mathbf{y}\end{split}$
(19)

and its conjugate transpose
$\mbox{\boldmath$K$}^{*}:=\overline{\mbox{\boldmath$K$}^{T}}$ is

${(K^{*})_{ij}}(z)=\int_{Q}2{\mu}(\mathbf{y};{z}){e(\mathbf{u}^{{j}}(z))}:\overline{e(\mathbf{u}^{{i}}(z))}d\mathbf{y}$
(20)

###  2.2 Weak solution of the Cell Problem (13)

The weak formulation of the cell problem (13) is

$\int_{Q_{1}\cup
Q_{2}}2\mu(\mathbf{y};z)e(\mathbf{u}^{k}):\overline{e(\mathbf{v})}d\mathbf{y}=\int_{Q_{1}\cup
Q_{2}}\mathbf{e}_{k}\cdot\bar{\mathbf{v}}d\mathbf{y},\qquad\forall\mathbf{v}\in
H(Q)$ (21)

From this, we see that the solutions satisfy the following symmetry

$\mathbf{u}^{k}(\mathbf{y};\overline{z})=\overline{\mathbf{u}^{k}(\mathbf{y};z)}$

Define the sesquilinear form on $H(Q)$

$a(\mathbf{u},\mathbf{v})=\int_{Q_{1}\cup
Q_{2}}2\mu(\mathbf{y};z)e(\mathbf{u}^{k}):\overline{e(\mathbf{v})}d\mathbf{y}$
(22)

It is clear that $a(\mathbf{u},\mathbf{v})$ is bounded in $H(Q)$. To check the
coercivity, assume $\mathbf{u}^{k}\neq 0$ and define the parameter

$\lambda:=\frac{\int_{Q}2\mu_{1}\chi_{2}e(\mathbf{u}^{k}):\overline{e(\mathbf{u}^{k})}d\mathbf{y}}{\int_{Q}2\mu_{1}e(\mathbf{u}^{k}):\overline{e(\mathbf{u}^{k})}d\mathbf{y}}$
(23)

then $0\leq\lambda\leq 1$. We note that

$\displaystyle\frac{a(\mathbf{u}^{k},\mathbf{u}^{k})}{\int_{Q}2\mu_{1}e(\mathbf{u}^{k}):\overline{e(\mathbf{u}^{k})}d\mathbf{y}}=\lambda
z+(1-\lambda)\cdot 1$ (24)

and hence as long as 0 is not on the line segment joining $z$ and 1, there
exist $\alpha(z):=\min_{0\leq\lambda\leq 1}|\lambda z+1-\lambda|>0$ such that

$\left|a(\mathbf{u}^{k},\mathbf{u}^{k})\right|\geq\alpha(z)\int_{Q}2\mu_{1}e(\mathbf{u}^{k}):\overline{e(\mathbf{u}^{k})}d\mathbf{y}=\alpha\|\mathbf{u}^{k}\|_{Q}^{2}$
(25)

Therefore for $z\in{\mathbb{C}}\backslash\left\\{\Re{z}\leq 0\right\\}$, by
the Lax-Milgram Lemma [13, Chapter 2], there exists a unique weak solution
$\mathbf{u}^{k}\in H(Q)$ to the cell problem (13) and with the solution
$\mathbf{u}^{k}$, we can construct $p^{k}\in L^{2}(Q)/\mathbb{C}$.

Since $\alpha(z)$ is a continuous function in $z$, the coercivity of the
sesquilinear form can be applied to conclude that $\mathbf{u}^{k}$ is analytic
in $z$ and its $m$-th derivative, $m\geq 1$, satisfies the following recursive
equation

$\int_{Q_{1}\cup
Q_{2}}2\mu(\mathbf{y};z)e\left(\frac{d^{m}\mathbf{u}^{k}}{dz^{m}}\right):\overline{e(\mathbf{v})}d\mathbf{y}=-\int_{Q_{2}}2{m}\mu_{1}e\left(\frac{d^{m-1}\mathbf{u}^{k}}{dz^{m-1}}\right):e(\overline{\mathbf{v}})d\mathbf{y},\,\forall\mathbf{v}\in
H(Q)$ (26)

As a result, $\mbox{\boldmath$K$}(z)$ is also analytic for
$z\in{\mathbb{C}}\backslash\left\\{\Re{z}\leq 0\right\\}$. To relate the two-
fluid problem with $\mbox{\boldmath$K$}^{(D)}$, we adapt the method used in
[16] to study the behavior of $\mbox{\boldmath$K$}(z)$ near $z=\infty$ in the
following section.

###  2.3  Analyticity of the Solution for large $|z|$

Let $w:=\frac{1}{z}$ and consider $Q$-periodic solution in the series form
near $w=0$

$\mathbf{u}_{\infty}(\mathbf{y};\mathbf{e},w):=\sum_{k=0}^{\infty}\mathbf{u}_{k}(\mathbf{y};\mathbf{e})w^{k}\mbox{
and
}p_{\infty}(\mathbf{y};\mathbf{e},w):=\sum_{k=0}^{\infty}p_{k}(\mathbf{y};\mathbf{e})w^{k}$
(27)

where the $\mathbf{e}$ is an arbitrary constant unit vector. To set up the
notation, we denote the restrictions of $\mathbf{u}_{k}$, $p_{k}$ in $Q_{2}$
(inclusion) and $Q_{1}$ as $\mathbf{u}^{in}_{k}$, $p^{in}_{k}$ and
$\mathbf{u}^{out}_{k}$, $p^{out}_{k}$ respectively and define

$\displaystyle\mathbf{u}_{\infty}^{in}(\mathbf{y};\mathbf{e},w):=\sum_{k=0}^{\infty}\mathbf{u}^{in}_{k}(\mathbf{y},\mathbf{e})w^{k},\qquad\mathbf{u}_{\infty}^{out}(\mathbf{y};\mathbf{e},w):=\sum_{k=0}^{\infty}\mathbf{u}^{out}_{k}(\mathbf{y},\mathbf{e})w^{k}$
(28)

By substituting (27) into (13) with the viscosity defined in (9), taking into
account the additional two interface conditions $\mathbf{u}\cdot\mathbf{n}=0$
and $\llbracket\mathbf{u}\rrbracket=0$, followed by equating terms of the same
order with respect to $w$, we arrive in the following equations in $Q_{1}$:

$\displaystyle
O(w^{0}):\qquad\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{0}^{out})-p_{0}^{out}\mathbf{I}\right)=-\mathbf{e}$
(29) $\displaystyle
O(w^{k}):\qquad\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)=\mathbf{0}\mbox{
for }k\geq 1$ (30)

and in $Q_{2}$:

$\displaystyle
O(w^{-1}):\qquad\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{0}^{in})\right)=\mathbf{0}$
(31) $\displaystyle
O(w^{0}):\qquad\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{1}^{in})-p_{0}^{in}\mathbf{I}\right)=-\mathbf{e}$
(32) $\displaystyle
O(w^{k}):\qquad\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{k+1}^{in})-p_{k}^{in}\mathbf{I}\right)=\mathbf{0}\mbox{
for }k\geq 1$ (33)

and the following interface conditions on $\Gamma$

$\displaystyle O(w^{-1}):$ $\displaystyle\
2\mu_{1}(e(\mathbf{u}_{0}^{in})\cdot\mathbf{n})|_{\Gamma}=C(\mathbf{y})\mathbf{n}\text{
for some function }C(\mathbf{y})$ (34) $\displaystyle O(w^{k}),\,k\geq 0:$
$\displaystyle\
\left(\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}_{k+1}^{in})-p_{k}^{in}\mathbf{I}\right)\right)\mathbf{n}$
(36)
$\displaystyle=\left\\{\left[\left(\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}_{k}^{in})-p_{k}^{in}\mathbf{I}\right)\right)\mathbf{n}\right]\cdot\mathbf{n}\right\\}\mathbf{n},$
$\displaystyle\mathbf{u}^{in}_{k}\cdot\mathbf{n}=\mathbf{u}^{out}_{k}\cdot\mathbf{n}=0\mbox{
and }\mathbf{u}^{in}_{k}=\mathbf{u}^{out}_{k}$

We introduce the following spaces, $i=1,2$

$\displaystyle H(Q_{1})=\left\\{\mathbf{v}:\mathbf{v}\in
H^{1}(Q_{1})^{n}\biggr{|}\;\text{div}_{\mathbf{y}}\mathbf{v}=0,\mathbf{v}\cdot\mathbf{n}=0\text{
on }\Gamma,Q\text{-periodic}\right\\}$ $\displaystyle
H(Q_{2})=\left\\{\mathbf{v}:\mathbf{v}\in
H^{1}(Q_{2})^{n}\biggr{|}\;\text{div}_{\mathbf{y}}\mathbf{v}=0,\mathbf{v}\cdot\mathbf{n}=0\text{
on }\Gamma,\,(\mathbf{v},\mathcal{R}(Q_{2}))_{H^{1}(Q_{2})}=0\right.$
$\displaystyle\left.\qquad\qquad,Q\text{-periodic}\right\\}$
$\displaystyle\mathring{H}(Q_{i})=\left\\{\mathbf{v}:\mathbf{v}\in
H^{1}(Q_{i})^{n}\biggr{|}\;\text{div}_{\mathbf{y}}\mathbf{v}=0,\mathbf{v}|_{\Gamma}=\mathbf{0},Q\text{-periodic}\right\\}\subset
H(Q_{i}),\,$ $\displaystyle L(Q_{i})/\mathbb{C}=\left\\{p:p\in
L^{2}(Q_{i}),\int_{Q_{i}}p(\mathbf{y})d\mathbf{y}=0,Q\text{-periodic},\right\\}$

Note that $H(Q_{1})\cap\mathcal{R}(Q_{1})=\\{\mathbf{0}\\}$ because $\partial
Q\subset\partial Q_{1}$. For $H(Q_{2})$, the boundary condition
$\mathbf{u}\cdot\mathbf{n}=0$ implies
$H(Q_{2})\cap\mathcal{R}(Q_{2})=\\{\mathbf{0}\\}$ because of the extra
condition $(\mathbf{v},\mathcal{R}(Q_{2}))_{H^{1}(Q_{2})}=0$[34]. Therefore,
$H(Q_{i})$ and $\mathring{H}(Q_{i})$ are equipped with inner product
$(\mathbf{u},\mathbf{v})_{Q_{i}}=\int_{Q_{i}}2\mu_{1}e(\mathbf{u}):\overline{e(\mathbf{v})}d\mathbf{y}$
and Korn’s inequalities are valid in $H(Q_{i})$, $i=1,2$.

######  Lemma 2.1.

Let $Q_{2}$ be a connected, open bounded set such that $\partial
Q_{2}\cap\partial Q=\emptyset$ and $\partial Q_{2}$ is in
$\mathcal{C}^{k,\sigma}$ , $k,\sigma\geq 0$, $k+\sigma\geq 2$. For any vector
field $\mathbf{u}^{in}\in H(Q_{2})$, there exists a unique weak solution
$\mathbf{u}^{out}(\mathbf{y};\mathbf{f}^{out})\in H(Q_{1})$ that satisfies the
following system

$\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}^{out})-p^{out}\mathbf{I}\right)&=\mathbf{f}^{out}\quad\text{
in }Q_{1}\\\ \mathbf{u}^{out}&=\mathbf{u}^{in}\quad\text{ on
}\Gamma\end{split}\right.$ (37)

where in our context, $\mathbf{f}^{out}=\mathbf{0}\text{ or
}\mathbf{f}^{out}=-\mathbf{e}$, a constant unit vector. Moreover,

$\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}\leq\frac{1}{B_{1}(Q)}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}+2E_{1}\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}.$
(38)

where the positive constants $B_{1}(Q_{{1}})$ is defined in (17) and
$E_{1}\geq 1$ depends only on $Q_{1}$ and $Q_{2}$.

###### Proof.

To handle the inhomogeneous boundary condition, we proceed as follows. By [24,
Corollary 3.2], there exists a bounded, divergence free extension
$T\left(\mathbf{u}^{in}\right)$ of $\mathbf{u}^{in}$ to a small neighborhood
$O$ of $Q_{2}$ and vanishes at $\partial O\subset Q_{1}$ such that

$\left\lVert T\left(\mathbf{u}^{in}\right)\right\rVert_{Q}\leq
E_{1}\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}$ (39)

where $E_{1}\geq 1$ depends only on $Q_{1}$ and $Q_{2}$. Furthermore, the
extension $T\left(\mathbf{u}^{in}\right)$ on $O$ can be extended periodically
to $\mathbb{R}^{n}$ [18] since $T\left(\mathbf{u}^{in}\right)$ vanishes on
$\partial O$ and hence on $\partial Q$. We denote the restriction
$T(\mathbf{u}^{in})|_{Q_{1}}$ as $\tilde{\mathbf{u}}^{out}\in H(Q_{1})$ and
$\mathring{\mathbf{u}}^{out}\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\mathbf{u}^{out}-\tilde{\mathbf{u}}^{out}\in\mathring{H}(Q_{1})$
and (37) becomes

$\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathring{\mathbf{u}}^{out})-p^{out}\mathbf{I}\right)=\mathbf{f}^{out}-\mu_{1}\Delta\tilde{\mathbf{u}}^{out}\quad\text{
in }Q_{1}$ (40)

Consider the variational formulation: Find
$\mathring{\mathbf{u}}^{out}\in\mathring{H}(Q_{1})$ such that
$\forall\Phi\in\mathring{H}(Q_{1})$,

$\int_{Q_{1}}2\mu_{1}e(\mathring{\mathbf{u}}^{out}):\overline{e(\Phi)}d\mathbf{y}=-\int_{Q_{1}}\mathbf{f}^{out}\cdot\overline{\Phi}d\mathbf{y}-\int_{Q_{1}}2\mu_{1}e(\tilde{\mathbf{u}}^{out}):\overline{e(\Phi)}d\mathbf{y},\quad$
(41)

The right hand side of (41) can be bounded as follows

$\displaystyle\left|\int_{Q_{1}}\mathbf{f}^{out}\cdot\overline{\Phi}d\mathbf{y}+\int_{Q_{1}}2\mu_{1}e(\tilde{\mathbf{u}}^{out}):\overline{e(\Phi)}d\mathbf{y}\right|\leq$
$\displaystyle\left(\frac{\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}}{B_{1}(Q_{{1}})}+\left\lVert\tilde{\mathbf{u}}^{out}\right\rVert_{Q_{1}}\right)\left\lVert\Phi\right\rVert_{Q_{1}}$

The sesquilinear form
$\int_{Q_{1}}2\mu_{1}e(\mathring{\mathbf{u}}^{out}):\overline{e(\Phi)}d\mathbf{y}$
is clearly bounded and coercive with constant 1. Hence by the Lax-Milgram
Lemmat there exists a unique weak solution
$\mathring{\mathbf{u}}^{out}\in\mathring{H}(Q_{1})$ to (40) such that
$\left\lVert\mathring{\mathbf{u}}^{out}\right\rVert_{Q_{1}}\leq(\frac{1}{B_{1}(Q_{{1}})}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}+\left\lVert\tilde{\mathbf{u}}^{out}\right\rVert_{Q_{1}})$.
In terms of $\mathring{\mathbf{u}}^{out}$, $\mathbf{u}^{out}$ can be expressed
as $\mathbf{u}^{out}=\tilde{\mathbf{u}}^{out}+\mathring{\mathbf{u}}^{out}$ and
satisfies the estimate

$\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}\leq\left\lVert\mathring{\mathbf{u}}^{out}\right\rVert_{Q_{1}}+\left\lVert\tilde{\mathbf{u}}^{out}\right\rVert_{Q_{1}}\leq\frac{1}{B_{1}(Q_{{1}})}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}+2E_{1}\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}$
(42)

To show the solution $\mathbf{u}^{out}$ is unique, suppose
$\mathbf{u}^{out}_{1}$ and $\mathbf{u}^{out}_{2}$ both solve (37) then the
difference
${\mathbf{w}}^{\text{diff}}=\mathbf{u}^{out}_{1}-\mathbf{u}^{out}_{2}\in\mathring{H}(Q_{1})$
must satisfy

$\int_{Q_{1}}2\mu_{1}e({\mathbf{w}}^{\text{diff}}):\overline{e(\Phi)}d\mathbf{y}=0,\qquad\forall\Phi\in\mathring{H}(Q_{1})$

Hence ${\mathbf{w}}^{\text{diff}}=0$ in $Q_{1}$ because
${\mathbf{w}}\in\mathring{H}(Q_{1})$. We note the two special cases:

$\displaystyle\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}\leq
2E_{1}\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}\mbox{ for
}\mathbf{f}^{out}=\mathbf{0}$ (43)
$\displaystyle\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}\leq\frac{1}{B_{1}(Q_{{1}})}\sqrt{|Q_{1}|}+2E_{1}\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}\mbox{
for }\mathbf{f}^{out}=-\mathbf{e}_{{j}},\,{j=1,\dots,n,}$ (44)

where $|Q_{1}|$ is the volume of $Q_{1}$.∎

######  Lemma 2.2.

Let $Q_{2}$ satisfy the same assumptions as those in Lemma 2.1. For any pair
of $\left(\mathbf{u}^{out},p^{out}\right)\in H(Q_{1})\times
L^{2}(Q_{1})/\mathbb{C}$ that satisfies (37), there exists a unique vector
$\mathbf{u}^{in}(\mathbf{y};\mathbf{f}^{in})\in H(Q_{2})$ that satisfies the
Stokes equation with continuity of tangential traction on $\Gamma$

$\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}^{in})-p^{in}\mathbf{I}\right)&=\mathbf{f}^{in}\quad\text{in
}Q_{2},\\\
\mathbf{n}\times\mathbf{n}\times\left[\left(\left(2\mu_{1}e(\mathbf{u}^{out})-p^{out}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}^{in})-p^{in}\mathbf{I}\right)\right)\cdot\mathbf{n}\right]&=\mathbf{0}\quad\text{on
}\Gamma,\end{split}\right.$ (45)

where in our context, $\mathbf{f}^{in}=\mathbf{0}\text{ or
}\mathbf{f}^{in}=-\mathbf{e}$. Moreover,

$\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}\leq\frac{1}{B_{1}(Q_{2})}\left\lVert\mathbf{f}^{in}\right\rVert_{L^{2}(Q_{2})}+\frac{E_{1}(Q_{1})}{B_{1}(Q_{1})}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}+E_{1}(Q_{1})\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}$

and $E_{1}(Q_{1})$, $B_{1}(Q_{1})$ and $B_{1}(Q_{2})$ depend only on $Q$ and
$\Gamma$.

###### Proof.

Take ${\Phi}\in H(Q_{2})$, the variation formulation for the PDE is

$\int_{Q_{2}}\left(\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}^{in})-p^{in}\mathbf{I}\right)\right)\cdot\bar{\Phi}d\mathbf{y}=\int_{Q_{2}}\mathbf{f}^{in}\cdot\bar{\Phi}d\mathbf{y}.$

For the ease of notation, let
$\mbox{\boldmath$\pi$}^{in}:=\mbox{\boldmath$\tau$}(\mathbf{u}^{in},\mu_{1},p^{in})$
and
$\mbox{\boldmath$\pi$}^{out}:=\mbox{\boldmath$\tau$}(\mathbf{u}^{out},\mu_{1},p^{out})$
where the stress function $\tau$ is defined in (8). Applying integration by
parts on the left hand side, followed by an application of the divergence
theorem leads to

$-\int_{\Gamma}({\mbox{\boldmath$\pi$}^{in}\cdot\bar{\Phi}})\cdot\mathbf{n}dS-\int_{Q_{2}}2\mu_{1}e(\mathbf{u}^{in}):\overline{e(\Phi)}d\mathbf{y}=\int_{Q_{2}}\mathbf{f}^{in}\cdot\bar{\Phi}d\mathbf{y}$
(46)

Let $\mathbf{t}$ denote the unit vector in the tangent plane such that
$\mathbf{t}\cdot\mathbf{n}=0$. The conditions on $\Gamma$ in (45) imply

$\left\\{\begin{split}&\Phi\cdot\mathbf{n}=0\Rightarrow\Phi=d(\mathbf{y})\mathbf{t}\mbox{
for some function }d(\mathbf{y}),\\\
&\mathbf{n}\times\mathbf{n}\times\left[\left(\mbox{\boldmath$\pi$}^{out}-\mbox{\boldmath$\pi$}^{in}\right)\cdot\mathbf{n}\right]=\mathbf{0}\Rightarrow\left(\mbox{\boldmath$\pi$}^{out}-\mbox{\boldmath$\pi$}^{in}\right)\cdot\mathbf{n}=C(\mathbf{y})\mathbf{n}\mbox{
for some function }C(\mathbf{y})\end{split}\right.$

With these observations and the fact that $\mbox{\boldmath$\pi$}^{in}$ is
symmetric, the first term in (46) can be expressed as

$\displaystyle-\int_{\Gamma}{\bar{\Phi}\cdot\mbox{\boldmath$\pi$}^{in}}\cdot\mathbf{n}\,dS$
$\displaystyle=\int_{\Gamma}\overline{d(\mathbf{y})}\mathbf{t}\cdot\left[\left(\mbox{\boldmath$\pi$}^{out}-\mbox{\boldmath$\pi$}^{in}\right)\cdot\mathbf{n}\right]-\int_{\Gamma}\overline{d(\mathbf{y})}\mathbf{t}\cdot(\mbox{\boldmath$\pi$}^{out}\cdot\mathbf{n})\,dS$
$\displaystyle=$
$\displaystyle-\int_{\Gamma}\bar{\Phi}\cdot\mbox{\boldmath$\pi$}^{out}\cdot\mathbf{n}\,dS$

and hence the variational form (46) becomes for all $\Phi\in H(Q_{2})$

$-\int_{Q_{2}}2\mu_{1}e(\mathbf{u}^{in}):\overline{e(\Phi)}d\mathbf{y}=\int_{Q_{2}}\mathbf{f}^{in}\cdot\bar{\Phi}d\mathbf{y}+\int_{\Gamma}\bar{\Phi}\cdot\mbox{\boldmath$\pi$}^{out}\cdot\mathbf{n}dS$
(47)

To bound the right hand side of (47), we first extend $\Phi\in H(Q_{2})$ by
the operator $T$ described in (39)

$\left\lVert T(\Phi)\right\rVert_{Q}\leq
E_{1}\left\lVert\Phi\right\rVert_{Q_{2}}$ (48)

$T(\Phi)$ rapidly decays to zero in a small neighborhood of $Q_{2}$ and stays
$0$ for the rest of $Q_{1}$. The restriction of $T(\Phi)$ in $Q_{1}$, denoted
by $\Phi^{out}$, has the following estimate

$\left\lVert\Phi^{out}\right\rVert_{Q_{1}}\leq\left\lVert
T(\Phi)\right\rVert_{Q}\leq E_{1}\left\lVert\Phi\right\rVert_{Q_{2}}$ (49)

Hence

$\begin{split}&\left|\int_{\Gamma}\bar{\Phi}\cdot\mbox{\boldmath$\pi$}^{out}\cdot\mathbf{n}dS\right|=\left|\int_{\Gamma}\bar{\Phi}^{out}\cdot\mbox{\boldmath$\pi$}^{out}\cdot\mathbf{n}dS\right|\\\
=&\left|\int_{Q_{1}}\left[\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}^{out})-p^{out}\mathbf{I}\right)\right]\cdot\overline{\Phi^{out}}d\mathbf{y}+\int_{Q_{1}}2\mu_{1}e(\mathbf{u}^{out}):\overline{e(\Phi^{out})}d\mathbf{y}\right|\\\
\leq&\left|\int_{Q_{1}}\mathbf{f}^{out}\cdot\overline{\Phi^{out}}d\mathbf{y}\right|+\left|\int_{Q_{1}}2\mu_{1}e(\mathbf{u}^{out}):\overline{e(\Phi^{out})}d\mathbf{y}\right|\\\
\leq&\frac{E_{1}}{B_{1}(Q_{{1}})}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}\left\lVert\Phi\right\rVert_{Q_{2}}+E_{1}\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}\left\lVert\Phi\right\rVert_{Q_{2}}\end{split}$

Therefore the right hand side of (47) is bounded by

$\left\lVert\Phi\right\rVert_{Q_{2}}\left(\frac{E_{1}\left\lVert\mathbf{f}^{in}\right\rVert_{L^{2}(Q_{2})}}{B_{1}(Q_{{2}})}+\frac{E_{1}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}}{B_{1}(Q_{{1}})}+E_{1}\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}\right)$

Finally, by the Lax-Milgram Lemma a unique solution $\mathbf{u}^{in}$ exists
such that

$\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}\leq\frac{E_{1}}{B_{1}(Q_{{2}})}\left\lVert\mathbf{f}^{in}\right\rVert_{L^{2}(Q_{2})}+\frac{E_{1}}{B_{1}(Q_{{1}})}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}+E_{1}\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}$
(50)

∎

The construction of the solution for $z$ with large magnitude will be carried
out using the following steps.

  1. 1.

$O(w^{-1})$: Consider the system of (31) and (34) for
$\mathbf{u}_{0}^{in}(\mathbf{y};\mathbf{e})\in H(Q_{2})$.

$\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{0}^{in})\right)&=\mathbf{0}\quad\text{in
}Q_{2}\\\
2\mu_{1}e(\mathbf{u}_{0}^{in})\cdot\mathbf{n}&=C(\mathbf{y})\mathbf{n}\quad\text{on
}\Gamma\end{split}\right.$ (51)

The variational formulation is

$-\int_{\Gamma}\bar{\mathbf{v}}\cdot
2\mu_{1}e(\mathbf{u}_{0}^{in})\cdot\mathbf{n}-\int_{Q_{2}}2\mu_{1}e(\mathbf{u}_{0}^{in}):\overline{e(\mathbf{v})}=0\quad\forall\mathbf{v}\in
H(Q_{2})$ (52)

The first term vanishes because of the boundary conditions. Hence

$\mathbf{u}_{0}^{in}(\mathbf{y};\mathbf{e})=\mathbf{0}\mbox{ in }Q_{2}\mbox{
{because $H(Q_{2})\perp\mathcal{R}(Q_{2})$.}}$

  2. 2.

$O(w^{0})$ in $Q_{1}$: Solve the system of (36) and (29) for
$\mathbf{u}_{0}^{out}(\mathbf{y};\mathbf{e})\in H(Q_{1})$.

$\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{0}^{out})-p_{0}^{out}\mathbf{I}\right)&=-\mathbf{e}\quad\text{in
}Q_{1}\\\ \mathbf{u}_{0}^{out}=\mathbf{u}_{0}^{in}&=\mathbf{0}\quad\text{on
}\Gamma\end{split}\right.$ (53)

An application of Lemma 2.1 and (44) leads to the following result

$\left\lVert\mathbf{u}_{0}^{out}\right\rVert_{Q_{1}}\leq\frac{\sqrt{|Q_{1}|}}{B_{1}(Q_{{1}})}+2E_{1}\left\lVert\mathbf{u}^{in}_{0}\right\rVert_{Q_{2}}=\frac{\sqrt{|Q_{1}|}}{B_{1}(Q_{{1}})}$
(54)

  3. 3.

$O(w^{0})$ in $Q_{2}$: Consider the system of (32) and (36) for
$\mathbf{u}_{1}^{in}\in H(Q_{2})$:

$\left\\{\begin{split}&\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{1}^{in})-p_{0}^{in}\mathbf{I}\right)=-\mathbf{e}\quad\text{in
}Q_{2}\\\
&\mathbf{n}\times\mathbf{n}\times\left[\left(\left(2\mu_{1}e(\mathbf{u}_{0}^{out})-p_{0}^{out}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}_{1}^{in})-p_{0}^{in}\mathbf{I}\right)\right)\cdot\mathbf{n}\right]=\mathbf{0}\quad\text{on
}\Gamma\end{split}\right.$

By applying Lemma 2.2 and (50) with $\mathbf{f}^{out}=-\mathbf{e}$ and
$\mathbf{f}^{in}=-\mathbf{e}$, we obtain

$\left\lVert\mathbf{u}_{1}^{in}\right\rVert_{Q_{2}}\leq{C_{1}E_{1}},\quad{C_{1}}:=\frac{\sqrt{|Q_{2}|}}{B_{1}(Q_{{2}})}+2\frac{\sqrt{|Q_{1}|}}{B_{1}(Q_{{1}})}$
(55)

  4. 4.

Induction step, $k\geq 1$: Given $\mathbf{u}_{k}^{in}\in H(Q_{2})$ and
$\mathbf{u}_{k-1}^{out}\in H(Q_{1})$, find
$\mathbf{u}_{k+1}^{in}(\mathbf{y};\mathbf{0})\in H(Q_{2})$ and
$\mathbf{u}_{k}^{out}(\mathbf{y};\mathbf{0})\in H(Q_{1})$.

    1. (a)

Applying Lemma 2.1 with $\mathbf{f}=\mathbf{0}$, we conclude that for a given
$\mathbf{u}_{k}^{in}\in H(Q_{2})$, $k\geq 1$, there exists a unique
$\mathbf{u}_{k}^{out}\in H(Q_{1})$ that solves the system of (30) and (36) and
assumes the estimate

$\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)&=\mathbf{0}\quad\text{in
}Q_{1}\\\ \mathbf{u}_{k}^{out}&=\mathbf{u}_{k}^{in}\quad\text{on
}\Gamma,\end{split}\right.$ (56)
$\left\lVert\mathbf{u}_{k}^{out}\right\rVert_{Q_{1}}\leq
2E_{1}\left\lVert\mathbf{u}_{k}^{in}\right\rVert_{Q_{2}}$ (57)

    2. (b)

By applying Lemma 2.2 with $\mathbf{f}^{in}=\mathbf{0}=\mathbf{f}^{out}$, we
see that for any given $\mathbf{u}_{k}^{out}\in H(Q_{1})$, $k\geq 1$ that
satisfies (56), there exists a unique solution $\mathbf{u}_{k+1}^{in}\in
H(Q_{2})$ to the system of equations (33) and (36)

$\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{k+1}^{in})-p_{k}^{in}\mathbf{I}\right)&=\mathbf{0}\quad\text{in
}Q_{2}\\\
\mathbf{n}\times\mathbf{n}\times\left[\left(\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}_{k+1}^{in})-p_{k}^{in}\mathbf{I}\right)\right)\mathbf{n}\right]&=\mathbf{0}\quad\text{on
}\Gamma\end{split}\right.$ (58)

Moreover,

$\left\lVert\mathbf{u}_{k+1}^{in}\right\rVert_{Q_{2}}\leq
E_{1}\left\lVert\mathbf{u}_{k}^{out}\right\rVert_{Q_{1}}$ (59)

Now we have found the coefficients
$\mathbf{u}^{in}_{n}(\mathbf{y};\mathbf{e})$ and
$\mathbf{u}^{out}_{n}(\mathbf{y};\mathbf{e})$ in (28) iteratively. We prove
the convergence of the series in the following theorem by taking into account
the fact that $\mathbf{u}_{0}^{in}=\mathbf{0}$.

######  Theorem 2.1.

Define the partial sums

$\mathbf{S}_{q}^{in}(\mathbf{y};\mathbf{e},w):=\sum_{k=0}^{q}\mathbf{u}^{in}_{k+1}(\mathbf{y},\mathbf{e})w^{k+1},\,\mathbf{S}_{q}^{out}(\mathbf{y};\mathbf{e},w):=\sum_{k=0}^{q}\mathbf{u}^{out}_{k}(\mathbf{y};\mathbf{e})w^{k}.$

Let $R\in(0,1)$, in the disk $|w|\leq\frac{R}{2E_{1}^{2}}$, the series
$\mathbf{S}_{q}^{in}(\mathbf{y};\mathbf{e},w)$ and
$\mathbf{S}_{q}^{out}(\mathbf{y};\mathbf{e},w)$ converge uniformly to
$\mathbf{u}^{in}_{\infty}(\mathbf{y};\mathbf{e},w)\in H(Q_{2})$ and
$\mathbf{u}^{out}_{\infty}(\mathbf{y};\mathbf{e},w)\in H(Q_{1})$,
respectively. Therefore,
$\mathbf{u}_{\infty}(\mathbf{y};\mathbf{e},w)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\mathbf{u}^{in}_{\infty}(\mathbf{y};\mathbf{e},w)\chi_{2}+\mathbf{u}^{out}_{\infty}(\mathbf{y};\mathbf{e},w)\chi_{1}\in
H(Q)$ solves the cell problem (13) and is analytic for
$|w|<\frac{1}{2E_{1}^{2}}$.

###### Proof.

For each $q\in\mathbb{N}$, $\mathbf{S}_{q}^{in}(\mathbf{y};\mathbf{e},w)$ is a
polynomial function of $w$ and maps from $\mathbb{C}$ to the Hilbert space
$H(Q_{2})$. Similarly, $\mathbf{S}_{q}^{out}(\mathbf{y};\mathbf{e},w)$ maps
from $\mathbb{C}$ to $H(Q_{1})$. To show uniform convergence, we note that
(57) and (59) imply there exists a positive constant $E_{1}$ that depends only
on $Q_{1}$ and $Q_{2}$ such that
$\left\lVert\mathbf{u}^{in}_{k+1}\right\rVert_{Q_{2}}\leq
E_{1}\left\lVert\mathbf{u}^{out}_{k}\right\rVert_{Q_{1}}\leq
2E_{1}^{2}\left\lVert\mathbf{u}^{in}_{k}\right\rVert_{Q_{2}}$. Therefore,

$\left\lVert\mathbf{u}_{k}^{in}\right\rVert_{Q_{2}}\leq\left(2E_{1}^{2}\right)^{k-1}\left\lVert\mathbf{u}_{1}^{in}\right\rVert_{Q_{2}},k\geq
1$ (60)

Let $m>q>N$, and define $r:=2E_{1}^{2}\left|w\right|$. Then by (55) implies

$\displaystyle\left\lVert\mathbf{S}_{m}^{in}(w)-\mathbf{S}_{q}^{in}(w)\right\rVert_{Q_{2}}\leq\left\lVert\mathbf{u}_{1}^{in}\right\rVert_{Q_{2}}\left((2E_{1}^{2})^{q}|w|^{q+1}+\cdots+(2E_{1}^{2})^{m-1}\left|w\right|^{m}\right)$
$\displaystyle\leq\frac{r^{q+1}-r^{m+1}}{1-r}\frac{\left\lVert\mathbf{u}_{1}^{in}\right\rVert_{Q_{2}}}{2E_{1}^{2}}\leq\frac{r^{q+1}-r^{m+1}}{1-r}\left({\frac{C_{1}}{2E_{1}^{2}}}\right)\mbox{
{,\, $C_{1}$ is defined in \eqref{bound_on_u1_in}}}.$

Therefore, for $r\leq R<1$, i.e. $|w|\leq\frac{R}{2E_{1}^{2}}$, where $R$ is
any fixed number in $(0,1)$,

$\left\lVert\mathbf{S}_{m}^{in}(w)-\mathbf{S}_{q}^{in}(w)\right\rVert_{Q_{2}}\leq\left({\frac{C_{1}}{2E_{1}^{2}}}\right)\left(\frac{R^{N+1}}{1-R}\right),\,\,\forall
m>q>N.$ (61)

For $\mathbf{S}_{q}^{out}(\mathbf{y};\mathbf{e},w)$ we have
$\left\lVert\mathbf{u}_{k}^{out}\right\rVert_{Q_{1}}\leq\left(2E_{1}^{2}\right)^{k-1}\left\lVert\mathbf{u}_{1}^{out}\right\rVert_{Q_{1}}$.
By a similar procedure, for $m>q>N$ and $|w|\leq\frac{R}{2E_{1}^{2}}$ the
following estimate is valid

$\left\lVert\mathbf{S}_{m}^{out}(w)-\mathbf{S}_{q}^{out}(w)\right\rVert_{Q_{1}}\leq{C_{1}}\left(\frac{R^{N+1}}{1-R}\right)$

Therefore, for every fixed $w$ satisfying $|w|\leq\frac{R}{2E_{1}^{2}}$ for
any $0<R<1$, $\mathbf{S}_{q}^{in}(\mathbf{y};w)$ and
$\mathbf{S}_{q}^{out}(\mathbf{y};w)$ converge uniformly to
$\mathbf{u}^{in}_{\infty}(\mathbf{y};w)\in H(Q_{2})$ and
$\mathbf{u}^{out}_{\infty}(\mathbf{y};w)\in H(Q_{1})$, respectively. Since for
each $q$, $\mathbf{S}_{q}^{in}(\mathbf{y};w)$ and
$\mathbf{S}_{q}^{out}(\mathbf{y};w)$ are polynomials of $w$, hence analytic,
the uniform convergence implies that the limit functions
$\mathbf{u}^{in}_{\infty}(\mathbf{y};w)$ and
$\mathbf{u}^{out}_{\infty}(\mathbf{y};w)$ are also analytic in
$|w|<\frac{1}{2E_{1}^{2}}$ with values in $H(Q_{1})$ and $H(Q_{2})$,
respectively, by applying Morera’s theorem for Banach space valued analytic
functions [27] to the uniformly converging sequences. By construction, the
function $\mathbf{u}_{\infty}(\mathbf{y};\mathbf{e},w)$ defined in (27) solves
the cell problem (13) for all $w$ in the disk
$\\{w:|w|<\frac{1}{2E_{1}^{2}}\\}:=B_{0}(\frac{1}{2E_{1}^{2}})$. Moreover, the
uniqueness of the solution implies that
$\mathbf{u}_{\infty}(\mathbf{y};\mathbf{e}_{k},w)=\mathbf{u}^{k}(\mathbf{y};\frac{1}{w})\text{
in }H(Q)$ for $w\in
B_{0}(\frac{1}{2E_{1}^{2}})\cap\\{w\in\mathbb{C}\setminus(-\infty,0]\\}$. ∎

The following theorem shows the relation between the two-fluid self-
permeability $K$ in (12) and the Darcy permeability
$\mbox{\boldmath$K$}^{(D)}$ in (5)

######  Theorem 2.2.

In the case of large viscosity $|z|>2E_{1}^{2}$ (or
$|w|<\frac{1}{2E_{1}^{2}}$), we have

  1. 1.

$\mathbf{u}_{\infty}^{in}(\mathbf{y};\mathbf{e}_{i},0)=\mathbf{0}$ in $Q_{2}$

  2. 2.

As $w\rightarrow 0$, the solution
$\mathbf{u}_{\infty}^{out}(\mathbf{y};\mathbf{e}_{i},w)$ converges uniformly
in $\mathring{H}(Q_{1})$ to the solution $\mathbf{u}_{D}^{i}(\mathbf{y})$ of
the classical cell problem (4).

  3. 3.

For $w\in B_{0}(\frac{1}{2E_{1}^{2}})$, the difference between the self-
permeability $\mathbf{K}(\mathbf{y};\mathbf{e}_{i},w)$ and the classical
permeability tensor $\mathbf{K}^{(D)}(\mathbf{y};\mathbf{e}_{i})$ satisfies
$\lvert K_{ij}-(K^{(D)})_{ij}\rvert=O(|w|)$, hence
$\mathbf{K}\to\mathbf{K}^{(D)}$ uniformly as $|w|\to 0$.

###### Proof.

The uniform convergence allows passing the limit $w\rightarrow 0$ inside the
summation of (28) to obtain
$\mathbf{u}_{\infty}^{in}(\mathbf{y};\mathbf{e}_{i},0)=\mathbf{0}$. Similarly,
the uniform convergence allows passing the limit $w\rightarrow 0$ inside the
summation of (28) to obtain

$\mathbf{u}^{out}_{\infty}(\mathbf{y};\mathbf{e}_{{i}},0)=\mathbf{u}_{0}^{out}(\mathbf{y};\mathbf{e}_{i})$

Furthermore, $\mathbf{u}_{0}^{out}(\mathbf{y};\mathbf{e}_{i})\in H(Q_{1})$
satisfies (53) and in fact
$\mathbf{u}_{0}^{out}(\mathbf{y};\mathbf{e}_{i})\in\mathring{H}(Q_{1})$ since
$\mathbf{u}_{0}^{out}(\mathbf{y};\mathbf{e}_{i})|_{\Gamma}=\mathbf{0}$, which
is identical to the equation for $\mathbf{u}_{D}$ (4). The uniqueness of the
solution then ensures that
$\mathbf{u}_{0}^{out}(\mathbf{y};\mathbf{e}_{i})=\mathbf{u}_{D}^{i}(\mathbf{y})$.
Therefore the series
$\mathbf{u}_{\infty}^{out}(\mathbf{y};\mathbf{e}_{i},w)\to\mathbf{u}_{D}^{i}(\mathbf{y})$
uniformly as $|w|\to 0$ in $\mathring{H}(Q_{1})$. For (iii), we note that

$\begin{split}\left|K_{ij}(w)-K^{(D)}_{ij}\right|&=\left|\int_{Q}\left(\mathbf{u}^{i}-\chi_{1}\mathbf{u}_{D}^{i}\right)\cdot\mathbf{e}_{j}d\mathbf{y}\right|\leq\left\lVert\mathbf{u}^{i}-\chi_{1}\mathbf{u}^{i}_{D}\right\rVert_{L^{2}(Q)}\\\
&\leq\frac{1}{B_{1}(Q)}\left\lVert\sum_{k=1}^{\infty}\left(\mathbf{u}_{k}^{in}(\mathbf{y};\mathbf{e}_{i})\chi_{2}+\mathbf{u}_{k}^{out}(\mathbf{y};\mathbf{e}_{i})\chi_{1}\right)w^{k}\right\rVert_{Q}\end{split}$

From (57), (60) and (55), we have for $|w|<\frac{1}{2E_{1}^{2}}$, or
equivalently $|z|>2E_{1}^{2}$,

$\left|K_{ij}(w)-K^{(D)}_{ij}\right|\leq{C_{1}\left(\frac{E_{1}+1}{2E_{1}B_{1}(Q)}\right)}\frac{2E_{1}^{2}|w|}{1-2E_{1}^{2}|w|}$
(62)

∎

In the following section, we study the behavior of $\mbox{\boldmath$K$}(z)$
near $z=0$, i.e. the inclusion is an air bubble.

###  2.4 Analyticity of the solution for small $|z|$

Let $\mathbf{e}$ be a constant unit vector in $\mathbb{R}^{n}$. We seek
solutions of the following form

$\displaystyle\mathbf{u}_{null}^{in}(\mathbf{y};\mathbf{e},z)$
$\displaystyle=\sum_{k=0}^{\infty}\mathbf{u}^{in}_{k}(\mathbf{y};\mathbf{e})z^{k},\quad
p^{in}(\mathbf{y};\mathbf{e},z)=\sum_{k=0}^{\infty}p^{in}_{k}(\mathbf{y};\mathbf{e})z^{k}\text{
in }Q_{2},$ (63)
$\displaystyle\mathbf{u}_{null}^{out}(\mathbf{y};\mathbf{e},z)$
$\displaystyle=\sum_{k=0}^{\infty}\mathbf{u}^{out}_{k}(\mathbf{y};\mathbf{e})z^{k},\quad
p^{out}(\mathbf{y};\mathbf{e},z)=\sum_{k=0}^{\infty}p^{out}_{k}(\mathbf{y};\mathbf{e})z^{k}\text{
in }Q_{1}$ (64)

By a procedure similar to that in Section 2.3, the following equations are
obtained via collecting terms with respect to the order of $z$. The PDEs for
$Q_{1}$ are as follows.

$\displaystyle O(1)$ $\displaystyle:$
$\displaystyle\qquad\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{0}^{out})-p_{0}^{out}\mathbf{I}\right)=-\mathbf{e}$
(65) $\displaystyle O(z^{k}),\,k\geq 1$ $\displaystyle:$
$\displaystyle\qquad\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)=\mathbf{0}$
(66)

Similarly, the PDEs for $Q_{2}$ are

$\displaystyle O(1)$ $\displaystyle:$ $\displaystyle\qquad-\nabla
p_{0}^{in}=-\mathbf{e}$ (67) $\displaystyle O(z^{k}),\,k\geq 1$
$\displaystyle:$
$\displaystyle\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}_{k-1}^{in})-p_{k}^{in}\mathbf{I}\right)=\mathbf{0}$
(68)

The interface condition (36) remains the same for the small $|z|$ case while
(34) and (36) now read

$\displaystyle\mathbf{n}\times\mathbf{n}\times\left[\left(\left(2\mu_{1}e(\mathbf{u}_{0}^{out})-p_{0}^{out}\mathbf{I}\right)\cdot\mathbf{n}-(-p_{0}^{in}\mathbf{I}\right)\cdot\mathbf{n}\right]=\mathbf{0}$
(69)
$\displaystyle\mathbf{n}\times\mathbf{n}\times\left[\left(\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}_{k-1}^{in})-p_{k}^{in}\mathbf{I}\right)\right)\cdot\mathbf{n}\right]=\mathbf{0},\,k\geq
1$ (70)

The first equation to be solved is (67), whose solution is simply

$p_{0}^{in}(\mathbf{y})=\mathbf{e}\cdot\mathbf{y}+{c}-\int_{Q_{2}}(\mathbf{e}\cdot\mathbf{y}+{c})d\mathbf{y}\mbox{
in }Q_{2}$ (71)

where ${c}$ is a constant. The next problem is the system of (65) and (69).
Similar to the calculation in Lemma 2.2, the weak formulation of this system
is: Find $u_{0}^{out}\in H(Q_{1})$ such that for all
$\mbox{\boldmath$\Phi$}\in H(Q_{1})$ and
$\pi_{0}^{out}:=2\mu_{1}e(\mathbf{u}_{0}^{out})-p_{0}^{out}\mathbf{I}$

$-\int_{\Gamma}(\bar{\Phi}\cdot(\pi_{0}^{out}+p_{0}^{in}\mathbf{I})-p_{0}^{in}\mathbf{I})\cdot\mathbf{n}dS-\int_{Q_{1}}2\mu_{1}e(\mathbf{u}_{0}^{out}):\overline{e(\Phi)}d\mathbf{y}=\int_{Q_{1}}-\mathbf{e}\cdot\bar{\Phi}d\mathbf{y}$

Since $\mbox{\boldmath$\Phi$}\cdot\mathbf{n}=0$ and
$p_{0}^{in}\mathbf{I}\cdot\mathbf{n}$ is parallel to $\mathbf{n}$, (69)
implies the integral on $\Gamma$ vanishes. Hence by the Lax-Milgram lemma, we
have

$\|\mathbf{u}_{0}^{out}\|_{Q_{1}}\leq\frac{\sqrt{|Q_{1}|}}{B_{1}(Q_{{1}})}$
(72)

The system for $\mathbf{u}_{k-1}^{in}$, $k\geq 1$ (inner problem) is to find
$\mathbf{u}_{k-1}^{in}\in H(Q_{2})$ with given $\mathbf{u}_{0}^{out}\in
H(Q_{1})$ such that

$\left\\{\begin{split}\text{div}\left(2\mu_{1}e(\mathbf{u}_{k-1}^{in})-p_{k}^{in}\mathbf{I}\right)=\mathbf{0}\mbox{
in }Q_{2}\\\
\mathbf{u}_{k-1}^{in}|_{\Gamma}=\mathbf{u}_{k-1}^{out}|_{\Gamma}\end{split}\right.$
(73)

With an argument similar to the derivation of Lemma 2.1, the following
estimate can be derived for system (73)

######  Lemma 2.3.

Let $Q_{2}$ satisfy the same assumption in Lemma 2.1. For any given vector
field $\mathbf{u}^{out}\in H(Q_{1})$, there exists a unique weak solution
$\mathbf{u}^{in}(\mathbf{y})\in H(Q_{2})$ s.t.

$\displaystyle\left\\{\begin{split}\text{div}_{\mathbf{y}}\left(2\mu_{1}e(\mathbf{u}^{in})-p^{in}\mathbf{I}\right)&=\mathbf{f}^{in}\quad\text{
in }Q_{2}\\\ \mathbf{u}^{in}&=\mathbf{u}^{out}\quad\text{ on
}\Gamma\end{split}\right.$ (74)
$\displaystyle\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}\leq\frac{1}{B_{1}(Q_{{2}})}\left\lVert\mathbf{f}^{in}\right\rVert_{L^{2}(Q_{2})}+2E_{2}\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}.$
(75)

where $E_{2}>1$ is the constant associated with the extension operator $T$,
$\|T(\mbox{\boldmath$\Phi$})\|_{Q}\leq
E_{2}\|\mbox{\boldmath$\Phi$}\|_{Q_{{1}}}$ for all $\mbox{\boldmath$\Phi$}\in
H(Q_{{1}})$ and $T(\mbox{\boldmath$\Phi$})$ decays rapidly to 0 inside
$Q_{2}$. Note that the periodic condition of space $H(Q_{1})$ implies
$\int_{\Gamma}\mathbf{u}^{out}\cdot\mathbf{n}\,dS=0$.

The system for $\mathbf{u}_{k}^{out}$ and $p_{k}^{out}$ with given
$\mathbf{u}_{k-1}^{in}\in H(Q_{2})$ and $p_{k}^{in}$, $k\geq 1$ is

$\left\\{\begin{split}\text{div}\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)=\mathbf{0}\\\
\mathbf{n}\times\mathbf{n}\times\left[\left(\left(2\mu_{1}e(\mathbf{u}_{k}^{out})-p_{k}^{out}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}_{k-1}^{in})-p_{k}^{in}\mathbf{I}\right)\right)\cdot\mathbf{n}\right]=\mathbf{0}\end{split}\right.$
(76)

By an argument similar to the one for Lemma 2.2, the system above can be shown
to satisfy the following estimate.

######  Lemma 2.4.

Let $Q_{2}$ satisfy the same assumption in Lemma 2.1. For any given pair of
$\left(\mathbf{u}^{in},p^{in}\right)\in H(Q_{2})\times
L^{2}(Q_{1})/\mathbb{C}$ that satisfies (74), there exists a unique vector
$\mathbf{u}^{out}(\mathbf{y};\mathbf{f}^{out})\in H(Q_{1})$ solving the
following system

$\displaystyle\left\\{\begin{split}&\text{div}\left(2\mu_{1}e(\mathbf{u}^{out})-p^{out}\mathbf{I}\right)=\mathbf{f}^{out}\text{
in }Q_{1}\\\
&\mathbf{n}\times\mathbf{n}\times\left[\left(\left(2\mu_{1}e(\mathbf{u}^{in})-p^{in}\mathbf{I}\right)-\left(2\mu_{1}e(\mathbf{u}^{out})-p^{out}\mathbf{I}\right)\right)\cdot\mathbf{n}\right]=\mathbf{0}\text{
on }\Gamma\end{split}\right.$ (77)
$\displaystyle\left\lVert\mathbf{u}^{out}\right\rVert_{Q_{1}}\leq\frac{E_{2}}{B_{1}(Q_{{1}})}\left\lVert\mathbf{f}^{out}\right\rVert_{L^{2}(Q_{1})}+\frac{E_{2}}{B_{1}(Q_{{2}})}\left\lVert\mathbf{f}^{in}\right\rVert_{L^{2}(Q_{2})}+E_{2}\left\lVert\mathbf{u}^{in}\right\rVert_{Q_{2}}$
(78)

where $E_{2}$, $B_{1}$ depend only on $Q$ and $\Gamma$.

Equation (76), Lemma 2.3, Equation (73) and Lemma 2.4 imply that for all
$k\geq 0$, we have $\|\mathbf{u}_{k}^{in}\|_{Q_{2}}\leq
2E_{2}\|\mathbf{u}_{k}^{out}\|_{Q_{1}}$ and
$\|\mathbf{u}_{k+1}^{out}\|_{Q_{1}}\leq E_{2}\|\mathbf{u}_{k}^{in}\|_{Q_{2}}$.
Therefore,

$\displaystyle\|\mathbf{u}^{in}_{k}\|_{Q_{2}}\leq\frac{(2E_{2}^{2})^{k+1}}{E_{2}}\|\mathbf{u}_{0}^{out}\|_{Q_{1}}\leq{(2E_{2}^{2})^{k+1}}\left(\frac{|Q_{1}|}{E_{2}B_{1}(Q_{{1}})}\right)$
(79)
$\displaystyle\|\mathbf{u}^{out}_{k}\|_{Q_{1}}\leq(2E_{2}^{2})^{k}\|\mathbf{u}_{0}^{out}\|_{Q_{1}}\leq(2E_{2}^{2})^{k}\left(\frac{|Q_{1}|}{B_{1}(Q_{1})}\right)$
(80)

Therefore, the series in (63) and (64) converge uniformly in the disk
$|z|<\frac{1}{2(E_{2})^{2}}$ to an analytic function in $Q_{2}$ and $Q_{1}$,
respectively. The limit functions
$\mathbf{u}_{null}^{in}(\mathbf{y},\mathbf{e},z)$,
$\mathbf{u}_{null}^{out}(\mathbf{y},\mathbf{e},z)$ and the corresponding
permeability $K_{ij}(z)$ in (19) are analytic at $z=0$. Define the
permeability (’B’ for ’bubbles)

$K_{ij}^{(B)}:=\int_{Q}[\chi_{1}\mathbf{u}_{0}^{out}(\mathbf{y};\mathbf{e}_{i})+\chi_{2}\mathbf{u}_{0}^{in}(\mathbf{y};\mathbf{e}_{i})]\cdot\mathbf{e}_{j}\,d\mathbf{y}$
(81)

then the following estimate, valid for $|z|<\frac{1}{2E_{2}^{2}}$, holds

$\displaystyle|K_{ij}(z)-K^{B}_{ij}|\leq\frac{|Q_{1}|(1+2E_{2})}{{B_{1}(Q)B_{1}(Q_{1})}}\left(\frac{2E_{2}^{2}|z|}{1-2E_{2}^{2}|z|}\right)=O(|z|).$
(82)

In conclusion, $\mbox{\boldmath$K$}(z)$ in (12) is analytic for
$z\in\mathbb{C}\setminus[-2E_{1}^{2},-\frac{1}{2E_{2}^{2}}]$, $E_{1},E_{2}\geq
1$. In the next section, and integral representation formula (IRF) for
$\mbox{\boldmath$K$}(z)$ will be derived in two different ways.

##  3 Integral representation of permeability $\mbox{\boldmath$K$}(z)$

We first observe two properties of $K$ implied by (19).

######  Proposition 3.1.

$\displaystyle\frac{\mbox{\boldmath$K$}(z)-\mbox{\boldmath$K$}^{*}(z)}{z-\bar{z}}\leq
0\mbox{ if }Im(z)\neq 0$ (83) $\displaystyle\mbox{\boldmath$K$}(x)\geq 0\mbox{
for }x>0$ (84)

###### Proof.

Note that
$K_{ij}(z)-(K^{*})_{ij}(z)={2\mu_{1}}(\overline{z}-z)\int_{Q_{2}}{e(\mathbf{u}^{j}(z))}:\overline{e(\mathbf{u}^{i}(z))}\,d\mathbf{y}$.
Hence

$\frac{K_{ij}(z)-K^{*}_{ij}(z)}{z-\overline{z}}={-}{2\mu_{1}}\int_{Q_{2}}{e(\mathbf{u}^{j}(z))}:\overline{e(\mathbf{u}^{i}(z))}\,d\mathbf{y}={-}(\mathbf{u}^{j},\mathbf{u}^{i})_{Q_{2}}=:-A_{ij}$

The matrix $\boldsymbol{A}$ is obviously Hermitian and for any
$\boldsymbol{\xi}\in\mathbb{C}^{n}$, we have
$\overline{\xi_{i}}A_{ij}{\xi_{j}}=(\xi_{j}\mathbf{u}^{j},\xi_{i}\mathbf{u}^{i})_{Q_{2}}\geq
0$. This proves (83). Recall that
$K_{ij}(x)=\left((\mathbf{u}^{j},\mathbf{u}^{i})_{Q_{1}}+x(\mathbf{u}^{j},\mathbf{u}^{i})_{Q_{2}}\right)$.
With a similar argument, (84) follows. ∎

With these two properties and the fact that $K$ is holomorphic in
$\mathbb{C}\setminus(-\infty,0]$, the characterization theorem for matrix-
valued functions belonging to the Stieltjes class [20] implies that there
exists a monotonically increasing matrix-valued function
$\boldsymbol{\sigma}(t)$ such that the following integral representation
formula (IRF) holds for $z\in\mathbb{C}\setminus(-\infty,0]$

$\mbox{\boldmath$K$}(z)=\boldsymbol{A}+\frac{\boldsymbol{C}}{z}+\int_{+0}^{\infty}\frac{1}{z+t}d\boldsymbol{\sigma}(t)$

where $\boldsymbol{A}\geq 0$, $\boldsymbol{C}\geq 0$,
$\int_{+0}^{\infty}\frac{1}{1+t}d\boldsymbol{\sigma}(t){:=\displaystyle{\lim_{\epsilon\downarrow
0}\int_{\epsilon}^{\infty}}\frac{1}{1+t}d\boldsymbol{\sigma}(t)}<\infty$ and
$\boldsymbol{A}+\boldsymbol{C}+\int_{+0}^{\infty}\frac{1}{1+t}d\boldsymbol{\sigma}(t)>0$.
Since $\mbox{\boldmath$K$}(0)=\mbox{\boldmath$K$}^{(B)}$, we must have
$\boldsymbol{C}=\boldsymbol{0}$. Also,
$\boldsymbol{K}(\infty)=\boldsymbol{K}^{(D)}$ implies
$\boldsymbol{A}=\boldsymbol{K}^{(D)}$

$\mbox{\boldmath$K$}(z)=\boldsymbol{K}^{(D)}+\int_{\frac{1}{2E_{2}^{2}}}^{2E_{1}^{2}}\frac{1}{z+t}d\boldsymbol{\sigma}(t)$

Therefore, for real valued $z$, $\mbox{\boldmath$K$}(z)$ is decreasing as $z$
increases, i.e., $\mbox{\boldmath$K$}(x_{1})-\mbox{\boldmath$K$}(x_{2})$ is
negative semidefinite if $x_{1}>x_{2}$. To study how the measure
$d\boldsymbol{\sigma}$ is related to the microstructure, we derive the
spectral representation of $\mbox{\boldmath$K$}(z)$ by using the underlying
system (13).

###  3.1 Spectral representation of $\mbox{\boldmath$K$}(z)$

Adding
$\int_{Q_{2}}2\mu_{1}e(\mathbf{u}^{k}):\overline{e(\mathbf{v})}d\mathbf{y}$ to
both sides of (21), we have

$\int_{Q}2\mu_{1}e(\mathbf{u}^{k}):\overline{e(\mathbf{v})}d\mathbf{y}=-\frac{1}{s}\int_{Q}2\mu_{1}\chi_{2}e(\mathbf{u}^{k}):\overline{e(\mathbf{v})}d\mathbf{y}+\int_{Q}\mathbf{e}_{k}\cdot\bar{\mathbf{v}}d\mathbf{y}$
(85)

where the new variable $s$ is defined as

$s:=\frac{1}{z-1}$

Let $\Delta_{\\#}^{-1}$ be the operator that solves for
$\mathbf{w}(\mathbf{y};\mathbf{f})\in H(Q)$ in the following variational
formulation

$\int_{Q}2\mu_{1}e({\mathbf{w}}):\overline{e(\mathbf{v})}d\mathbf{y}=\int_{Q}\mathbf{f}\cdot\bar{\mathbf{v}}d\mathbf{y}$
(86)

where $\mathbf{f}\in L^{2}(Q)$ and $Q$-periodic. In other words, solution
${\mathbf{w}}(\mathbf{y})=\Delta_{\\#}^{-1}\mathbf{f}\in H(Q)$ is a weak
solution to the cell problem

$\left\\{\begin{split}-\mu_{1}\Delta{\mathbf{w}}&=\mathbf{f}\quad\text{in
}Q_{1}\cup Q_{2}\\\
\llbracket\mbox{\boldmath$\pi$}\rrbracket\mathbf{n}&=\left(\llbracket\mbox{\boldmath$\pi$}\mathbf{n}\rrbracket\cdot\mathbf{n}\right)\mathbf{n}\text{
on }\Gamma\end{split}\right.$ (87)

In order to get the spectral representation, we apply $\Delta_{\\#}^{-1}$ on
both sides of (85) and symbolically represent the resulted equations as

${\mathbf{w}}_{1}={-\frac{1}{s}}{\mathbf{w}}_{2}+{\mathbf{w}}_{3}$

Then clearly, we have ${\mathbf{w}}_{1}=\mathbf{u}^{k}$ and
${\mathbf{w}}_{3}=\Delta_{\\#}^{-1}\mathbf{e}_{k}$. Observe that
${\mathbf{w}}_{2}$ solves

$\int_{Q}2\mu_{1}e({\mathbf{w}}_{2}):\overline{e(\mathbf{v})}d\mathbf{y}=\int_{Q}2\mu_{1}\chi_{2}e(\mathbf{u}^{k}):\overline{e(\mathbf{v})}d\mathbf{y}\mbox{
for all }\mathbf{v}\in H(Q)$ (88)

Define the operator $\Gamma_{\chi}$ such that
${\mathbf{w}}_{2}=\Gamma_{\chi}\mathbf{u}^{k}$ and (88) can be expressed as

$(\Gamma_{\chi}\mathbf{u}^{{k}},\mathbf{v})_{Q}=\int_{Q}2\mu_{1}\chi_{2}e(\mathbf{u}^{{k}}):\overline{e(\mathbf{v})}d\mathbf{y}\mbox{
for all }\mathbf{v}\in H(Q).$ (89)

The subscript $\chi$ is used to signify the dependence of $\Gamma_{\chi}$ on
$\chi_{2}$, the characteristic function of $Q_{2}$. Clearly, $\Gamma_{\chi}$
is self-adjoint with respect to the inner product $(\cdot,\cdot)_{Q}$ because

$(\Gamma_{\chi}\mathbf{u},\mathbf{v})_{Q}=\overline{\int_{Q}2\mu_{1}\chi_{2}e(\mathbf{v}):\overline{e(\mathbf{u})}d\mathbf{y}}=\overline{(\Gamma_{\chi}\mathbf{v},\mathbf{u})_{Q}}={(\mathbf{u},\Gamma_{\chi}\mathbf{v})_{Q}.}$

Formally, we have
$\Gamma_{\chi}\mathbf{u}=\triangle_{\\#}^{-1}(\nabla\cdot\chi_{2}e(\mathbf{u}))$.
Now (85) becomes

$\mathbf{u}^{k}=-\frac{1}{s}\Gamma_{\chi}\mathbf{u}^{k}+\Delta_{\\#}^{-1}\mathbf{e}_{k}\Leftrightarrow\left(I+\frac{\Gamma_{\chi}}{s}\right)\mathbf{u}^{k}=\Delta_{\\#}^{-1}\mathbf{e}_{k}$
(90)

######  Proposition 3.2.

The self-adjoint operator $\Gamma_{\chi}$ defined in (89) is positive and
bounded with $\left\lVert\Gamma_{\chi}\right\rVert\leq 1$.

###### Proof.

It can be proved by choosing $\mathbf{v}=\mathbf{u}$ in (89) and observe that
$0\leq\int_{Q_{2}}2\mu_{1}e(\mathbf{u}):\overline{e(\mathbf{u})}d\mathbf{y}\leq\int_{Q}2\mu_{1}e(\mathbf{u}):\overline{e(\mathbf{u})}d\mathbf{y}=(\mathbf{u},\mathbf{u})_{Q}$.
∎

######  Theorem 3.1.

For $|s|>1$, the solution $\mathbf{u}^{k}\in H(Q)$ admits a series
representation

$\mathbf{u}^{k}(\mathbf{y};s)=\sum_{m=0}^{\infty}\left(-\frac{1}{s}\right)^{m}\left(\Gamma_{\chi}\right)^{m}\Delta_{\\#}^{-1}\mathbf{e}_{k}$
(91)

and the components of $K$ can be represented by the following IRF

$K_{kl}(s)=s\int_{0}^{1}\int_{Q}\frac{\left(\tilde{M}(d\lambda)\Delta_{\\#}^{-1}\mathbf{e}_{k}\right)_{l}}{s+\lambda}\,d\mathbf{y},\quad
k,l=1,\dots,n,$ (92)

for some projection-valued measures $\tilde{M}(d\lambda)$ and a series
representation

$K_{kl}(s)=\int_{Q}\left(\Delta_{\\#}^{-1}\mathbf{e}_{k}\right)_{l}d\mathbf{y}+\sum_{m=1}^{\infty}\frac{\tilde{\lambda}_{kl}^{m}}{(-s)^{m}}\quad\mbox{with
}\tilde{\lambda}_{kl}^{m}:=\int_{Q}\left(\left(\Gamma_{\chi}\right)^{m}\Delta_{\\#}^{-1}\mathbf{e}_{k}\right)_{l}\,d\mathbf{y}.$

###### Proof.

From (90), since $\Gamma_{\chi}$ is self-adjoint with norm bounded by 1, for
$|s|>1$, the spectral theory for self-adjoint operator implies the existence
of a projection-valued measure $\tilde{M}$ such that

$\mathbf{u}^{k}(\mathbf{y};s)=\left(I+\frac{\Gamma_{\chi}}{s}\right)^{-1}\Delta_{\\#}^{-1}\mathbf{e}_{k}={s}\int_{0}^{1}\frac{\tilde{M}(d\lambda)\left(\Delta_{\\#}^{-1}\mathbf{e}_{k}\right)}{s+\lambda}$
(93)

Hence the $kl$-the element of permeability $K$ has the following IRF

$K_{kl}(s)=\int_{Q}(\mathbf{u}^{k})_{l}d\mathbf{y}=s\int_{0}^{1}\int_{Q}\frac{\left(\tilde{M}(d\lambda)\Delta_{\\#}^{-1}\mathbf{e}_{k}\right)_{l}}{s+\lambda}d\mathbf{y}$
(94)

On the other hand, for $|s|>1$, the geometric expansion of the middle termnear
$s=\infty$ in (93) results in the following expression

$K_{kl}(s)=\int_{Q}\left[\sum_{m=0}^{\infty}\left(-\frac{1}{s}\right)^{m}\left(\Gamma_{\chi}\right)^{m}\Delta_{\\#}^{-1}\mathbf{e}_{k}\right]\cdot\mathbf{e}_{l}d\mathbf{y}=\sum_{m=0}^{\infty}\frac{\tilde{\lambda}_{kl}^{m}}{(-s)^{m}}$
(95)

where $\tilde{\lambda}_{kl}^{m}$ is defined as
$\tilde{\lambda}_{kl}^{m}:=\int_{Q}\left(\left(\Gamma_{\chi}\right)^{m}\Delta_{\\#}^{-1}\mathbf{e}_{k}\right)_{l}d\mathbf{y}.$
∎

For the three-dimensional space $n=3$, the expansion (94) can be cast in the
matrix form

$\mbox{\boldmath$K$}(s)=\sum_{m=0}^{\infty}\frac{\tilde{\mbox{\boldmath$\Lambda$}}_{m}}{(-s)^{m}}$
(96)

with the matrix-valued moments defined as

$\tilde{\mbox{\boldmath$\Lambda$}}_{m}:=\begin{pmatrix}\int_{Q}\left(\Gamma_{\chi}\right)^{m}\Delta_{\\#}^{-1}\mathbf{e}_{1}d\mathbf{y}&\int_{Q}\left(\Gamma_{\chi}\right)^{m}\Delta_{\\#}^{-1}\mathbf{e}_{2}d\mathbf{y}&\int_{Q}\left(\Gamma_{\chi}\right)^{m}\Delta_{\\#}^{-1}\mathbf{e}_{3}d\mathbf{y}\end{pmatrix}$
(97)

###  3.2 Relationships between two representations and characterization of the
microstructral information on permeability

The calculations in the previous section reveal that the variable
$s:=\frac{1}{z-1}$ is the natural one to use. Because of this, we will
consider $K$ as a function of $s$. Note that $s$ maps $(-\infty,0]$ on the
$z$-plane to $[-1,0]$ on the $s$-plane.The following properties of
$\mbox{\boldmath$K$}(s)$ can be easily deduced from the results in Proposition
3.1 .

  1. 1.

$\mbox{\boldmath$K$}(s)$ is holomorphic in
$\mathbb{C}\setminus{[-\frac{2E_{2}^{2}}{1+2E_{2}^{2}},-\frac{1}{1+2E_{1}^{2}}]}.$

  2. 2.

$\frac{\mbox{\boldmath$K$}(s)-(\mbox{\boldmath$K$}(s))^{*}}{s-\overline{s}}\geq
0$ for all $Im{(s)}\neq 0$

  3. 3.

$\mbox{\boldmath$K$}(s)\geq 0\mbox{ for }\mathbb{R}\ni s>0$ because $s>0$ iff
$\mathbb{R}\ni z>1$.

Then by the representation theorem in [20, Theorem 3.1], there exists a
monotonically increasing matrix-valued function $\boldsymbol{\sigma}(t)$,
matrices $\boldsymbol{A}\geq 0$ and $\boldsymbol{C}\geq 0$ such that
$\int_{+0}^{\infty}\frac{d\boldsymbol{\sigma}}{1+t}<\infty$,
$\boldsymbol{A}+\boldsymbol{C}+\int_{+0}^{\infty}\frac{d\boldsymbol{\sigma}}{1+t}>0$
and

$\mbox{\boldmath$K$}(s)=\boldsymbol{A}+{\boldsymbol{C}}{s}+\int_{+0}^{\infty}\frac{s}{s+t}d\boldsymbol{\sigma}(t),$
(98)

As $s\rightarrow\infty$, $z\rightarrow 1$ and hence
$\mbox{\boldmath$K$}\rightarrow\mbox{\boldmath$K$}(z=1)$. Therefore, we must
have $\boldsymbol{C}=\boldsymbol{0}$. Moreover,
$\boldsymbol{A}=\mbox{\boldmath$K$}(s=0)=\mbox{\boldmath$K$}^{(D)}$. Also,
since $\mbox{\boldmath$K$}(s)$ is holomorphic in
$\mathbb{C}\setminus{[-\frac{2E_{2}^{2}}{1+2E_{2}^{2}},-\frac{1}{1+2E_{1}^{2}}]}$,
we have

$\mbox{\boldmath$K$}(s)=\boldsymbol{K}^{(D)}+\int_{\frac{1}{1+2E_{1}^{2}}}^{\frac{2E_{2}^{2}}{1+2E_{2}^{2}}}\frac{s}{s+t}d\boldsymbol{\sigma}(t),$
(99)

which is valid for all
$s\in\mathbb{C}\setminus{[-\frac{2E_{2}^{2}}{1+2E_{2}^{2}},-\frac{1}{1+2E_{1}^{2}}]}\subset(-1,0)$.
To compare with (96), which is valid only for $|s|>1$, we expand (99) near
$s=\infty$ to obtain the following series expansion

$\mbox{\boldmath$K$}(s)=\mbox{\boldmath$K$}^{(D)}+\sum_{m=0}^{\infty}(-1)^{m}\left(\frac{1}{s}\right)^{m+1}\boldsymbol{\mu}^{\sigma}_{m}$
(100)

where $\boldsymbol{\mu}^{\sigma}_{m}$ is the $m$-th moment of the measure
$d\boldsymbol{\sigma}$. Equating the coefficients term by term with (96) leads
to the following relation between $\boldsymbol{\mu}^{\sigma}_{m}$ and the
’geometrical information’ coefficients in (97)

$\displaystyle\mbox{\boldmath$K$}^{(D)}+\boldsymbol{\mu}^{\sigma}_{0}=\tilde{\mbox{\boldmath$\Lambda$}}_{0}={\mbox{\boldmath$K$}(s=\infty)},\mbox{
i.e. }$
$\displaystyle\boldsymbol{\mu}^{\sigma}_{0}=\mbox{\boldmath$K$}(z=1)-\mbox{\boldmath$K$}^{(D)}$
(102)
$\displaystyle\boldsymbol{\mu}^{\sigma}_{m}=\tilde{\mbox{\boldmath$\Lambda$}}_{m},\,m\geq
1$

Recall that $K$ can be regarded as a function of $s$ as well as a function of
$z$, $s:=\frac{1}{z-1}$. In particular, the first moment
$\boldsymbol{\mu}^{\sigma}_{1}$ can be calculated explicitly as follows

$\displaystyle\tilde{\lambda}_{kl}^{1}=(\Gamma_{\chi}\mathbf{u}^{k}(\mathbf{y};1),\mathbf{u}^{l}(\mathbf{y};1))_{Q}=2\mu_{1}\int_{Q}\chi_{2}e(\mathbf{u}^{k}(\mathbf{y};1)):\overline{e(\mathbf{u}^{l}(\mathbf{y};1))}\,d\mathbf{y}$
(103)

##  4 Numerical verification

The computational domain with $Q=(0,1)^{2}$, $Q_{2}=[1/4,3/4]^{2}$ and
$Q_{1}=Q\setminus Q_{2}$ is illustrated in Figure 2. is chosen in the first
two numerical examples, (104) and (105).

$Q_{1}$$\tilde{\Gamma}$$Q:$$Q_{2}$${\bf{e}}^{(1)}$ Figure 2:  Computational
domain

We consider three cases: (1) $Q_{2}$ is a solid obstacle, (2) $Q_{2}$ is a
bubble, and (3) $Q_{2}$ is another fluid.

For case (1), we find $({\bf{u}}_{1},p_{1})\in{\bf{V}}_{1}\times P_{1}$, such
that

$\displaystyle\left\\{\begin{aligned}
(e({\bf{u}}_{1}),e({\bf{v}}))-(p_{1},\text{div}{\bf{v}})&=({\bf{e}}_{1},{\bf{v}})\quad\forall{\bf{v}}\in{\bf{V}}_{1},\\\
(q,\text{div}{\bf{u}}_{1})&=0\qquad\forall q\in P_{1},\end{aligned}\right.$
(104)

where

$\displaystyle{\bf{V}}_{1}$ $\displaystyle=\\{{\bf{v}}\in H^{1}(Q_{1})^{2}\
\mid{\bf{v}}|_{\partial Q_{2}}={\bf{0}},\text{\ ${\bf{v}}$ is
}Q\text{-periodic}\\},$ $\displaystyle P_{1}$ $\displaystyle=\\{q\in
L^{2}_{0}(Q_{1})\ \mid q=\text{div}{\bf{v}}\text{ \ for some
${\bf{v}}\in{\bf{V}}_{1}$ }\\}.$

For case (2), we find $({\bf{u}}_{2},p_{2})\in{\bf{V}}_{2}\times P_{2}$, such
that

$\displaystyle\left\\{\begin{aligned}
(e({\bf{u}}_{2}),e({\bf{v}}))-(p_{2},\text{div}{\bf{v}})&=({\bf{e}}_{1},{\bf{v}})\quad\forall{\bf{v}}\in{\bf{V}}_{2},\\\
(q,\text{div}{\bf{u}}_{2})&=0\qquad\forall q\in P_{2},\end{aligned}\right.$
(105)

where

$\displaystyle{\bf{V}}_{2}$ $\displaystyle=\\{{\bf{v}}\in H^{1}(Q_{1})^{2}\
\mid{\bf{v}}\cdot{\bf{n}}|_{\partial Q_{2}}=0,\text{\ ${\bf{v}}$ is
}Q-\text{periodic}\\},$ $\displaystyle P_{2}$ $\displaystyle=\\{q\in
L^{2}_{0}(Q_{1})\ \mid q=\text{div}{\bf{v}}\text{ \ for some
${\bf{v}}\in{\bf{V}}_{2}$ }\\}.$

For case (3), we set $\mu_{1}=1$ and $\mu_{2}=\mu$. We find
$({\bf{u}}_{3},p_{3})\in{\bf{V}}_{3}\times P_{3}$, such that

$\displaystyle\left\\{\begin{aligned} (\mu
e({\bf{u}}_{3}),e({\bf{v}}))-(p_{2},\text{div}{\bf{v}})&=({\bf{e}}_{1},{\bf{v}})\quad\forall{\bf{v}}\in{\bf{V}}_{3},\\\
(q,\text{div}{\bf{u}}_{2})&=0\qquad\forall q\in P_{3},\end{aligned}\right.$
(106)

where

$\displaystyle{\bf{V}}_{3}$ $\displaystyle=\\{{\bf{v}}\in H^{1}(Q)^{2}\
\mid{\bf{v}}\cdot{\bf{n}}|_{\partial Q_{2}}=0,\text{\ ${\bf{v}}$ is
}Q-\text{periodic}\\},$ $\displaystyle P_{3}$ $\displaystyle=\\{q\in
L^{2}_{0}(Q)\ \mid q=\text{div}{\bf{v}}\text{ \ for some
${\bf{v}}\in{\bf{V}}_{3}$ }\\}.$

The computation is done on square grids. The first level grid consists of 12
squares, for the first two cases. Each square is subdivided into 4 sub-squares
to get the next level grid, $\mathcal{T}_{h}=\\{T\\}$. We use the
$Q_{5,4}^{1,0}\times Q_{4,5}^{0,1}$ velocity finite element space with the
$Q_{4,4}^{0,0}$ pressure finite element space. Here $Q_{5,4}^{1,0}$ means the
space of polynomials of degree at most 5 in $y_{1}$ and of degree at most $4$
in $y_{2}$ which is $C^{1}$ in $y_{1}$-direction and $C^{0}$ in
$y_{2}$-direction. That is,

$\displaystyle Q_{5,4}^{1,0}$
$\displaystyle=\Big{\\{}u_{1}|_{T}=\sum_{i=0}^{5}\sum_{j=0}^{4}c_{ij}y_{1}^{i}y_{2}^{j}\
\Big{|}\ u_{1}\hbox{ and }\partial_{y_{1}}u_{1}\in C^{0}(Q_{1}),\hbox{and
}Q\hbox{-periodic}\Big{\\}},$ $\displaystyle Q_{4,4}^{0,0}$
$\displaystyle=\Big{\\{}p|_{T}=\sum_{i=0}^{4}\sum_{j=0}^{4}c_{ij}y_{1}^{i}y_{2}^{j}\
\Big{|}\ p\in C^{0}_{0}(Q),\hbox{and }Q\hbox{-periodic}\Big{\\}}.$

We note that $\text{div}(Q_{5,4}^{1,0}\times Q_{4,5}^{0,1})=Q_{4,4}^{0,0}$.
Therefore, the finite element velocity is also pointwise divergence-free. We
plot the velocity field of these two problems in Figure 3. We can see the
magnitude of the latter is much bigger, as the resistance from a slippery
bubble is much less.

Figure 3:  The velocity field ${\bf{u}}^{1}$ for a solid obstacle $Q_{2}$
(104), and for a slippery bubble $Q_{2}$ (105).

In Figures 4 and 5, we plot the two velocity fields of two-fluid flow (106)
for two viscosity coefficients $\mu_{2}$. When $\mu_{2}$ is big, the sticky
inner fluid flows less and drags the outer fluid near the interface. When
$\mu_{2}$ approaches infinity, the inner fluid stops and it posts a zero
Dirichlet boundary condition for on tangential velocity of the outer fluid at
the inner boundary $\tilde{\Gamma}=\partial Q_{2}$. The model of a solid
obstacle (104) is a limit case of the model of two-fluid (106) when
$\mu_{2}\to\infty$. We can compare the left chart of Figure 3 and the left
chart of Figure 4.

Figure 4:  The velocity field ${\bf{u}}_{3}$ for two-fluid flow (106) with
$\mu_{2}=10^{2}$ on $Q$ (left), on $Q_{2}$ (right, scaled by 200).

When $\mu_{2}$ approaches zero, the inner fluid flow freely which produces
little drag on the outer fluid. In theory, the force inside fluid $Q_{2}$ may
even push outer fluid somewhat. But due to the zero outflow boundary condition
on the velocity at $\partial Q_{2}$, such a force would be balanced by its
left portion and right portion of an edge of $\partial Q_{2}$. It is
equivalent to zero tangential stress boundary on the outer flow. That is,
model of a slippery bubble (105) is a limit model of two-fluid (106) with
$\mu_{2}\to 0$. We may compare the right chart of Figure 3 and the left chart
of Figure 5.

Figure 5:  The velocity field ${\bf{u}}_{3}$ for two-fluid flow (106) with
$\mu_{2}=10^{-2}$ on $Q$ (left), on $Q_{2}$ (right, scaled by 2).

The homogenized permeability tensor
$\mbox{\boldmath$K$}=\begin{pmatrix}k_{11}&k_{12}\\\
k_{21}&k_{22}\end{pmatrix}$ is computed by

$\displaystyle k_{11}$ $\displaystyle=\frac{1}{|Q|}\int_{Q\setminus
Q_{2}}{\bf{u}}_{1}\cdot{\bf{e}}_{1}d{\bf{y}}\quad\text{ for }\eqref{rough},$
(107) $\displaystyle k_{11}$ $\displaystyle=\frac{1}{|Q|}\int_{Q\setminus
Q_{2}}{\bf{u}}_{2}\cdot{\bf{e}}_{1}d{\bf{y}}\quad\text{ for
}\eqref{slippery},$ (108) $\displaystyle k_{11}$
$\displaystyle=\frac{1}{|Q|}\int_{Q}{\bf{u}}_{3}\cdot{\bf{e}}_{1}d{\bf{y}}\quad\text{
for }\eqref{porous}.$ (109)

Due to the symmetry, in all our examples we have $k_{11}=k_{22}$ and
$k_{12}=k_{21}=0.$

Table 1: Computed permeability $k_{11}$ by (107)-(109).  level | (104) | (106) | (105)
---|---|---|---
|  | $\mu_{2}=10^{4}$ | $\mu_{2}=1$ | $\mu_{2}=10^{-4}$ |
1 | 0.0105 | 0.0105 | 0.0122 | 0.0140 | 0.0140
2 | 0.0119 | 0.0119 | 0.0144 | 0.0181 | 0.0181
3 | 0.0125 | 0.0125 | 0.0154 | 0.0209 | 0.0209
4 | 0.0128 | 0.0128 | 0.0159 | 0.0228 | 0.0228
5 | 0.0129 | 0.0129 | 0.0161 | 0.0240 | 0.0240

To verify the convergence results stated in (62) and (82), we solve the two-
fluid problem (106) with $\mu_{1}=1$ and $\mu_{2}=10^{-4},\,1,\,10^{4}$. In
Table 1, this model is between the two ‘limiting’ models (104) and (105).

To see how viscosity $\mu_{2}$ influences the flow, we plot
$({\bf{u}}_{3})_{1}$ in Figure 6 for two different $\mu_{2}$ with $\mu_{1}=1$.

Figure 6:  The first component of velocity ${\bf{u}}_{3}$, from (106), for
$\mu=10^{2}$ and $\mu_{2}=10^{-2}$.

Though the magnitude of $({\bf{u}}_{3})_{1}$ is way larger than that of
$({\bf{u}}_{3})_{2}$, their corresponding stress are about the same size. In
Figure 7, we plot them for a comparison. We plot the stress intensity
$|\nabla{\bf{u}}^{1}|$ in Figure 8.

Figure 7:  The stress $\nabla(u_{3})_{1}$, and $\nabla(u_{3})_{2}$ for (106)
with $\mu_{2}=10^{2}$.  Figure 8:  The stress intensity
$|e(({\bf{u}}_{3})_{1})|$ in (106) with $\mu_{2}=10^{2}$, $\mu_{2}=1$,
$\mu_{2}=10^{-2}$.

Finally we compute the energy of the two-fluid flow,

$\displaystyle E(Q_{2})$
$\displaystyle=\int_{Q_{2}}\mu({\bf{y}})e({\bf{u}}_{3}):e({\bf{u}}_{3})d{\bf{y}},$
(110) $\displaystyle E(Q)$
$\displaystyle=\int_{Q}\mu({\bf{y}})e({\bf{u}}_{3}):e({\bf{u}}_{3})d{\bf{y}}.$
(111)

The homogenized permeability can also be computed by the energy,

$\displaystyle k_{ij}$
$\displaystyle=\frac{1}{|Q|}\int_{Q}\mu({\bf{y}})e({\bf{u}}_{3}):e({\bf{u}}_{3})d{\bf{y}}.$
(112)

In Table 2, we demonstrate the equivalence of these two definitions for
$k_{11}$.

Table 2: Computed permeability $k_{11}$ both ways and energy.  level |  $k_{11}$ (109) |  $k_{11}$ (112) |  $E(Q_{2})$ (110) |  $E(Q_{2})/E(Q)$ (111)
---|---|---|---|---
| For ${\bf{u}}_{3}$ in (106) with $\mu_{2}=10^{2}$
1 | 0.107E-01 | 0.107E-01 | 0.952E-04 | 0.888E-02
2 | 0.121E-01 | 0.121E-01 | 0.811E-04 | 0.670E-02
3 | 0.127E-01 | 0.126E-01 | 0.676E-04 | 0.534E-02
4 | 0.129E-01 | 0.129E-01 | 0.604E-04 | 0.468E-02
5 | 0.130E-01 | 0.130E-01 | 0.570E-04 | 0.438E-02
| For ${\bf{u}}_{3}$ in (106) with $\mu_{2}=1$
1 | 0.122E-01 | 0.122E-01 | 0.752E-03 | 0.615E-01
2 | 0.144E-01 | 0.144E-01 | 0.135E-02 | 0.936E-01
3 | 0.154E-01 | 0.154E-01 | 0.175E-02 | 0.113E+00
4 | 0.159E-01 | 0.159E-01 | 0.200E-02 | 0.125E+00
5 | 0.162E-01 | 0.162E-01 | 0.215E-02 | 0.132E+00
| For ${\bf{u}}_{3}$ in (106) with $\mu_{2}=10^{-4}$
1 | 0.140E-01 | 0.140E-01 | 0.432E-06 | 0.308E-04
2 | 0.181E-01 | 0.181E-01 | 0.109E-05 | 0.602E-04
3 | 0.209E-01 | 0.209E-01 | 0.183E-05 | 0.875E-04
4 | 0.228E-01 | 0.228E-01 | 0.254E-05 | 0.111E-03
5 | 0.240E-01 | 0.240E-01 | 0.316E-05 | 0.131E-03

##  5 Conclusion and future work

In this paper, we show that the permeability of a porous material [40] and
that of a bubbly fluid [29] are limiting cases of the complexified version of
the two-fluid models posed in [29]. We assume the viscosity of the inclusion
fluid is $z\mu_{1}$ and the viscosity of the hosting fluid is $\mu_{1}$,
$z\in\mathbb{C}$. The proof is carried out by construction of solutions for
large $|z|$ and small $|z|$ by an iteration process similar with the one used
in [16, 21] and analytic continuation. Moreover, we also show that for a fixed
microstructure, the permeabilities of these three cases share the same
integral representation formula (IRF) (99) with different values of $s$, as
long as the ’contrast parameter’ $s:=\frac{1}{z-1}$ is not in the interval
$[-\frac{2E_{2}^{2}}{1+2E_{2}^{2}},-\frac{1}{1+2E_{1}^{2}}]$, where the
constants $E_{1}$ and $E_{2}$ are the extension constants that depend on the
geometry of $Q_{1}$, $Q_{2}$ and $Q$. For the mixture with bubbles, $s=-1$ and
thus

$K^{(B)}=\mbox{\boldmath$K$}^{(D)}+\int_{\frac{1}{1+2E_{1}^{2}}}^{\frac{2E_{2}^{2}}{1+2E_{2}^{2}}}\frac{1}{1-t}d\boldsymbol{\sigma}(t)$
(113)

Also, we note that the matrix-valued measure in (92) has a Dirac measure
sitting at $\lambda=0$ with strength equal to $\mbox{\boldmath$K$}^{(D)}$. The
permeability $\mbox{\boldmath$K$}^{(D)}$ is related to the measure in the
sense that the zero-th moment of the measure is equal to
$\mbox{\boldmath$K$}(z=1)-\mbox{\boldmath$K$}^{(D)}$.

Clearly, the positive matrix-valued measure $d\boldsymbol{\sigma}$ is
independent of $s$ and it characterizes how the geometry influences the
permeability. We have shown that this measure is related to the projection
measure of the self-adjoint operator $\Gamma_{\chi}$ and its moments can be
computed by equation (102).

Because the IRF is valid for most of $s$ on the complex plane, the IRF will be
useful in the study of two-fluid mixture with complex viscosities such as
dehomogenization for these fluid. Also, the integration limits in the IRF
should imply bounds on the permeability tensors. We will explore the results
of this paper in these direction in the future.

Acknowledgement The work of CB and MYO was partially sponsored by the US
National Science foundation via grants NSF-DMS-1413039 and NSF-DMS-1821857.

## References

  * [1] Allaire, G. Homogenization of the Stokes flow in a connected porous medium.  Asymptotic Analysis 2, 3 (1989), 203–222.
  * [2] Allaire, G. Continuity of the darcy’s law in the low-volume fraction limit.  Annali della Scuola Normale Superiore di Pisa-Classe di Scienze 18, 4 (1991), 475–499.
  * [3] Allaire, G. Homogenization and two-scale convergence.  SIAM Journal on Mathematical Analysis 23, 6 (1992), 1482–1518.
  * [4] Allaire, G. Homogenization in porous media.  CEA-EDF-INRIA school on homogenization, 2010.
  * [5] Auriault, J., Borne, L., and Chambon, R. Dynamics of porous saturated media, checking of the generalized law of Darcy.  The Journal of the Acoustical Society of America 77 (1985), 1641\.
  * [6] Avellaneda, M., and Majda, A. J. Stieltjes integral representation and effective diffusivity bounds for turbulent transport.  Physical review letters 62, 7 (1989), 753.
  * [7] Avellaneda, M., and Majda, A. J. An integral representation and bounds on the effective diffusivity in passive advection by laminar and turbulent flows.  Communications in Mathematical Physics 138, 2 (1991), 339–391.
  * [8] Avellaneda, M., and Torquato, S. Rigorous link between fluid permeability, electrical conductivity, and relaxation times for transport in porous media.  Physics of Fluids A: Fluid Dynamics 3, 11 (1991), 2529–2540.
  * [9] Beliaev, A. Y., and Kozlov, S. Darcy equation for random porous media.  Communications on pure and applied mathematics 49, 1 (1996), 1–34.
  * [10] Bensoussan, A., Lions, J.-L., and Papanicolaou, G. Asymptotic analysis for periodic structures, vol. 374.  American Mathematical Soc., 2011.
  * [11] Bergman, D. J. The dielectric constant of a composite material—a problem in classical physics.  Physics Reports 43, 9 (1978), 377–407.
  * [12] Biot, M. A. Mechanics of deformation and acoustic propagation in porous media.  Journal of applied physics 33, 4 (1962), 1482–1498.
  * [13] Brenner, S., and Scott, R. The mathematical theory of finite element methods, vol. 15.  Springer Science & Business Media, 2007.
  * [14] Brinkman, H. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles.  Flow, Turbulence and Combustion 1, 1 (1949), 27.
  * [15] Bruno, O. P. The effective conductivity of strongly heterogeneous composites.  Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 433, 1888 (1991), 353–381.
  * [16] Bruno, O. P., and Leo, P. H. On the stiffness of materials containing a disordered array of microscopic holes or hard inclusions.  Archive for rational mechanics and analysis 121, 4 (1993), 303–338.
  * [17] Cioranescu, D., Donato, P., and Ene, H. I. Homogenization of the Stokes problem with non-homogeneous slip boundary conditions.  Mathematical Methods in the Applied Sciences 19, 11 (1996), 857–881.
  * [18] Conca, C. On the application of the homogenization theory to a class of problems arising in fluid mechanics.  J. Math. Pures Appl 64, 1 (1985), 31–75.
  * [19] Darcy, H. P. G. Les Fontaines publiques de la ville de Dijon. Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d’eau, etc.  V. Dalamont, 1856.
  * [20] Dyukarev, Y., and Katsnelson, V. Multiplicative and additive classes of Stieltjes analytic matrix valued functions, and interpolation problems associated with them.  American Mathematical Society Translations 131 (1986), 55–70.
  * [21] Golden, K., and Papanicolaou, G. Bounds for effective parameters of heterogeneous media by analytic continuation.  Communications in Mathematical Physics 90, 4 (1983), 473–491.
  * [22] Johnson, D. L., Koplik, J., and Dashen, R. Theory of dynamic permeability and tortuosity in fluid-saturated porous media.  Journal of fluid mechanics 176 (1987), 379–402.
  * [23] Kantor, Y., and Bergman, D. J. Elastostatic resonances—a new approach to the calculation of the effective elastic constants of composites.  Journal of the Mechanics and Physics of Solids 30, 5 (1982), 355–376.
  * [24] Kato, T., Mitrea, M., Ponce, G., and Taylor, M. Extension and representation of divergence-free vector fields on bounded domains.  Mathematical Research Letters 7, 5 (2000), 643–650.
  * [25] Keller, J. B. Darcy’s law for flow in porous media and the two-space method.  Tech. rep., STANFORD UNIV CA, 1980.
  * [26] Lévy, T. Fluid flow through an array of fixed particles.  International Journal of Engineering Science 21, 1 (1983), 11–23.
  * [27] Limaye, B. V. Banach space-valued analytic functions.  In Spectral perturbation and approximation with numerical experiements (Canberra AUS, 1987), Centre for Mathematical Analysis, The Australian National University, pp. 44–60.
  * [28] Lions, J.-L. Some methods in the mathematical analysis of systems and their control(book).  Beijing, Science Press (1981).
  * [29] Lipton, R., and Avellaneda, M. Darcy’s law for slow viscous flow past a stationary array of bubbles.  Proceedings of the Royal Society of Edinburgh Section A: Mathematics 114, 1–2 (1990), 71–79.
  * [30] Lundgren, T. S. Slow flow through stationary random beds and suspensions of spheres.  Journal of Fluid Mechanics 51, 2 (1972), 273–299.
  * [31] Milton, G. W. The theory of composites.  The Theory of Composites, by Graeme W. Milton, pp. 748. ISBN 0521781256\. Cambridge, UK: Cambridge University Press, May 2002. (2002), 748\.
  * [32] Neuman, S. P. Theoretical derivation of darcy’s law.  Acta Mechanica 25, 3-4 (1977), 153–170.
  * [33] Nguetseng, G. A general convergence result for a functional related to the theory of homogenization.  SIAM Journal on Mathematical Analysis 20, 3 (1989), 608–623.
  * [34] O.A. Oleinik, A.S. Shamaev, G. Y. Mathematical Problems in Elasticity and Homogenization., 1 ed., vol. 26 of Studies in mathematics and its applications.  North Holland, 1992.
  * [35] Ou, M.-J. Y. On reconstruction of dynamic permeability and tortuosity from data at distinct frequencies.  Inverse Problems 30, 9 (2014), 095002.
  * [36] Poreh, M., and Elata, C. An Analytical Derivation of Darcy Law.  Publication (Tekhniyon - Makhon tekhnologi le-Yisra’el. ha-Fakultah le-handasah ezrahit). Technion-I.I.T., Faculty of Civil Engineering, 1965.
  * [37] Saffman, P. G. On the boundary condition at the surface of a porous medium.  Studies in applied mathematics 50, 2 (1971), 93–101.
  * [38] Sanchez-Palencia, E. Non-homogeneous media and vibration theory.  Lecture notes in physics. Springer, 1980.
  * [39] Tam, C. K. The drag on a cloud of spherical particles in low reynolds number flow.  Journal of Fluid Mechanics 38, 3 (1969), 537–546.
  * [40] Tartar, L. Incompressible fluid flow in a porous medium-convergence of the homogenization process.  Appendix of Non-homogeneous media and vibration theory (1980).
  * [41] Tice, I. From stokes flow to darcy’s law.  CNA Working Group on Homogenization, 2014.