1. business and industrial
2. chemicals industry
3. plastics and polymers

Downloaded from http://rspa.royalsocietypublishing.org/ on December 3, 2014
A di® use interface model of two-phase ° ow in
porous media
By Ve n k a t G a n e s a ny a n d H o w a rd B r e n n e r
Department of Chemical Engineering, Massachusetts Institute of Technology,
Cambridge, MA 02139-4307, USA
Received 20 April 1998; revised 7 April 1999; accepted 26 April 1999
A single-phase, two-component mixture Darcyscale model characterized by a di¬use
interface is proposed as a rational alternative to the conventional singular interface Darcyscale empirical model currently employed for analysing two-phase ®ows of
immiscible ®uids through porous media. The proposed scheme possesses conceptual
as well as computational and analytical advantages. On the conceptual side, we are
able to clearly de­ ne the macroscale concept of capillary pressure, which in contemporary literature is incorrectly confounded with the well-known microscale, pore-level
Laplace boundary condition at the curved interfaces between the immiscible phases.
In this same context we o¬er insights into fundamental issues that have a°icted mixture theories involving `interpenetrating continua’. This is accomplished in part by
clarifying the distinction between whether the continua being referred to constitute
phases (consisting of `immiscible’ multiphase continua) or species (present in an inhomogeneous single-phase continuum). In particular, we show that the issue devolves
upon the scale at which the phenomenon is viewed. On the computational and analytic sides our scheme o¬ers the advantages of dealing with continuous ­ elds rather
than with discontinuous ­ elds, the latter necessitated in contemporary literature by
the existence of singular (mobile) phase boundaries, namely interfaces. Also on the
analytical side, our scheme permits a formal physicomathematical transition from the
di¬use microscale to the coarser singular-surface scale view, by using singular perturbation techniques to achieve the requisite change in scale. Within this framework
we formulate, in a rigorous manner, de­ nitions of macroscale quantities, following
which we identify the pertinent phase-speci­ c Darcy’s laws. Moreover, one can, in
principle, calculate the phenomenological coe¯ cients appearing therein in terms of
quadratures of the prescribed microscale data. Our macroscale framework treats the
Darcyscale continuum as a multicomponent mixture (referred to as the `di¬use Darcyscale’ view) rather than adopting the more commonly employed coarse-grained
picture embodied in interpenetrating continua models or singular-surface Darcyscale
models. This ­ ne-scale viewpoint is directed towards the foundational aspects of
two-phase ®ows at both the micro- and macroscales|especially with regard to the
existence and interpretation of phase-speci­ c Darcy’s laws, including capillary pressure as a macroscale ­ eld variable. Eventually, we implement this framework with
a simpli­ ed linear example to substantiate our thesis. Our study embodies several
distinct investigations, logically intertwined by the common objective of erecting
y Present address: Department of Chemical Engineering, University of California, Santa Barbara, CA
93106, USA.
Proc. R. Soc. Lond. A (2000) 456, 731{803
® c 2000 The Royal Society
731
732
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
rational foundations for describing and quantifying two-phase ®ows through porous
media.
Keywords: two-phase ° ows in p orous media; di® use interface mo del;
Darcy’s laws; capillary pressur e hysteresis
1. Introduction
(a) Flows in porous media
A plethora of natural phenomena and engineered processes involve the ®ow of immiscible ®uids through porous media. For example, the movement of water and air
through soils, relevant to the ­ eld of hydrology (Philip 1970), has long been studied
from both theoretical and experimental perspectives beginning with the pioneering
work of Buckingham (1907). In the ­ eld of petroleum engineering, crude oil is often
displaced from the interstices of the rock by injecting water. Examples of bi- and
multi-phase ®ows relevant to industrial applications are far too numerous and well
known to require review here. An overall introduction relevant to applications is
provided by Marle (1981).
Theoretical analysis of biphasic ®ows through a porous matrix, which constitutes
the main focus of this paper, is extremely complicated owing, inter alia, to the complex geometric con­ guration of the pore space. Realistic porous-medium geometries
constitute a virtually intractable con­ guration, even for single-phase ®ows. In fact,
complete experimental characterization of the skeletal structure of an actual porous
medium is itself an extremely di¯ cult geometric proposition (Adler 1992).
Given these circumstances, and motivated by both hydrological and petroleum
engineering applications, it has proved expedient heretofore to formulate theoretical
analyses of ®uid ®ow phenomena through porous media at an observable scale (the
Darcyscale or macroscale) through the use of averaged pore-scale (i.e. microscale)
­ elds. Fundamental to this notion is the assumption that to each macroscale `point’
there corresponds a microscale volume element encompassing a representative sample
of the local (instantaneous or time-averaged) contents of the porous medium in the
neighbourhood of that point (Bear & Bachmat 1984). This dual scale `micro{macro’
modelling is justi­ ed under the assumption that there exists a wide disparity in
linear dimensions between the pore and observable scales, whence microscale lengths
appear as in­ nitesimals at the macroscale. The porous rock formations present in
naturally occurring oil recovery wells exemplify such disparate length-scales. While
on the one hand the pore level is characterized by length-scales of the orders of
millimetres, with pore sizes of the order of microns, on the other hand the Darcyscale
or macroscale is characterized by changes occurring on the length-scales of metres|
thereby facilitating such a `macro’ approach.
Under such a micro{macro modelling approach, and for single-phase ®ows, Darcy’s
law governs the slow, monophasic ®ow of a Newtonian ®uid through a porous medium
under the animating e¬ects of gravity and/or externally imposed pressure gradients (Neumann 1977). The essence of Darcy’s law is the relationship between the
macroscopic pressure gradient rp and the macroscopic (super­ cial) velocity vector
v (notation for the macroscopic quantities is clari­ ed in a later section):
K
(1.1)
(rp ¡ » g);
v=¡
·
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
733
where · and » , respectively, denote the viscosity and the density of the interstitial
®uid, g is the gravity vector, and K the permeability tensor. The latter phenomenological coe¯ cient encapsulates all of the macroscopically unobservable microscale
phenomena. For a macroscopically inhomogeneous ®ow the above relationship is
conventionally assumed to hold at each macroscale point R in the porous medium.
In modelling multiphase ®ows, Darcy’s law has been interpreted as applying to
each phase separately, this multiphase analogue of Darcy’s law being proposed on an
empirical basis (Marle 1981). Explicitly, with v ¬ the super­ cial or seepage velocity
vector of phase ¬ relative to the ­ xed pores, the two-phase form of Darcy’s law is
taken to be
1
(¬ = 1; 2);
(1.2)
KKr¬ (rp¬ ¡ » ¬ g)
v¬ = ¡
· ¬
where Kr¬ represents the relative permeability of phase ¬ . The latter is de­ ned as
Kr¬ = K¬ =K;
with K¬ the permeability of the porous medium to phase ¬ , and K a norm of the
single-phase permeability of the porous medium (assumed to be a strictly geometrical
quantity, possessing the same value for each phase) (cf. (1.1)); » ¬ denotes the macroscopic density of phase ¬ (i.e. the mass of phase ¬ per unit of macroscopic interstitial
volume). The presence of more than one phase is conventionally quanti­ ed through
the saturation s¬ (s1 + s2 = 1), the latter denoting the macroscopic volume fraction
of phase ¬ . Kr¬ is assumed to be dependent only on the saturations. Furthermore,
p1 ¡ p2 = Pc [s1 ] de­ nes the capillary pressure, assumed to be a function only of the
saturation.
(b) Review and critique of prior works
The three major aspects of multiphase ®ows towards which prior researches have
primarily been directed are as follows.
(i) Many studies (Leverett 1941; Marle 1981; deGennes 1983) have focused on
the dynamics of displacement fronts in two-phase systems, modelling the resulting
phenomena using the multiphase form (1.2) of Darcy’s law. Conservation equations
for the respective saturations s1 and s2 , namely
@s2
@s1
(1.3)
+ r (s2 v2 ) = 0;
+ r (s1 v1 ) = 0;
@t
@t
supplement (1.2)|with closure of the system of equations e¬ected by the relationships Pc = Pc [s1 ] and Kr¬ = Kr¬ [s1 ], the explicit functional relationship being
assumed known from experiments.
(ii) Percolation theories (Broadbent & Hammersly 1957) have been used as simplistic models of the complex pore-level ®ow phenomena (Larson et al . 1977). These
theories, which concentrate on universal behaviour invariant to the quantitative geometrical details of the system (as well as to the physical properties of the ®owing
®uids), are usually used to model transitions occurring between gross types of system behaviour. However, the exact relationship between such discrete models and
the actual ®ow phenomena taking place in porous media has not been rationally
established (Adler & Brenner 1988).
Proc. R. Soc. Lond. A (2000)
734
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
(iii) A third aspect of research, pursued primarily by continuum mechanicians
(Marle 1977; Gray & Hassanizadeh 1990), attempts to establish the credibility of the
macroscale equations (1.2) to (1.3), by averaging the `pore-level’ Navier{Stokes and
continuity equations underlying them, albeit generally by di¬erent schemes speci­ c
to each investigator.
Our own contribution falls into the third category, although from a unique perspective which involves viewing the ®uid{®uid interfaces in a porous medium as di¬use,
continuous transition regions, rather than as singular surfaces separating immiscible
phases.y
It has long been argued (Hassanizadeh & Gray 1993) that empirical approaches to
modelling multiphase ®ows based on assumed analogies with single-phase ®ows might
be invalidated either by theory, experiments, or both. We expound below on what
we believe to be some of the accepted shortcomings of such an empirical extension.
Eventually, this critique provides the motivation for our own work.
(a) Whereas Darcy’s law for single-phase ®ow can be proved rigorously (Neumann
1977; Keller 1985) beginning with a microscale model of the pertinent ®ow phenomena, the rational extension of this approach to multiphase ®ows has not yet been
achieved. Even the choices for the appropriate variables to describe such multiphase
®ows (saturation s1 , capillary pressure Pc ) have yet to rigorously justi­ ed (Adler &
Brenner 1988). Furthermore, the fundamental assumption that the seepage velocity
of each phase should depend functionally only upon the pressure gradient existing
in that phase has been called into question (deGennes 1983; Kalaydjian 1993).
(b) A natural fallout arising from the traditional choice of independent variables
is the presence of the capillary pressure Pc , which is usually branded as an `experimental’ quantity (Marle 1981; deGennes 1983). The physical interpretation of this
variable has proved a major deterrent to many theoretical analyses of multiphase
®ows for manifold reasons|the primary obstacle being the commonly observed hysteretic phenomena occurring in experimental capillary pressure versus saturation
curves (­ gure 1), depending on whether drainage (displacement of the wetting ®uid)
or imbibition (displacement of the non-wetting ®uid) is involved. Though normally
represented by two distinct curves, it is not universally recognized that the capillary pressure may also take on any value intermediate between the two. Yet to be
explained is the fundamental reason for expecting such a `phase equilibrium’ type of
relationship to exist between the macroscopic pressures in the two ®owing phases|
neither of which phases exists in a state of equilibrium, much less with one another.
In addition, the experimental capillary pressure curves are obtained by allowing the
multiphase mixture to attain static equilibrium. It is questionable whether use of the
same static data is appropriate for modelling the inherently dynamical phenomena
occurring during ®ow through porous media, though that is what is frequently done.
Our major criticism of prior approaches to multiphase ®ow modelling lies in their
failure to clearly distinguish conceptually between what constitutes a `phase’ and
what constitutes a `component’ (or `species’) (Bedford & Drumheller 1983). The
confusion is compounded by loose usage of the concept of interpenetrating continua
y It should be noted that even these past singular interface e¬orts (Marle 1977; Whitaker 1986;
Pavone 1989; Gray & Hassanizadeh 1990), aimed at rationally e¬ecting this micro{macro transition,
have generally resulted in only partial (and in some cases even contradictory) answers.
Proc. R. Soc. Lond. A (2000)
Pc
1
2
735
irreducible saturation in
wetting fluid
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
residual saturation in
non-wetting fluid
0
SNM
1
Figure 1. Capillary pressure versus saturation relationship. S N M represents the saturation of
the non-wetting ° uid (adapted from Marle (1981)).
(Gray & Hassanizadeh 1990). Such interpenetrating continuum (IPC)-type theories of immiscible mixtures were originally formulated in the context of ®uid{solid
particle mixtures (Nigmatulin 1979), wherein such a viewpoint seemed intuitively
appropriate, especially when no material was exchanged between the solid and the
®uid phases. Among other distinguishing features of this special class of mixtures is
the fact that which of the two phases constitutes the continuous phase is never in
doubt. However, mixture theories for immiscible ®uid{®uid systems possess manifestly contrasting features (Bedford & Drumheller 1983; Dobran 1984), primarily as
a result of the interfacial boundary conditions as well as of the question of which
of the two phases is the continuous phase (not even addressing the issue of `bicontinuous’ mixtures (Davis 1995), where neither phase can be uniquely classi­ ed as
constituting the continuous or discontinuous phase). Recently, the interfacial boundary conditions were themselves also averaged (Gray & Hassanizadeh 1990), thereby
creating additional ­ eld equations involving new physical variables. The additional
`macroscale’ interfacial equations arising from this interfacial averaging scheme then
need to be solved conjointly with the above singular Darcyscale equations. However,
the physical interpretation to be ascribed to such additional interfacial ­ eld variables
seems obscure in the context of the IPC picture. Moreover, as indicated below, the
existence of two ad hoc tenets underlying the philosophy of mixture theories has
diminished the practical utility of its predictions.
(a) The spirit of the mixture theories is to accept the macroscale equations (obtained by suitable averaging of the microscale equations) as fundamental, and to
propose constitutive relations for the various macroscale ­ elds appearing therein as
functions of a select set of independent variables (chosen so as to satisfy suitable
general criteria) (Gray & Hassanizadeh 1990). The sentiment is to compare the
resulting mathematical predictions (after complete quantitative calculations of the
requisite macroscale ­ eld variables) of such a model with experimental observations
so as to authenticate the model|a step rarely achieved fully in practice owing to
the complexity of the governing equations as well as of the experimental di¯ culties
in e¬ecting the requisite measurements.
Proc. R. Soc. Lond. A (2000)
736
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
(b) The manner in which mixture theories treat constitutive equations as fundamental entities has eventually obscured the physical interpretation of the singular Darcyscale ­ elds in terms of their precursor singular-interface microscale ­ elds.
For instance, Hassanizadeh & Gray (1993) state that: `The pressures present in the
capillary pressure relationships are not the averages of the (respective) microscopic
pressures but instead are to be determined from a macroscopic constitutive equation
for the Helmholtz free energy A.’ However, the precise relationship of the macroscale
free energy A to its microscale counterpart A is never clari­ ed.
(c) Outline and summary of article
What follows is not a comprehensive programme aimed at reformulating the theory of mixtures as a whole based on the di¬use interface model. Rather, attention is
focused explicitly on the more limited goal of modelling slow biphasic ®ows through
porous media. Nevertheless, some of the concepts advanced are believed to possess
relevance within the more general context of mixture theory. The philosophy used in
our modelling approach di¬ers from that underlying conventional theories of multiphase ®ow through porous media in the following ways.
(a) The central notions of immiscibility, species, and phases at both the microand macro-levels are clari­ ed, invoking a viewpoint di¬erent from that employed in
traditional (singular Darcy) approaches to multiphase ®uid ®ow phenomena. Simultaneously, those aspects involving interfacial balance conditions are accounted for
in a consistent manner. Purely formal axiomatic continuum-mechanical theories are
eschewed in favour of the physical aspects of the phenomena. For example, rational
mechanical treatments of mixtures establish appropriate thermodynamic identities
through use of the equilibrium version of the Coleman{Noll inequality (cf. Gray &
Hassanizadeh (1990) and the references cited therein). In contrast, here it is assumed
that those classical intensive thermodynamic identities which are locally valid pointwise for inhomogeneous microscale continua hold equally well at the macroscale
(di¬use Darcyscale). Given the simplistic nature of the description a¬orded by our
treatment, we believe that such a hypothesis would be borne out from the Coleman{
Noll inequality too.
(b) Our identi­ cation of appropriate macro­ elds is based upon their rigorous physical, scale-invariant de­ nitions rather than upon simple ad hoc volume averaging
of the corresponding micro­ elds. For instance (see x 3 e), the continuum-level Darcyscale mass-average seepage velocity vector v at a given point R, say, is de¯ned
as that ­ eld which when scalarly multiplied by an arbitrarily directed Darcyscale
surface element ds situated at R together with the Darcyscale mass density » gives
the physical mass rate of ®ow,
_ = » ds v;
dm
(1.4)
at which ®uid crosses ds in the direction of n (the normal to ds). This continuummechanical de­ nition of v contrasts with the usual de­ nition of v as being the
volume-average of the micro­ eld v taken over some representative volume element
V centred at R:
1
v dV:
v=
V V
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
737
Similarly, the Darcyscale continuum stress ­ eld P is de­ ned so as to satisfy the
physical, Cauchy, continuum-mechanical de­ nition
dF = ds P ;
(1.5)
where (with the usual convention) dF is the vector force exerted across ds by the
®uid situated on one side of the surface upon the ®uid lying on the opposite side. At
appropriate points in the subsequent text, similar physically based macroscale de­ nitions will be enunciated relevant to the other pertinent ­ elds there being introduced.
Such considerations are of prime importance since they relate directly to the manner in which an experimental probe, with its appropriately sized aperture, properly
registers the local values of the dynamical variables it purports to measure. Conventional theories, which use a mathematical volume averaging of the microscale ­ elds
to de­ ne their respective macroscale counterparts, generally fail to provide such a
physical interpretation of the macroscale ­ elds (Whitaker 1986, 1999).
(c) The di¬use Darcyscale ­ elds appearing in the macroscale conservation equations governing the pertinent transport phenomena are explicitly de­ ned through
their microcontinuum counterparts, as they are, after all, a manifestation of the
underlying microscale phenomena viewed on a di¬erent length-scale. In any rational
theory, operational de­ nitions of the macroscale quantities calculable in principle
(from solutions of the governing microscale equations) need to be explicitly provided
for those willing to expend the computational resources necessary to literally perform
these calculations based on these `recipes’. This philosophy is similar to that underlying the statistical mechanics of transport processes (the foundations of which were
legitimatized and rendered rigorous by Kirkwood (1967)). In particular, that scheme
provides microscopic, molecularly based expressions for the magnitudes of macroscopic transport coe¯ cients like viscosity, etc., despite the fact that the constitutive
equations themselves were already known on the basis of macroscopic experiments.
In our analogy we are in some sense, albeit imperfectly, likening: (i) the discrete
molecular, microscale view with the discrete nature of the particle- (or equivalent
pore-) level view of the porous medium; (ii) the continuum-level manifestation of
this discrete molecular view with our continuum Darcyscale view.
Important conceptual di¬erences exist between conventional two-component systems and two-phase systems, especially with respect to the temporal con­ gurational
evolution of the pertinent ­ elds, where the concept of a `micro{macro’ time-scale
needs to be identi­ ed and quanti­ ed. Our contribution appears to be the ­ rst wherein
a rigorous foundational basis is indicated for a time-averaging scheme (see also Ishii
(1975)). Eventually, rigorous, unambiguous, physically based de­ nitions are adopted
for spatially and time-averaged di¬use Darcyscale ­ elds in terms of their microscale
­ eld counterparts. Subsequently, these macroscale de­ nitions are used to derive the
Darcyscale ­ eld conservation equations from their pore-level precursors.
Within the proposed framework of a multicomponent mixture description of ®ow
at the Darcyscale, the concept of phase-speci¯c pressures is identi­ ed with speciesspeci¯c pressures, and a relationship for the capillary pressure derived in terms
of appropriately de­ ned microcontinuum ­ elds. Such an analysis serves to clarify
the relationship between our ­ ne-scale, di¬use Darcyscale model and conventional,
coarser-scale, singular-surface Darcyscale models (see x 3 g). Discussion following the
derivation of the capillary pressure relationship in x 3 i leads to several signi­ cant
Proc. R. Soc. Lond. A (2000)
738
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
observations. As an example illustrating the application of our model in practice, in
x 3 m we use a phenomenological set of pore-level constitutive equations to show how
the previously empirical two-phase Darcy’s law formulation can be derived from ­ rst
principles via a set of linearizing assumptions. Explicitly, the di¬use two-component
Darcyscale ­ eld equations are subjected to a singular perturbation analysis, enabling
us to derive the conventional, otherwise empirical, singular-surface, Darcyscale equations commonly employed in practice. This illustrates the fact that the novel framework erected in this paper can serve a dual purpose|namely, to provide insights
into the forms of macroscopic constitutive equations as well as to derive the more
generic macroscale transport equations into which such constitutive relations are to
be embedded.
It is important to delineate the major assumptions involved in the subsequent
analysis so as to both identify the scope of our work and to spur future extensions
thereof. We concentrate on a simple geometric model of a porous medium, namely a
spatially periodic skeletal array (Brenner 1972; Brenner & Edwards 1993) (described
in x 3 c). Such a de­ nite geometrical construct, which may at ­ rst appear overly
restrictive, has been shown to serve as a good theoretical model for a number of
phenomena to date, especially single-phase ®ows through porous media (Bensoussan
et al . 1978). Deterministic models of this nature (in contrast to statistical models)
enable one to de­ ne in a completely rigorous manner all of the macroscopic ­ elds
appearing in the Darcyscale description of the several transport processes, while
simultaneously furnishing de­ nitive error estimates arising in the eventual continuum
description resulting from the spatial and temporal homogenization of the discrete
microscale system. Describing the ®ow and transport in such model porous media
leads to a consideration of what are termed `homogeneous multiphase ®ows’, de­ ned
as follows (Quintard & Whitaker 1988; Whitaker 1999).
A porous medium is homogeneous with respect to a given process and
averaging volume when the e¬ective transport coe¯ cients in the volumeaveraged equations are independent of position. If the porous medium is
not homogeneous, it is heterogeneous.
Results based on strictly spatially periodic ­ elds can also be expected to apply in
circumstances wherein a slow spatial modulation is superposed on the spatial periodicity (common models for homogenization theories (Bensoussan et al . 1978; Auriault 1987)). One might also hope that results gleaned from this modulated model
could also be used for the macroscopic modelling of displacement fronts|by arguing
that the macroscopic front is, in reality, a di¬use transition region extending over
a signi­ cant microscale distance (embodying many unit cells of the spatially periodic medium, but nevertheless still small compared with the overall length of the
porous medium in the direction of net ®ow). The phenomenon of solid{liquid contact lines (Dussan 1979), both static and dynamic, as well as issues of wettability
of the porous medium, is omitted in this initial communication. Though very relevant to practical ®ows in porous media, this area remains extremely controversial.
Its role in multiphase ®ow phenomena is still widely debated to the extent that even
the boundary-condition equations governing the microscale ®uid motions are not yet
fully established. We hope to return to these issues in future.y
y Issues pertaining to wettability can perhaps be incorporated into our work by modelling the presence
of the solid surface through an interaction potential (see Cahn 1977; Teltezke et al . 1982).
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
739
Finally, with regard to the thermodynamics of these di¬use `multiphase’ systems,
we make the crucial assumption in the subsequent development that it is permissible
to employ standard thermodynamic identities at both the micro- and macroscales.
This assumption appears reasonable in view of the multicomponent mixture model
(see also Marle 1977; Pavone 1989; Kalaydjian 1993) subsequently proposed in x 2 as
the conception underlying di¬use Darcyscale phenomena. For instance, we hypothesize that a macroscale thermodynamic chemical potential · can be de­ ned within
our di¬use Darcy framework as
·
d ef
=
@F
;
@»
(1.6)
given the existence of a macroscopic free energy density function F and a macroscopic
density » , each being derived (in a manner to be discussed) from the respective
comparable continuum microscale ­ elds F and » by well-de­ ned spatial and temporal
averages of the latter.
While a complete explanation of the panoply of multiphase ®ow phenomena accompanying ®ow through a porous medium of arbitrary con­ guration is the eventual goal
of our di¬use interface approach, such an e¬ort is overly ambitious at present. Porous
rock formations encountered in oil recovery wells are neither homogeneous nor possess a periodic microstructure. Further, the wetting properties of the liquids play
an important role in determining the microscale ®ow patterns. Moreover, ®ows in
real rock formations may be characterized by highly non-equilibrial phenomena, like
Haynes jumps, etc., representative of high Reynolds number ®ows. Our assumptions,
outlined in the previous paragraph and subsequently elaborated in the text, eschew
these complications with the aim of developing a mathematically tractable model of a
porous medium so as to uncover key insights which would hopefully prove applicable
to real ®ows through the irregular formations encountered in nature.
In essence, the paper is divided into two distinct complementary parts, represented
by xx 2 and 3. Section 2 provides a self-contained generic introduction to the di¬use
interface model without reference to its subsequent application to porous medium
®ows. Section 3 entails speci­ c applications of the generic framework erected in x 2
to situations involving ®ows through porous media.
Section 2 reviews the di¬use interface model embodying our own perspective for
describing immiscible ®ows, using what we have chosen to call the `multicomponent
mixture’ framework. This approach successfully overcomes many of the di¯ culties
outlined above that have been encountered in prior analyses. Section 2 a begins with
an introduction to the equations and boundary conditions encountered in the conventional or the singular interface model for describing two-phase ®ows. Section 2 b
reviews the spirit of the di¬use interface model, along with a repertoire of comparable conservation equations encountered within this framework. Section 2 c outlines
the constitutive forms for the various ­ elds present in the conservation equations
outlined in the previous section. Section 2 d is devoted to clarifying issues relating to
immiscibility, and to the manner in which it is supported within our di¬use interface
framework. Finally, in x 2 e we illustrate (using singular perturbation techniques)
the precise manner in which the di¬use interface equations reduce to the singularinterface multiphase ®ow equations. This section also serves simultaneously to clarify
the generic mathematical framework within which the di¬use interface model needs
to be used in order to e¬ect analytical calculations of practical ®uid-mechanical probProc. R. Soc. Lond. A (2000)
740
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
lems, namely, through a sequential, successively improvable, singular perturbation
approach.
Section 3 begins with a brief introduction to the geometrical characteristics of
the spatially periodic model of a porous medium, placing special emphasis on the
notation and de­ nitions subsequently used in the text. Sections 3 b,c consider the
appropriate microscale conservation equations applicable within the binary mixture
framework model. Using physically based de­ nitions outlined in x 3 d we execute a
rigorous transition from the di¬use microscale to the di¬use Darcyscale to derive
the Darcyscale equations in x 3 e. Concepts drawn from irreversible thermodynamics
are used in x 3 f to suggest and illustrate possible forms for the macroscale constitutive relations. Section 3 g de­ nes phase-speci­ c quantities based upon di¬use
Darcyscale ­ elds. In x 3 i we propose a rational de­ nition of capillary pressure, arising as a natural consequence of the transition from the di¬use Darcyscale to the
singular Darcyscale. Up until x 3 j the analysis completely eschews use of the constitutive equations outlined in x 2 c. In x 3 j and x 3 k we use the constitutive equations
outlined in x 2 c to provide a number of new insights into the observed hysteresis
phenomena occurring during capillary pressure versus saturation experiments. This
discussion is supplemented with an exposition of results gleaned from the constitutive equations governing the non-equilibrium conservation equations in x 3 l. Section
3 m applies singular perturbation arguments to the di¬use-interface microscale equations, enabling us to correlate our multiphase ®ow, porous medium results with more
traditional, singular interface approaches. Finally, x 4 summarizes our work, pointing
out possible future research directions. Several appendices and footnotes, elaborating
details that are glossed over in the text, complete the paper.
2. Part I. Di® use interface model
An interface has traditionally been viewed (Gibbs 1957) as a singular surface separating two immiscible bulk ®uids or phases. In equilibrium situations these bulk
phases are regarded as being separately homogeneous, although each will generally
be inhomogeneous in non-equilibrium states. Each bulk phase is separately characterized by its own material properties which, while continuous within the individual
bulk phases, may display a sharp jump in magnitude across the interface. This singular interfacial surface is itself endowed with and characterized by its own material
properties, such as interfacial tension, shear and dilatational viscosities, etc. (Scriven
1960). Traditional microscale descriptions of two-phase ®ows (in the absence of the
porous medium) involve conservation ­ eld equations written for each of the two bulk
phases, together with corresponding matching conditions imposed on the respective
values of the bulk ­ elds at the interfacial surface|the latter matching conditions
being derived by considering the pertinent transport phenomena occurring across
and within the interface (cf. (2.1) and (2.2)) entailing new, interfacial ­ elds (Scriven
1960; Edwards et al . 1991).
This two-dimensional, singular view of the interface has, however, long been recognized as representing only an asymptotic mathematical approximation of the true
physical state of a¬airs (Kirkwood & Bu¬ 1949; Tolman 1949). In reality, when
viewed on a su¯ ciently ­ ne length-scale the interfacial region between two immiscible ®uids is a highly inhomogeneous three-dimensional transition region (even at
equilibrium) over which rapid changes in continuum-mechanical ­ elds and concomiProc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
741
Table 1. Identi¯cation of the volumetric and areal ° ux density terms appearing in the generic
conservation equation (2.1) for the transport property P
(P, physical property; Ã , volumetric density; J , areal-° ux density; ¼ , volumetric rate of production; ³ , volumetric rate of supply by external sources.)
P
mass
species i
linear momentum
Ã
J
¼
³
»
0
JD
¡P
0
0
0
0
0
» F
»
i
» v
tant material properties may occur (Eliassen 1963; Goodrich 1981). In this threedimensional di¬use view of the interface, the relevant continuum-mechanical phenomenological properties vary continuously throughout the entire ®uid domain, with
some ­ elds experiencing extremely steep gradients within the interfacial region in a
direction normal to the `interface’. Examples of such rapidly varying ­ elds include
the pressure ­ eld and, when a surfactant is `adsorbed’ within the interfacial region,
the three-dimensional mass density ­ eld » i of the surfactant species i (Slattery 1964).
Although the traditional singular-surface view of an interface proves specially convenient in many practical contexts, some applications require a ­ ner-scale continuum
view of the interfacial region. Such microscale perspectives have primarily been used
for two di¬erent purposes.
(a) Theoretical investigation of transport processes involving moving and deforming interfaces (Eliassen 1963; Goodrich 1981; Mavrovouniotis & Brenner 1993): these
works concentrate on the distinction between the true, di¬use picture of the interface and its approximate representation as a singular surface. Such theories attempt
a reconciliation of the continuous, di¬use-interface view, which admits large inhomogeneities in material properties over a narrow transition region, with the more
conventional singular-surface view of the interface, which lumps the inhomogeneities
into so-called `surface-excess’ quantities. One such application, emphasized in the
work of Mavrovouniotis & Brenner (1993), deals with the distinction between material and non-material interfaces.
(b) More recently, the di¬use view of interfaces has gone beyond simply clarifying the physical meaning of surface-excess quantities. Phenomenological constitutive
equations have been proposed that enable one to work directly with the ­ ner-scale
continuous, three-dimensional, inhomogeneous ®uid. This is in lieu of working separately with the bulk ®uids lying on either side of the singular interface together with
the corresponding matching conditions imposed at that interface. These approaches,
which owe their origin to the seminal work of Cahn & Hilliard (1958), have been
revived and popularized by Caginalp (1986), Caginalp & Fife (1988), and others. The
present investigation entails one further, albeit novel, application of this approach.
(a) Singular interface equations
The singular interface approach exempli­ es the discrete nature of the bulk ®uids
or phases of a two-phase system. Such a discrete structure is inherent in the governing microscale conservation equations (Edwards et al . 1991), which can be written
Proc. R. Soc. Lond. A (2000)
742
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
separately for each of the two phases in the generic form (refer to table 1 for speci­ c
identi­ cations),
@Ã
+ r (vÃ) + r J = ¼ + ³;
@t
(2.1)
to which are appended matching conditions at the interface embodying the generic
interfacial transport equation:
¯ sÃs
¡
¯ t
u rs à s ¡
rs (Is vÃ)s + rs (Is J s ) ¡
¼s ¡
Ãs
= ¡ n [[(v ¡
u)Ã + J ]]: (2.2)
In the latter the superscript s denotes surface-excess quantities ascribed to the singular interface, whereas n represents the unit normal to the interface, Is ² I ¡ nn
the surface idemfactor, rs ² Is r the surface gradient operator, and u the velocity
of the interface (which in the general case of non-material interfaces can di¬er from
the bulk velocities v (i.e. v + ; v ¡ ) on either side of the interface); [[: : : ]] denotes the
jump in the value of the argument across the interfacial surface. Quanti­ cation of the
surface-excess mass balance is achieved through the surface-excess mass density » s ,
which is conventionally assumed to be zero. Adsorbed surfactant species at the interface correspond to non-zero values of » si (i = 1; 2; : : : ). Furthermore, ¯ s =¯ t denotes
the interfacial convected time derivative (Edwards et al . 1991; Mavrovouniotis &
Brenner 1993).
(b) Di® use interface model
Theoretical analyses which treat the interface as a region characterized by a continuous variation of physical properties owe their genesis to the phenomenological
constitutive model proposed by Cahn & Hilliard (1958). The latter can be also be
viewed as a variant of the dynamical Landau{Ginzburg theory of critical phenomena
(Hohenberg & Halperin 1977). This model allows for the coexistence of two phases by
representing the free energy of the interfacial system as a sum of two terms: (a) a free
energy function (dependent on the component concentrations) describing the miscibility behaviour; (b) a term dependent on the gradients (commonly truncated at the
second order) of the component concentrations, serving to model interfacial tension
and other capillary e¬ects (cf. x 2 c). This approach is widely used to model spinodal
decomposition and nucleation phenomena (Cahn 1965; Cahn & Hilliard 1959). While
it is generally cautioned that such models intrinsically possess a restricted range of
validity (Cahn & Hilliard 1958), being strictly applicable only near the critical point
(or the consolute point, in the case of binary liquid mixtures), great success has
nevertheless been achieved in using this model for studying non-critical behaviour in
several other contexts, including the equilibrium thermodynamics of microstructures
(Davis & Scriven 1982) as well as the dynamics of solid melts (Caginalp 1986). It
is generally believed that the approximations inherent in using this model for noncritical mixtures do not result in excessively large errors when applied to binary
liquid mixtures (Turski & Langer 1980).
In this section we outline the general form of the two-phase equations, valid for
incompressible and immiscible two-phase mixtures. Discussion of the constitutive
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
743
equations required therein, for which we use variants of the Cahn{Hilliard model, is
postponed until the next section.
It is imperative before proceeding further to imbibe the spirit of the di¬use interface model by recognizing that at this ­ ne-scale level of description there no longer
exists a discontinuous or singular surface separating two bulk phases. Rather, the
system is envisioned as being an inhomogeneous single-phase multicomponent mixture, more precisely a solution (the number of components being chosen here as two),
locally described by conventional, spatially varying composition variables (Cahn &
Hilliard 1958). The governing conservation transport equations are then simply the
classical mass conservation equations for each of the two species, together with the
momentum conservation equation for the mixture as a whole; explicitly,
@» 2
@» 1
+ r (» 2 v2 ) = 0;
+ r (» 1 v1 ) = 0;
@t
@t
@(» v)
» i Fie ;
+ r (» vv) = ¡ rp + r ¿v + r ¿c +
@t
(2.3)
(2.4)
where » i , vi and Fie (i = 1; 2) respectively denote the species-speci­ c mass density,
velocity, and external force ­ elds.y As usual, p denotes the pressure ­ eld, ¿v the viscous deviatoric stress tensor, » the mass density of the mixture (» = » 1 + » 2 ), and v
the mass-average velocity of the mixture (cf. (2.7)). The continuous microscale vector
­ eld r ¿c ² Fp , say, represents the precursor of the macroscale physicochemical
capillary or interfacial forces, arising as a consequence of the steep species concentration gradients @» i =@n existing within the di¬use interfacial region in a direction
n normal to the `interface’. The implicit macroscale interfacial tension and related
non-equilibrium capillary e¬ects traditionally associated with the presence of a singular surface (which includes interfacial viscositiesz) are thus explicitly accounted
for in our model by the appearance of this force ­ eld. As such, the physicochemical
volumetric force density ­ eld Fp constitutes a central feature of our di¬use interface
model, and the ultimate success of this model depends upon our ability to provide
a rational constitutive equation for Fp |an issue deferred until the next section (cf.
(2.18)).
We assume at the outset that the mixture as a whole is incompressible. This
assumption allows the microscale pressure ­ eld p to be treated as an independent
dynamical variable rather than as a dependent variable, functionally dependent upon
temperature, composition, etc. Further, upon supposing that the law of additive
volumes applies to the mixture, the relation between the density » of the mixture
and the densities » 0i of the pure components is given simply by the linear relation
(Antanovskii 1995),
» = » 01 c + » 02 (1 ¡
c);
(2.5)
y The presence of a body-force ­ eld F i e , e.g. a gravity ­ eld, ultimately precludes the scenario of
spatially periodic microscale ­ elds envisioned in later sections in cases where the immiscible phases
possess di¬erent densities. However, to maintain generality of the exposition, we eschew that constraint
until x 3 m.
z In the singular interface model, interfacial viscosities are associated with the adsorption of surfactant
species on the `interfacial surface’. In our di¬use model of the interface, the presence of a surfactant
species in the two-phase mixture would constitute a ternary system. Though most of our arguments in
subsequent sections remain valid in this ternary case, we do not formally delve into the complexities of
such situations in this preliminary report, wherein we con­ ne ourselves exclusively to binary systems.
Proc. R. Soc. Lond. A (2000)
744
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
where c denotes the volume fraction of component 1. The latter concentration is
henceforth taken as the relevant composition variable when characterizing our binary
system. The assumption of incompressibility shows that at the microscale only one
composition variable is relevant, namely, the volume fraction c (since » 1 = c» 01 , and
» 2 = (1 ¡ c)» 02 ). This observation also makes the number of scalar equations de­ ned
by (2.3){(2.4), namely ­ ve, consistent with the corresponding number of unknown
scalar variables embodied in c; p and v.
Equations (2.3) may be rewritten as
@» 1
+ r (» 1 v) = ¡ r J1 ;
@t
@» 2
+ r (» 2 v) = ¡ r J2 ;
@t
(2.6)
wherein the mass-average velocity v is de­ ned by the relation,
» v = » 1 v1 + » 2 v2 ;
(2.7)
and the respective di¬usion ®uxes J1 and J2 as
J1 = » 1 (v1 ¡
v);
J2 = » 2 (v2 ¡
v):
(2.8)
(In view of the de­ nition (2.7) of the mixture velocity v, the two di¬usion ®uxes are
related by the identity J1 + J2 = 0, whence we can choose to denote J1 = ¡ J2 ² J ,
say, and so by adding equations (2.3) obtain the conventional equation of continuity
for the mixture:
@»
(2.9)
+ r (» v) = 0;
@t
together with a di¬usion equation for one of the species, say 1:
@c
+ r (cv) = ¡ r J :
@t
(2.10)
While this index-free notation representation will be the form employed in the
next section when dealing with constitutive equations, considerations of invariance
a¬orded by the interchangeability of the arbitrary labels identifying species `1’ and `2’
suggest that for the present we opt for the individual species-speci­ c forms, namely
(2.6).)
Equations (2.6) and (2.4) together with appropriate boundary conditions dictated
by geometrical and kinematical considerations constitute a complete set of conservation equations serving to determine the di¬use microscale ­ elds c; p; v. These equations are identical to the conventional multicomponent equations used to describe the
®ow of miscible, single-phase, binary mixtures (solutions) except for the appearance
of the additional term Fp , representing the microscale precursor of the macroscopic
capillary forces. In the next section it will be shown that the sole di¬erence between
the respective miscibility and immisciblilty cases arises in toto from the forms of the
constitutive equations chosen for Fp and J (= J1 ² ¡ J2 ). This seemingly subtle
feature is the main thrust of our hypothesis. It shows how even immiscible two-phase
systems may be alternatively regarded as two-component, single-phase systems when
the system description is appropriately supplemented by constitutive equations for
the di¬usive ®ux J and capillary force ­ eld Fp (Fixman 1962). The fact that such
constitutive equations can be provided and embedded in physical context renders our
work novel beyond its classi­ cation as a mere academic exercise. Section 2 c, which
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
745
follows, thereby constitutes in a fundamental sense a prescription for performing
explicit two-phase ®ow calculations using the di¬use interface model. Furthermore,
it shows in an explicit manner how the traditional equations governing the singular
interface model can be recovered from the di¬use model through the natural apparatus available for identifying and e¬ecting such microscale/macroscale equivalences,
namely, singular perturbation analysis (Mavrovouniotis & Brenner 1993). Though
some related aspects of the subsequent analysis have been dealt with in prior studies
by others, we take pains to clarify the subtleties involved in such an exercise, enabling
us thereby to transcend the purely numerical or mathematical details of the scheme.
(c) Microscale constitutive equations
This section uses the respective forms (2.9) and (2.10) of the di¬usion and continuity equations in place of the precursor pair (2.6). Moreover, constitutive expressions
for J and Fp are discussed.
(i) Constitutive equation for J
The assumption of linear response theory together with the assumed absence of
coupling between the ®uxes requires that
J =¡
¤
r· ;
kT
(2.11)
where · denotes the di¬erence in chemical potential between the two species, and kT
is the Boltzmann factor. The mobility coe¯ cient ¤ is, in general, a second-rank tensor and may also depend upon c. However, consistent with the nature of the present
work, stressing only on the simplest physical elements, we eschew both generalizations. Thus, we assume ¤ to be a concentration-independent scalar constant. For the
chemical potential · we have by de­ nition that
· =
¯ F
;
¯ c
(2.12)
with F is the free energy density of the mixture and ¯ =¯ c the functional derivative
with respect to c. The constitutive form of F is assumed to be that proposed by
Cahn & Hilliard (1958), namely
F =
d3 r [f (c) + 12 K(rc)2 ];
(2.13)
where f (c) denotes the excess volumetric free energy density of the mixturey and
K is a dimensionless phenomenological coe¯ cient, related to the equilibrium direct
correlation function between the two species (Fleming et al . 1976; Davis & Scriven
1982). Though K may be a function of concentration, it is here assumed to be a
constant. In this work we use a modi­ ed form of the above equation to explicitly
y
The actual free energy density of the mixture is of the form:
f 0 = f (c) + fid +
1
2
K (rc)2 ;
where fid = cf1 (pure) + (1 ¡ c)f2 (pure) denotes the free energy density of an ideal mixture, and f 0
denotes the microscale volumetric free energy density. However, in most applications we are concerned
only with r· , rather than · itself. As such, fid proves irrelevant.
Proc. R. Soc. Lond. A (2000)
746
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
incorporate the interfacial `thickness’ ¹ , the latter representing the characteristic
linear dimension of the interfacial transition zone:
F =
1
¹
d3 r [f (c) + 12 K¹ 2 (rc)2 ]:
(2.14)
Thus, J takes the form
J =¡
@f
¤ 1
¡
r
@c
kT ¹
K¹ 2 r2 c :
(2.15)
Substitution into the di¬usion equation (2.10) yields
@f
1 ¤
@c
¡
r2
+ r (vc) =
@c
¹ kT
@t
K¹ 2 r2 c
(2.16)
for the microscale convective-di¬usion equation (see also Fixman (1962) and Antanovskii (1995)). The constitutive form of f (c) is taken up subsequently.
(ii) Constitutive equation for Fp
An appropriate constitutive relationship for Fp has been suggested by a number
of researchers, frequently starting from the virtual work principle (Fixman 1962;
Starovitov 1994; Antanovskii 1995), as well as other arguments.y We adopt the following constitutive expression for ¿c :
¿c = ¹ K[(rc)(rc) ¡
jrcj2 I];
(2.17)
whence from the de­ nition r ¿c = Fp of the physicochemical force density we obtain
Fp = ¹ Kr [(rc)(rc) ¡
jrcj2 I]:
(2.18)
(iii) Constitutive equation for ¿v
In view of the fact that the interfacial transition region is strongly inhomogeneous
in a direction normal to the interface, the viscous stress tensor is expected to be
transversely isotropic (Eliassen 1963; Goodrich 1981). However, consistent with our
desire to focus attention primarily on the physicochemical (rather than rheological)
y Though conventionally derived from a variational formulation, we can heuristically justify the
expression (2.17) for ¿ c in the following manner: De­ ne a pressure
d ef
pc = c
¯ F
¡F
¯ c
in analogy to the thermodynamic de­ nition
p=c
@f
¡ f:
@c
This enables a separation of rp c into the sum rp + r ¿ c , with
rp = r c
wherein ¿ c is given as above.
Proc. R. Soc. Lond. A (2000)
@f
¡f ;
@c
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
747
aspects of the dynamics, we adopt the usual simple isotropic Newtonian form for ¿v ,
namely
1
(I
3
¿v = 2· [D ¡
: D)I];
(2.19)
with
D = 12 [(rv) + (rv)y ]
(2.20)
the deformation gradient and · ² · (c) the local viscosity of the mixture (solution),
which by de­ nition satis­ es the relationy
· ¡ !
·
·
0
1
0
2
as c ¡ ! 1;
as c ¡ ! 0:
(2.21)
Upon using the above identi­ cations, the momentum equation (2.4) adopts the
form
@(» v)
+ r (» vv) = ¡ rp + r ¿v + ¹ Kr [(rc)(rc) ¡
@t
jrcj2 I] + » g:
(2.22)
(iv) Form of f (c)
Conventional theories of dynamical critical phenomena (wherein our model is
equivalent to that of Model H of Hohenberg & Halperin (1977)) adopt a double-well
potential for the free energy density f (Á), where Á represents an order parameter. In
a binary liquid system the order parameter is usually taken to denote the di¬erence
in composition between the two species. Thus, we choose
f (Á) ¹
Á 2 )2 ;
(1 ¡
(2.23)
where
jÁj = j1 ¡
2cj:
Consequently, the following form is adopted here for f (c):
f (c) = AK(1 ¡
c)2 c2 ;
(2.24)
where A is a dimensionless constant and K is the same phenomenological constant
introduced earlier. Note that f (c) ! 0 as c ! 0 and 1; this dual property will later
prove crucial in identifying the compositions of coexisting phases.
(d ) Immiscibility
This section is devoted to clarifying the notion of immiscibility in the present context of species rather than phases, especially the manner in which immiscibility is
supported by the above equations. To achieve this objective we start with a brief
y Though we will have no need here for an explicit functional relationship between · and c, if necessary
one could, in order to satisfy (2.21), simplistically assume · to be linearly related to the respective
viscosities · 01 and · 02 of the pure species by the expression
· =·
Proc. R. Soc. Lond. A (2000)
0
1c
+·
0
2 (1
¡ c):
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
748
(b)
0
diffusivity, D(c)
free energy density, f (c)
(a)
0
0.25
0.5
0.75
concentration, c
0
0
0.25
0.75
0.5
concentration, c
1
1
Figure 2. Double-well potential form for (a) free-energy density f (c); and (b) di® usivity D(c).
discussion pertinent to the existence of phase-separated solutions of the modi­ ed
di¬usion equation (2.16). Eventually, we perform a singular perturbation analysis
of the governing equations, thereby explicitly demonstrating that our di¬use interface model does indeed furnish an appropriate pair of `bulk-phase’ equations (outer
equations) as well as yielding `interfacial’ transport equations (matching conditions)
connecting these two bulk ­ elds at the `interface’. And it is these bulk ­ eld regions
that constitute the `immiscible phases’ in the macroscale view.
In the interests of simplicity we temporarily ignore hydrodynamical considerations
as well as e¬ects arising from the presence of K in the expression (2.14) for F (c).
Figure 2 displays the general con­ gurations of the concentration-dependent forms of
the excess free energy density function f (c) appearing in (2.24) as well as the di¬usion coe¯ cient D(c), the latter being the term appearing in the resulting `di¬usion’
equation,
@c
= Dr2 c;
@t
(2.25)
and de­ ned as
d ef
D =
¤ @2 f
:
kT @c2
(2.26)
Since ¤ > 0 as a consequence of the Second Law of Thermodynamics, we observe from
(2.25) and (2.24) that there exists a range of concentrations over which D(c) < 0
(representing a thermodynamically unstable region, wherein an initially homogeneous mixture (solution) of composition 0 < c < 1 spontaneously separates into two
distinct `phases’ of respective compositions c = 0 and 1). In contrast, for a completely
miscible mixture, D(c) > 0 8 c, whence an initially homogeneous mixture remains
homogeneous for all time! This observation forms the basis of our two-component
mixture model as supporting phase-separated stable solutions of the modi­ ed di¬usion equation.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
749
This preceding discussion was meant to provide a brief overview for readers unfamiliar with the work of Cahn & Hilliard (1958). Excellent reviews and details can
be found in Langer (1974) and Gunton et al . (1983). However, our work is not concerned with the dynamics of phase separation, but rather assumes at the outset that
a phase-separated mixture already exists at the microscale, one whose dynamics can
be described via the above system of equations.
Some comments are in order regarding the form of F set forth in (2.14).
(i) The factor 1=¹ appearing as the multiplier of F is obtained by normalizing
the volume over which F is non-zero; ¹ is expected to be small for immiscible
systems.
(ii) While the above equations appear to be super­ cially similar to those possibly
encountered for miscible systems, the contrasting phase behaviour of the two
systems arises from the di¬erent forms of the excess free energy function f
present in F above. While the form of f illustrated in equation (2.23) supports
a phase-separated equilibrium scenario, the corresponding functional form of
f for a miscible system would describe an equilibrium scenario with partial
miscibility.
(iii) The factor K and the second gradient in F account for interfacial tension e¬ects
(Cahn & Hilliard 1958).
The next section is devoted to e¬ecting a singular perturbation analysis of the
preceding system of microscale equations so as to rationally derive obtain the conventional macroscale two-phase equations (2.1) and (2.2), together with appropriate
matching conditions. This section also serves to clarify the precise manner in which
analytical calculations of practical ®uid-mechanical problems need to be e¬ected
when using the di¬use interface model. Such a procedure necessarily involves the use
of singular perturbation techniques to determine the governing equations and their
solutions at every order of the perturbation framework.
(e) Singular perturbation analysis
Below, we non-dimensionalize the variables appearing in (2.16) and (2.22). To
avoid a proliferation of symbols, the same unadorned symbols will be used to denote
both dimensional quantities and their non-dimensional counterparts, imagining that
the preceding dimensional symbols had previously been augmented by an asterisk:
v=
» =
¤
»
»
;
o
v¤
;
U
p=
r=
p¤ L
;
· oU
r¤
;
L
¿v =
r = Lr¤ ;
¿v¤ L
;
· oU
¤ =
¤
¤
½ D
:
L2
Here, U , L, » o and · o respectively represent a characteristic velocity, length, density
and viscosity, whereas ½ D represents a characteristic di¬usion time-scale. In terms of
Proc. R. Soc. Lond. A (2000)
750
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
these dimensionless variables, equations (2.9), (2.16) and (2.22) respectively become
@»
+ r (» v) = 0;
@t
Re
(2.27)
@(» v)
+ r (» vv) = ¡ rp + r ¿v + Ca¡1 ¯ r [(rc)(rc) ¡
@t
jrcj2 I] + » g
(2.28)
and
¯ ½
D
@f
U @c
¡
+ r (vc) = ¤ r2
@c
L @t
¯ 2 r2 c ;
(2.29)
where the following non-dimensional parameters have been introduced:
Re =
» oU L
;
·
¯ =
¹
;
L
Ca =
· U
;
K
¤ =
¤ K
;
LkT
and where f is now given by the expression
f (c) = A(1 ¡
c)2 c2 :
(2.30)
Re and Ca are, of course, the Reynolds and capillary numbers.
For the immiscible systems contemplated, ½ D is expected to be very small compared
with the convective time-scale L=U . In this section, circumstances will be established
as to how the above system of equations reduce to the conventional system of twophase, singular interface equations. The motivation behind such an exercise is to
explain the spirit of the di¬use interface model. A description of the actual dynamics
of phase separation (Pego 1989) is beyond the scope of this work; in lieu of this we
make an a priori assumption as to the existence of an interface characterized by a
transition region with steep variations in physical properties.
The natural mathematical apparatus for extracting the bulk equations and matching interfacial conditions from the above system of equations is singular perturbation
analysis, with the interfacial region representing a thin boundary layer sandwiched
between two outer bulk regions (phases). To the best of our knowledge the only
such interfacial analysis previously e¬ected via this approach is embodied in the
work of Starovitov (1994) (see also Caginalp (1986) for an analysis with similar
objectives, albeit in a di¬erent context). However, that analysis was restricted to a
one-dimensional interface; moreover, a number of aspects of that work appear controversial. Here, we treat a moving material interface in three-dimensional space,
characterized by a special metrical geometry specially chosen to avoid notational
complexities. In this context the analysis of systems containing mobile interfaces is
facilitated by parametrizing three-dimensional space via a semiorthogonal surface­ xed coordinate system (cf. ­ gure 3), namely (q1 ; q2 ; q3 ), with q3 lying normal to the
`parent surface’ (Mavrovouniotis & Brenner 1993; Edwards et al. 1991). Since our
aim is merely to illustrate the spirit of the di¬use interface model without performing
an exhaustive and comprehensive investigation to rigorously justify the applicability of our results under very general circumstances, we consider a special kind of
interface|one in which the interfacial transition region can be characterized by an
orthogonal coordinate system possessing the metrical coe¯ cient h3 = 1 (the `parallel
surface’ representation of Eliassen (1963)). Additionally, the interfacial curvature will
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
x(q1, q2, n)
parent
surface
P
n
surface-fixed
origin
751
q2
n=0
q1
q3 = q 03
O
space-fixed origin
Figure 3. Parent surface and representation of the surface-¯ xed coordinate system; (q1 ; q2 ; n)
represent the coordinates of a point in the outer region as measured from axes ¯xed on the
parent surface.
be assumed to be much greater than the interfacial thickness (and of O(1) relative
to the small perturbation parameter ¯ ).
As is conventional with di¬use interface models we identify the macroscopic interface with a so-called `parent surface’, one whose precise location in the transition
region to within an error of O(¯ ) is not critical to the subsequent ¯ ½ 1 asymptotic analysis. This parent surface is parametrized by the normal coordinate q3 = q30 ,
whence the variable
d ef
q3
h3 dq3 = q3 ¡
n =
q30
q30
(2.31)
denotes an algebraically signed distance measured normal to interface. By de­ nition,
n = 0 on the parent surface. Given the presence of the small perturbation parameter ¯
in the equations governing the two-component transport process, one can anticipate
the existence of disjoint inner and outer regions. Length-scales in the two outer
regions lying on either side of the `interface’ are taken to be macroscopic, with the
magnitudes of physical quantities and their gradients in these regions assumed to be
everywhere of O(1).
In the outer regions we assume, subject to a posteriori veri­ cation, that each physical variable (here denoted generically by Ã) admits a regular perturbation expansion
of the form
Ã(q1 ; q2 ; n; t; ¯ ) = Ã 0 (q1 ; q2 ; n; t) + ¯ Ã1 (q1 ; q2 ; n; t) +
:
(2.32)
In contrast, in the inner region, whose thickness is of O(¯ ), these physical quantities
may display steep gradients. Accordingly, distances in the inner region are normalized
Proc. R. Soc. Lond. A (2000)
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
752
with ¯ , thereby magnifying them. Explicitly, the independent inner variable,
d ef
q30 )
(q3 ¡
;
(2.33)
¯
represents an algebraically signed scaled distance in the inner region, normal to the
~ and r
~ 2 need to be explicitly rewritten
`interface’. Concomitant expressions for r
in terms of (q1 ; q2 ; n
~ ). Furthermore, subject to a posteriori veri­ cation, the generic
physical variables are each assumed to possess an inner expansion of the form
n
~ =
~ 1 ; q2 ; n
Ã(q
~ ; t; ¯ ) = Ã~0 (q1 ; q2 ; n
~ ; t) + ¯ Ã~1 (q1 ; q2 ; n
~ ; t) + : : : :
(2.34)
Matching conditions imposed on the respective outer and the inner ­ elds, Ãi and
Ã~i , have been established by Caginalp & Fife (1988), albeit in a di¬erent context, as
lim
j~
nj¡!
1
Ã~0 (q1 ; q2 ; n
~ ; t) = Ã0 (q1 ; q2 ; n = 0; t)
(2.35)
and
lim
j~
nj¡!
1
@Ã0 (q1 ; q2 ; n; t)
Ã~1 (q1 ; q2 ; n
~ ; t) = Ã1 (q1 ; q2 ; n = 0; t) + n
~
@n
:
(2.36)
n= 0
Higher-order matching conditions may also be easily derived; however, the need for
such additional conditions does not arise in the present work as a consequence of the
fact that our theory is strictly asymptotic (Mavrovouniotis & Brenner 1993) rather
than being serially sequential. As such, only the leading- and ­ rst-order terms of the
respective expansions are involved in the subsequent analysis.
(i) Outer equations
This section addresses the leading-order outer ­ elds. A¯ xes § will be used to
denote the respective regions lying on either side of the interface. Upon substituting
the outer expansion into (2.29), the O(1) di¬usion equation is found to be of the
form
r2
§
@f
@c
= 0;
(2.37)
0
where the subscript zero denotes the zeroth-order ­ eld. Since a soliton solution of
the above equation is anticipated, we assume subject to a posteriori veri­ cation that
the appropriate solution of the above equation is
@f
@c
§
= 0:
0
As noted earlier, f (c) ! 0 and @f =@c ! 0 as c ! (0; 1). Consequently, we ­ nd
that c0 = (0; 1)8 n lying on either side of the interface. This is consistent with our
physical picture of the two outer regions as constituting immiscible bulk phases. For
de­ niteness, we arbitrarily choose c+0 = 1 and c¡
0 = 0.
At O(1) the equation of continuity becomes
@» §
0
+ r (»
@t
Proc. R. Soc. Lond. A (2000)
§ §
0 v0 )
= 0:
(2.38)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
753
In view of the preceding solution for c0 on either side of the interface, it follows that
r v0§ = 0:
(2.39)
The O(1) momentum equation here takes the form
§
D(» 0 v0 )
Dt
Re
§
= ¡ rp§
0 + r ¿v0 + »
§
i0 g:
(2.40)
Equations (2.39) and (2.40) represent the outer equations in the respective regions
lying on either side of the interface. These equations accord with the usual bulk-phase
equations employed in conventional singular surface, two-phase models of interfacial
phenomena. Considerations pertaining to the inner equations will be shown to lead
to appropriate matching conditions imposed on these outer ­ elds at their common
interface.
(ii) Inner equations
In the inner region the normal-distance coordinate n needs to be rescaled via the
transformation (2.33). Relevant rescalings of the operators thereby obtained, along
with the leading- and the ­ rst-order terms that arise, are indicated in Appendix B.
Here, we note that the dynamical evolution of the interface leads to the following
rescaling of the time-derivative (Caginalp & Fife 1988):
@
@t
=
n
@
@t
q_30
¯
¡
n
~
@
@~
n
(q1 ; q2 = const:);
(2.41)
t
where q_30 denotes the normal velocity of the parent surface, n
~ = 0. Furthermore,
~ 2 (q1 ; q2 ; n
the metrical coe¯ cients ~
h1 (q1 ; q2 ; n
~ ) and h
~ ) admit inner expansions of the
respective forms,
~ 1 (q1 ; q2 ; n
h
~ ) = h10 (q1 ; q2 ; q30 ; t) ¡
~ 2 (q1 ; q2 ; n
h
~ ) = h20 (q1 ; q2 ; q 0 ; t) ¡
3
¯ h10 (q1 ; q2 ; q30 ; t)µ1 (q1 ; q2 ; n
~ ; t)~
n+
;
(2.42)
h20 (q1 ; q2 ; q30 ; t)µ2 (q1 ; q2 ; n
~ ; t)~
n
;
(2.43)
¯
+
where the de­ nitions (B 4) and (B 5) has been used. Higher-order contributions prove
unnecessary in the subsequent analysis.
At the leading order the continuity equation becomes
¡ q_30
@(~
» 0 v~30 )
@ » ~0
= 0:
+
@~
n
@n
~
(2.44)
Integration of the latter at constant (q 1 ; q 2 ) gives
» ~0 (~
v30 ¡
q_30 ) = g(q1 ; q2 );
(2.45)
where g(q1 ; q2 ) is an integration function. However, since v~30 = q_30 at n
~ = 0 owing to
the material nature of the interface, it follows that
v~30 = q_30
8n
~:
(2.46)
Since q_30 is independent of n
~ , this implies that
@~
v30
= 0:
@n
~
Proc. R. Soc. Lond. A (2000)
(2.47)
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
754
Moreover, the matching conditions (2.35) require that
§
v30
(n = 0) = q_30 ;
(2.48)
[[v30 ]] = 0;
(2.49)
whence
where we have used the jump operator [[: : : ]] of Edwards et al. (1991) to denote the
discontinuity, if any, in the value of the argument across the parent surface.
At the leading order, the inner form of the di¬usion equation yields
@f
¡
@c
@ 2 c~0
= 0:
@n
~2
(2.50)
Multiplying the above equation by @~
c0 =@ n
~ and integrating with respect to n
~ gives
f¡
c0
1 @~
~
2 @n
2
= 0;
(2.51)
where the facts that both @~
c0 =@ n
~ and f ! 0 as n
~ ! §1 have been used. Into the
above expression one can substitute the explicit form of f given in (2.30), subsequently integrating the resulting equation to obtain the leading-order inner concentration pro­ le as
p
exp(2~
n A)
p :
(2.52)
c~0 =
1 + exp(2~
n A)
(Note that c~0 ! 0 as n
~ ! ¡ 1 and c~0 ! 1 as n
~ ¡ ! 1:)
The leading-order inner components appearing in the momentum equation are
¡
@ v~10
@
@ p~¡1
·
i3 + i1
@n
~
@n
~
@n
~
+ i2
@ v~20
@
·
@n
~
@n
~
+ i3
@
@n
~
4
·
3
@~
v30
@~
n
= 0;
(2.53)
where (i1 ; i2 ; i3 ) respectively denote unit vectors along the coordinate axes (q 1; q 2; q 3 ).
Additionally, p~¡1 denotes the O(¯ ¡1 ) inner pressure ­ eld. Using (2.47) together with
the fact that v~10 and v~20 are ­ nite as n
~ ! §1 gives
v~10 = ¸ (q1 ; q2 );
v~20 = ² (q1 ; q2 );
p~¡1 = 0;
(2.54)
(2.55)
(2.56)
in which ¸ and ² denote integration functions. The ­ rst pair of the preceding equations require that
@~
v20
@~
v10
= 0;
=
@n
~
@n
~
(2.57)
[[v10 ]] = [[v20 ]] = 0:
(2.58)
whence
In combination with (2.49) we thereby obtain
[[v0 ]] = 0;
Proc. R. Soc. Lond. A (2000)
(2.59)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
755
which is equivalent to (B 13) (see Appendix B). As such, our considerations of the
leading-order terms of the inner equations lead to the expected jump condition (2.59)
imposed on the velocity ­ eld for the case of material interfaces (Edwards et al . 1991).
At the next order the inner form of the continuity equation requires that
v~31
@ » ~0
+ » ~0 À = 0;
@n
~
(2.60)
where (2.46) and (2.47) have been used, and wherein À is as de­ ned in (B 11).
Consider the n
~ ! §1 limit of the above equation. Since, as a consequence of (2.5)
and (2.52), @ » ~0 =@ n
~ ! 0 exponentially as n
~ ! §1, it follows that
À !0
as n
~ ! §1:
(2.61)
Using (B 10) and (B 12) the momentum equation at ­ rst order reads
@
·
@n
~
@
·
@n
~
¡
@ v~31
@
@ p~0
2·
+
@n
~
@~
n
@n
~
¡
@~
v30
@~
v11
¡
+ h10
@q1
@n
~
@~
v30
@~
v21
¡
+ h20
@q2
@n
~
v~10 µ1
= 0;
(2.62)
v~20 µ2
= 0;
(2.63)
= 0:
(2.64)
2
@~
c0
2 @(· À )
+ Ca¡1 (µ1 + µ2 )
@n
~
n
3 @~
Integrate (2.62) and (2.63), and use the matching condition
@Ã0 (q1 ; q2 ; n = 0; t)
@ Ã~1 (q1 ; q2 ; n
~ ; t)
=
@n
@n
~
lim
n
~ ¡!
1
(2.65)
together with (2.54) and (2.55) to obtain the jump conditions
·
h10
@v
@v 30
+ 10 ¡
@n
@q1
v10 µ1
=0
(2.66)
h20
@v20
@v 30
¡
+
@n
@q2
v20 µ2
= 0:
(2.67)
and
·
Integration of (2.64) yields
¡ p~0 + 2·
@~
v31
@n
~
¡
2
·
3
À + Ca¡1 (µ1 + µ2 )
d~
n
@~
c0
@~
n
2
= ± (q1 ; q2 );
where the right-hand term represents an integration function.
Consider the limit as n
~ ! §1 together with the fact that À
This yields
¡ p0 + 2·
Proc. R. Soc. Lond. A (2000)
@v 30
@n
= ¡ 2¼ H;
(2.68)
! 0 as n
~ ! §1.
(2.69)
756
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
with the mean curvature H de­ ned as in (B 16), and the macroscopic interfacial
tension ¼ de­ ned (non-dimensionally) in terms of the continuous microscale ­ eld c~0
as
¼ = Ca¡1
+ 1
d~
n
¡1
@~
c0
@~
n
2
:
(2.70)
The latter expression represents a classical result, one which has been derived innumerable times for equilibrium situations (Cahn & Hilliard 1958; Fleming et al. 1976),
albeit by perhaps less formal procedures. However, our analysis appears to be the ­ rst
to systematically treat the complete dynamical situation case, so as to also derive
the interfacial matching conditions imposed on the stress and velocity ­ elds. Though
we have ignored possible surface rheological e¬ects arising from the presence of an
adsorbed surfactant, it would be straightforward to incorporate these into the above
model by using appropriate constitutive equations and scalings, and by extending
our analysis to include ternary systems (with the surfactant representing the third
species).y Moreover, by retaining higher-order terms in the perturbation scheme it is
possible to proceed along similar lines to derive the well-known e¬ects of curvature
on interfacial tension (Bu¬ & Saltsburg 1957; Ono & Kondo 1960). However, in the
interests of focusing only on the fundamentals of di¬use interfaces we do not embark
on such extensions in the present contribution.
In summary, the di¬use interface model provides a complete physicochemical identi­ cation of the bulk equations and interfacial boundary conditions employed in conventional singular interface models, albeit one that eschews both interfacial rheology
and large curvature e¬ects. The present section was devoted to rationalizing the
admissibility of incorporating two-phase immiscible systems into the framework of
a two-component mixture model of interpenetrating continua. The microscale constitutive equations introduced were shown to reduce to the conventional singular
interface model through the use of physical scaling arguments accompanied by an
asymptotic singular perturbation scheme. In accomplishing this demonstration, we
have also simultaneously illustrated the manner in which analytical solutions of problems can be e¬ected using the di¬use interface model. Such a procedure entails use of
singular perturbation techniques concurrently with prescribed constitutive equations
to obtain solutions accurate to the required degree of approximation.
Henceforth, the di¬use interface model will be used exclusively in the subsequent
analysis to treat multiphase ®ows at the microscale level, appropriate to pore-level
interstitial transport processes, in lieu of the traditional singular interface pore-level
model. However, until the stage where x 3 j is reached, the subsequent analysis will
be seen to be invariant to the explicit choice of constitutive equations.
Sections 3 c{m serve to apply the preceding single-phase (di® use), two-component, interstitial microscale description of macroscopically multiphase systems to
problems involving their ®ow through porous media. Such two-phase ®ow processes
have conventionally been regarded as governed at the (singular surface) microscale
level by equations (2.1) and (2.2) together with appropriate boundary conditions
imposed at the bed-particle surfaces. Our alternative di® use interface model, however, requires that we instead concern ourselves with the continuous equations (2.4)
y For instance, the viscosity would no longer be of O (1) within the interfacial zone, but might rather
display very large values, eventually leading to non-zero interfacial viscosities (Eliassen 1963; Goodrich
1981; Edwards et al . 1991; Mavrovouniotis et al . 1993).
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
757
and (2.6) together with appropriate bed-particle boundary conditions imposed on
the ­ elds appearing in these equations. These equations form the basis of our analysis of biphasic ®ows in porous media. This multicomponent description at the di¬use
microscale leads naturally to a multicomponent description at the di¬use macroscale,
and hence ultimately to an assignation of distinct macroscale pressures to each of
the two `species’ (i.e. `bulk phases’). The relationship between these two bulk-phase
pressures is shown, in turn, in x 3 i to be the formal equivalent of the usual capillary
pressure relationship connecting these `phase-speci­ c’ pressures. The entire analysis
which follows is performed in a manner such as to admit completely general functional forms for the constitutive equations governing J and Fp . However, in x 3 j
we use the constitutive equations outlined above to provide a number of new and
fundamental insights into the physical basis of capillary pressure as well as into the
dynamical equations currently used in the modelling of two-phase ®ows in general.
3. Part II. Applications to biphasic ° ows in porous media
(a) Geometry of unbounded spatially periodic porous media
(i) Microscale
A spatially periodic ­ xed bed-particle con­ guration and a corresponding (instantaneously) spatially and time-periodic mobile interstitial ®uid velocity ­ eld is proposed as a tractable geometrical model of multiphase ®ow through a porous medium.
Owing to its deterministic con­ guration, this speci­ c geometrical microcontinuum
model permits a rigorous mathematical analysis of the pertinent physical phenomena
using standard mathematical tools, including Fourier series, Bloch functions (Brillouin 1946), and homogenization techniques (Bensoussan et al. 1978). Such models
have already been used for computational purposes with regard to emulsions, suspensions, etc. (cf. Nitsche & Brenner (1989) and the references cited therein), albeit
at the singular surface microscale level. We are explicitly concerned here with macroscopically homogeneous multiphase ®ows, where time-averaged macroscale quantities
like saturation, velocity, pressure gradient, etc., are spatially and temporally uniform
throughout the porous medium. However, we will also brie®y indicate the more general forms of equations that would result for circumstances in which a slow spatial
macroscale variation is superposed on the otherwise spatially periodic con­ guration
(for which our rigorous macroscopic de­ nitions hold locally ). While the model of a
spatially periodic medium might appear unnecessarily restrictive in view of the fact
that the entire porous medium is assumed generated by replication of the unit cell
contents (in contrast with the disordered media usually encountered in practice),
it is to be pointed out here that such media do not, however, exclude a complex
disordered unit cell as the template (Nitsche 1989).
Reviews of the properties of spatially periodic media are available elsewhere (Brenner 1972, 1980; Nitsche 1989; Adler 1992; Brenner & Edwards 1993) and will not be
repeated here. Rather, we restrict ourselves in what follows to a brief summary of
the notation to be used. Moreover, attention is focused exclusively on unit cells characterized by planar faces. Extension to curvilinear cells is straightforward (Nitsche
& Brenner 1989).
Spatially periodic systems are characterized by an associated lattice system embodying the translational symmetries of the medium. The complete lattice is generProc. R. Soc. Lond. A (2000)
758
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
t
o
a
lattice points
t
a’
f
l1
particles
l2
r
Rn
origin, O
Figure 4. Spatially periodic medium (two-dimensional projection) comprised of particles of arbitrary shape. O refers to the lattice origin, R n the centroid of the cell, and r a point within a
unit cell as measured from the centroid; a, a0 denote representative lattice points, whereas ½ f ; ½ o
respectively denote the domains of the interstitial ° uid and the unit cell.
ated from a single lattice point (the origin) by discrete displacements of the form
R n = n 1 l1 + n 2 l2 + n 3 l3 ;
(3.1)
where fn1 ; n2 ; n3 g ² n is a triplet of integers (n1 ; n2 ; n3 = 0; §1; §2; §3; : : : ); and
(l1 ; l2 ; l3 ) a triad of non-coplanar vectors, termed basic lattice vectors. As in the
references cited in the preceding paragraph, the position vector R of a point denotes
the displacement of the point from a space-­ xed origin, whereas r = R ¡ Rn denotes
the intracellular position vector measured from the lattice point Rn , the latter being
taken as a local intracell origin. Furthermore, as regards a unit cell, @½ o denotes the
areal domain encompassing the totality of all six faces bounding the cell, whereas ½ f
and ½ o denote the respective volumetric domains of the interstitial ®uid and super­ cial cell (­ gure 4). Within a unit cell the solid particles comprising the porous medium
are stationary. Fluid velocities are to be measured with respect to these ­ xed particles. The wetted surfaces of the ­ xed bed particles within a cell are designated
collectively as sp .
(ii) Macroscale
In order that a continuum macroscale description of the microscale phenomena be
valid, there needs to exist a wide disparity in the magnitudes of the characteristic
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
759
lengths (and times) separating the micro- and macroscales. Most theories of multiphase ®ow in porous media (Marle 1977; Whitaker 1986; Gray & Hassanizadeh 1990;
Nitsche & Brenner 1989) assume the existence of an intermediate length-scale lying
between the two|thereby serving to identify a macroscale `point’ with a non-zero
microscale skeletal-space volume. Similarly, a macroscale `di¬erential volume’ element is in­ nitesimal only in macroscale terms; in reality it needs to be su¯ ciently
large to enclose a `statistically representative’ sample of the local skeletal geometry (as well as its contents). Quantitative details regarding the errors involved in
homogenizing the microscale phenomena by this scheme are set forth in the work of
Nitsche & Brenner (1989), with special emphasis on the kinematical aspects of the
®ow phenomena.
Since interest centres in this paper only on the physical aspects of two-phase ®ows
in porous media, we concentrate on a rudimentary description of the micro{macro
transition, as expounded by Brenner (1972) (see also Brenner & Edwards (1993)).
In this description a `macroscopic point’ (whose macroscopic position vector will be
denoted by R) refers to a conveniently chosen reference lattice point located within
the microscale unit cell Rn (e.g. at its centroid):
R ² Rn :
(3.2)
While the lattice points Rn denote a discrete, in­ nitely denumerable set of points
at the microscale, the macroscopic description treats R as a continuum of vectors,
with the `in­ nitesimal’ displacement vector de­ ned as
d R = lk ;
(3.3)
wherein lk (k = 1; 2; 3) denotes a basic lattice vector. More explicitly, to remove any
confusion between the equivalence of the continuous vector ­ eld R and the discrete
set of vectors Rn , in place of (3.2) we heuristically write R = Rn + (klk), where l
denotes the characteristic linear dimension of a unit cell.
The macroscopic di¬erential directed areal surface element is identi­ ed with any
one of the six directed surface areas of the cell faces:
d s¸ = s k ;
(3.4)
with s¡k = ¡ sk (k = §1; §2; §3), and where the index ¸ connotes the direction
of the normal vector to the face. Henceforth, however, we will resort to the notation ds to represent the macroscopic areal element, unless a speci­ c need exists for
explicitly identifying the particular face in question. (Proof of the facts that dR and
ds respectively constitute acceptable de­ nitions for the in­ nitesimal displacement
vector and areal element at the macroscale (explicitly, that both can be chosen to
be of arbitrary magnitude and orientation in a macrosense) is discussed by Brenner
(1972) and Nitsche & Brenner (1989).)
The macroscopic in­ nitesimal volume element is taken to be the super­ cial volume
of the unit cell:
d3 R = ½ o ;
a fact consistent with (3.3) in conjunction with the relation ½ o = l1 £ l2 l3 : It follows
from (3.3) and (3.4), together with the identity
½
(valid for k = 1; 2 or 3), that
Proc. R. Soc. Lond. A (2000)
o
= lk s k
d3 R = ds dR:
760
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
(b) Microscale ¯elds
Due to the spatially periodic geometry of the porous medium, the instantaneous
microscale ­ elds » ; c and v can be expected to be spatially periodic with periods
fl1 ; l2 ; l3 g along the three lattice directions. In such a scenario we have, using (3.1),
the generic relation
Ã(r + Rn ; t) = Ã(r; t)
8 r 2 ½ f;
(3.5)
where à collectively denotes the ­ elds (c; » ; v) as well as » i . As a speci­ c application
of the above formula we have that:
kÃk = 0;
(3.6)
where kÃk denotes the (algebraically signed) jump in the value of à at equivalent
points lying on opposite cell faces.y The microscale spatial periodicity condition
embodied in (3.6) eventually enables us to de­ ne macroscopic quantities as averages
over a unit cell, a procedure carried out in the next section.
Under the assumption that ¿v and ¿c are functions at most only of (i) the velocity
v, (ii) the concentration c, and (iii) their respective gradients, we have
¿v (R + Rn ; t) = ¿v (R; t)
(3.7)
¿c (R + Rn ; t) = ¿c (R; t):
(3.8)
and
In view of the above facts it follows from (3.12) that the pressure gradient is also
spatially periodic:
rp(R + Rn ; t) = rp(R; t):
(3.9)
The time-dependence of the various microscale ­ elds requires special comment.
On the whole, microscale issues regarding the temporal evolution of two-phase ®ow
­ elds are rarely discussed in the literature, especially in circumstances where the
mean ­ elds are themselves time independent. The implicit assumption of quasistationarity thereby implied tends to be all pervasive. Ad hoc `volume-averaging’
is usually applied to some `random quasi-static con­ guration’, the latter taken to
denote some `representative’ time-averaged geometric con­ guration. However, at the
(singular surface) microscale level, two-phase ®ows of the respective ­ elds inherently
involve time evolution of the pertinent ­ elds due to the relative motion of the twophases with respect to one another, as well as with respect to the ­ xed particles comprising the porous medium (see also Adler et al . (1985), Nitsche & Brenner (1989),
and the comments therein). The corresponding macroscale ®ow may either be steady
(the common assumption under which the multiphase Darcy’s law is employed) or
unsteady. In the former case, the manner in which a microscopically unsteady motion
can result in a macroscopically steady motion leads in our analysis to the formal concept of a `micro{macro’ time-scale. To the best of our knowledge this feature has not
previously been explicitly identi­ ed in the context of biphasic ®ows, though various
y This unit cell face operator k : : : k should be carefully distinguished from the interfacial jump operator [[: : : ]] de­ ned following equation (2.2).
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
761
arguments have been proposed for operationally performing the requisite time averaging (Ishii 1975). Further concerns regarding the evolution of the system in relation
to the ergodicity of the spatial averages, often used to justify the time-averaging
hypothesis, are indicated in the work of Nitsche & Brenner (1989). However, even
such measure-preserving ®ows as those we envision in a unit cell (equivalent to ®ows
on a torus) can exhibit a wide range of droplet behaviour, including periodic evolution, almost periodic evolution, and even chaotic motion (Ottino 1989). Accordingly,
it is virtually impossible to state anything of general validity about the time dependence of the singular surface microscale ­ eld variables in spatially periodic porous
media (whose time dependence arises from the gross motions of the droplets relative
to ­ xed bed particles). Here, we make the following plausible assumption regarding the time dependence of the di¬use microscale ­ eld variables for a ®ow which is
speci­ ed to be macroscopically steady (albeit unsteady at the microscale level).
Assumption 3.1. Let f be any bounded function inside the unit cell, i.e. on a
suitable norm,
kf k C;
where C is a constant. We assume that 9 a T such that
1
kf (t + T ) ¡ f (t)k O(T ) 8 t:
T
(3.10)
Though formally stated above in a perhaps ponti­ cal aphysical manner, the fact
that the functions f are assumed to be dependent upon time only through the instantaneous location and con­ guration of the `droplets’ requires the latter to return to,
or at least approach very closely, their original con­ gurations after a time T . (Of
course, in using coarse-scale words like `droplets’ in our di¬use microscale description we are speaking loosely, since such terminology is strictly appropriate only at the
singular microscale view.) This time periodic (or almost time periodic (Bohr 1947))
assumption represents the dynamical constraint analogous to the geometrical constraint imposed by spatial periodicity. However, while an explicit length-scale exists
for the porous medium (namely L) to compare with l, the variable T presents no
such explicit choice in view of the steadiness of the macroscale ­ elds. Thus, T can be
considered in­ nitesimal compared with the relaxation time of the measuring instruments employed at the macroscale. Accordingly, the above assumption assures the
requisite smoothness (on the time-scales of measuring instruments) of the macroscale
temporal derivatives.
The quantity T , which is of course ­ nite at the microscale, will be identi­ ed later as
constituting an in¯nitesimal macroscale time element dt, whence the group of terms
appearing on the left-hand side of equation (3.10) will be identi­ ed as representing
the de­ nition of the macroscopic time derivative. Furthermore, the identi­ cation of
T de­ nes the time interval over which the physical variables need to be time averaged
in order to complement the spatial averaging. Both are needed to derive expressions
for the macroscale ­ elds from their microscale counterparts.
(c) Microscale conservation equations
This section discusses the conservation equations governing the microscale-level
(di¬use interface) transport processes.
Proc. R. Soc. Lond. A (2000)
762
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
At the pore-space (microscale) level of description our biphasic ®ows are modelled from a di¬use interface viewpoint. The conservation equations governing species
transport at the microscale correspond to the convective-di¬usion and linear momentum equations outlined in x 2 b. They are summarized here, along with the corresponding microscale transport equations for both energy and moment-of-momentum.
Di® usion equations:
@» 1
+ r (» 1 v) = ¡ r J1 ;
@t
@» 2
+ r (» 2 v) = ¡ r J2 :
@t
(3.11 a)
(3.11 b)
Momentum equations:
(a) Linear momentum.y
@(» v)
+ r (» vv) = ¡ rp + r ¿v + r ¿c +
@t
2
» i gi ;
(3.12)
i= 1
where gi is the external body force per unit mass exerted on species i.
d ef
(b) Moment-of-momentum. De­ ne ¿ = ¿v + ¿c and assume the ®uid to be nonpolar. In such circumstances we have that (Aris 1962)
¿ = ¿y :
(3.13)
Energy equation (deGroot & Mazur 1984):
@[» (e + 12 v 2 )]
+ r [» (e + 12 v 2 )v + q + vp ¡
@t
2
v ¿] ¡
(» i v + Ji ) gi = 0; (3.14)
i= 1
where e denotes the volumetric internal energy densityz and q the heat-®ux vector.
The pore-level boundary conditions are dictated by the geometry of the system.
At the surfaces of the bed particles we have (i) the no-slip condition,
v=0
on
sp ;
(3.15)
and (ii) the condition of no net normal ®ux, which together with the vanishing of
the convective ®ux on the bed particles as implied by (3.15) requires that
º Ji = 0
on
sp
(i = 1; 2);
(3.16)
with º the unit normal on the bed particles. Consistent with the di¬use interface
model, no interfacial boundary conditions exist.
y To be consistent with the spatially periodic interstitial-®uid ­ eld model of a porous medium here
envisioned, gravity forces cannot strictly be included within this framework since the concentration of
the individual species (and hence » i ) vary over the unit cell. However, as will be made clear during
the subsequent exposition, the spatially periodic model serves as a backbone for constructing rigorous,
physically based de­ nitions of the macroscale ­ elds in terms of the microscale ­ elds. The ­ nal forms of
the macroscale transport equations, namely (3.53){(3.56), are expected to be independent of this strict
periodicity assumption. In line with such an approach, we ignore the inconsistency arising from the
inclusion of gravity forces within the spatially periodic model. Alternatively, the inconsistency can be
removed by setting the gravity terms to zero.
z For liquid mixtures as considered in this article, the di¬erence between the internal energy density
and the free energy density is expected to be negligible, thereby requiring that the volumetric density
» e be identical to f 0 , whose de­ nition was outlined earlier in the footnote on p. 745.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
763
(d) Macroscale de¯nitions
(i) Time-scales
As already observed in connection with the temporal evolution of the microscale
­ elds accompanying the ®ow, the relative motion of the two `phases’ leads to an
inherently unsteady phenomenon at the di¬use (as well as the singular) microscale|
this despite the fact that the macroscale ®ow may itself be steady, that is time
independent! In such circumstances one can identify the concept of a micro{macro
time interval in accordance with the assumption of time-periodic or almost timeperiodic ®ows, with T (assumption 3.1) constituting an in­ nitesimal time-interval dt,
say, at the macroscale. Consequently, time-intervals of order T will be regarded as
being di¬erentials at the macroscale. The macroscale manifestation of the underlying
microscale phenomena is observed only in a time-averaged sense at the macroscale,
consistent with the continuum description. Therefore, we write
dt = T;
(3.17)
where t constitutes the macroscale time. In analogy with the previous de­ nition of
macroscale lengths following (3.3), one can heuristically write t = t + O(T ). Furthermore, all subsequent de­ nitions of macroscale ­ elds will involve a time averaging of
the comparable microscale ­ elds of the generic form
1
T
t+ T =2
(: : : ) dt:
(3.18)
t¡T =2
For simplicity of notation, we henceforth use the succinct abbreviation
1
T
t+ T =2
(: : : ) dt ²
t¡T =2
(
)
(3.19)
t
to denote the time-averaging operator.
The subsequent analysis focuses primarily on ®ows that are steady at the macroscale (i.e. satisfy assumption 3.1). The ensuing macroscale equations are then generated by appropriately space- and time-averaging the microscale equations. In view of
the assumed spatial- and time-periodicity (or, less stringently, almost time-periodicity) of the microscale ­ elds, space-averaging is equivalent to averaging over a unit
cell ½ o , and time-averaging to averaging over one period T . Quantities of O(klk) and
O(T ) are subsequently neglected,y consistent with the identi­ cation of such quantities as being in­ nitesimal at the macroscopic level. While the main focus of our work
will be on ®ows that are steady at the macroscale, we will also indicate the manner
in which the de­ nitions of the macroscale ­ elds can be extrapolated to include timedependent macroscale ­ elds as well, albeit only those macroscale ­ elds which vary on
y
In the case of almost time-periodic functions, the de­ nition (Bohr 1947)
lim
T0 ¡ !
1
1
T0
t+T 0
d ef
f (t) dt = f^;
t
of f^ is such that the convergence to f^ occurs uniformly in t. As such, one can identify a ­ nite timeinterval T su¯ cient to de­ ne a mean value of f that is independent of t to within an error of O (T ). This
permits one to de­ ne macroscopically steady values for almost time-periodic microscopic ­ elds, since
quantities of O (T ) are considered in­ nitesimal at the macroscale.
Proc. R. Soc. Lond. A (2000)
764
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
a much longer time-scale than their microscale counterparts. In such circumstances
the above de­ nitions can be assumed to hold at a particular macroscopic instant
t, with the averages allowed to be functions of the macrotime t. Similar generalizations may be applied to slowly varying spatial phenomena, thereby allowing for the
possibility of spatially inhomogeneous macroscale ­ elds.
Notation. We distinguish macroscale ­ elds and operators from their microscale
counterparts by a¯ xing an overbar to the latter. Examples include the macroscale
position vector R, the phase-speci­ c pressures p1 , p2 , and the macroscale gradient
operator r ² @=@ R.
(ii) Macroscale physical quantities
Macroscale physical quantities are de­ ned `pointwise’ at the position R and `instantaneously’ at the time t as being the respective volume averages of comparable microscale physical quantities. With à representing the volumetric density of a
generic macroscale physical property, the de­ nition of the corresponding interstitial
macroscopic density is (cf. (3.48))
1
d ef
à =
½
d3 r Ã:
f
t
½
(3.20)
f
By integrating only over ½ f we have implicitly assumed that the bed particles are
impermeable to transport through their interiors of the physical property in question.
The de­ nition (3.20) is consistent with the usual viewpoint adopted by a macroscopic
observer in representing the total amount, t ½ f d3 r Ã, of a `physical quantity’ instantaneously present in a macroscopic di¬erential (interstitial) volume as the product
of a macroscopic volumetric density à and the magnitude ½ f of the `di¬erential’ volume element. However, in view of the fact that only a fraction of the pore space is
occupied by interstitial ®uid, the volumetric density can also be de­ ned equivalently
to (3.20) as
"Ã =
1
½
d3 r Ã;
o
t
½
(3.21)
f
where " ² ½ f =½ o represents the interstitial void fraction or porosity of the porous
medium. In the standard parlance of ®ow through porous media (Brenner 1972;
Brenner & Edwards 1993), Ã represents the interstitial volumetric density and "Ãy
the super­ cial volumetric density of the macroscale property in question.
The generic de­ nitions outlined in x 3 d lead us to identify the following physical
di¬use Darcyscale densities.
Mass density:
1
d ef
"» =
½
o
d3 r » :
t
½
(3.22)
f
y In the case where à is assumed to be spatially non-uniform (possessing weak gradients) the above
de­ nition of the macroscale physical quantities, together with the de­ nition of their respective macroscale
gradients, cf. equation (3.49), requires l to be small compared with the length-scale of the second derivatives in order to justify the neglect of second derivatives.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
765
Mass density of species 1:
"»
1
d ef
1
=
½
o
d3 r » 1 :
t
½
(3.23)
f
The latter enables us to de­ ne a macroscale volume fraction (or saturation) of phase 1
as
1
d ef
"c =
so that »
1
= » 01 c and »
2
= » 02 (1 ¡
½
o
d3 r c;
t
½
c) ² » 02 (1 ¡
» =»
=»
(3.24)
f
c). Hence, we ­ nd from (2.5) that
1+» 2
0
0
1 c + » 2 (1
¡
c);
(3.25)
consistent with expectations. In view of the assumed incompressibility of the mixture
at both the micro- and the macroscales, it is not surprising that the conventional
de­ nition of saturation appears as the natural choice for the physical composition
variable in our continuous multicomponent description (of what is normally regarded
as a two-phase system).
The de­ nition of macroscale momentum density leads to the following de­ nition
of the macroscale interstitial velocity v of the binary mixture (with " v representing
the so-called super­ cial or seepage velocity of the mixture):
1
d ef
"» v =
½
o
d3 r » v:
t
½
(3.26)
f
The velocity v appearing in the above expression also constitutes the mass-average
velocity (cf. Appendix D). In contrast with conventional two-phase theories of averaging, involving a singular interface Darcyscale description of the pertinent phenomena
(Marle 1977; Whitaker 1986), it will prove unnecessary for us to de­ ne individual
phase velocities v 1 and v 2 at the macroscale. This innovation is consistent with our
di¬use Darcyscale model, which uses a multicomponent (rather than multiphase)
viewpoint accounting for the several species present, and where accordingly the only
relevant velocity is the mixture velocity. Nevertheless, for the sake of comparison
with existing theories of multiphase ®ow, individual macroscale di¬usion ®uxes of
the two species will be de­ ned in the next section, enabling us to indirectly identify
species velocities (equivalent to so-called phase velocities). In turn, this will enable us
to relate our ­ ne-scale di¬use Darcyscale description to its more traditional, coarserscale, singular Darcyscale counterpart.
It is important to emphasize the fact that the velocity v plays a dual role as both
(i) the momentum density (per unit mass) in a dynamical context; and (ii) the mass
®ux vector in a kinematical context. Hence, our subsequent generic de­ nition of the
macroscopic ®uxes in (3.32), when applied to the mass ®ux, should (and does) furnish
an expression identical to (3.26) for the macroscopic velocity (see Appendix D).
Angular momentum density
For a polar microscale interstitial ®uid continuum the total angular momentum
density (per unit mass) will generally include an intrinsic pseudovector contribution
Proc. R. Soc. Lond. A (2000)
766
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
l from the internal angular momentum in addition to the usual contribution arising
from the moment of the linear momentum (Dahler & Scriven 1963). However, it is
a consequence of our assumption regarding the non-polar nature of the ®uid at the
microscopic level, that l = 0. Nevertheless, the macroscale intrinsic angular momentum density will generally embody within it a non-zero contribution arising from
external moment-of-momentum considerations in relation to the screw-symmetric
nature of the particles comprising the skeletal porous medium (Brenner 1970, 1972).
Thus,
1
" » (l + R £ v) =
½
d3 r » R £ v:
o
t
½
(3.27)
f
Together with the de­ nitions R = Rn + r and Rn ² R, the latter reduces to the
following expression whereby l can be calculated from the microscale data as
1
d ef
"» l =
½
d3 r » r £ v:
o
t
½
(3.28)
f
From the latter de­ nition of l one can discern that it is O(klk). Accordingly, to
maintain consistency with the subsequent analysis l will henceforth be neglected.
The macroscale energy density e per unit volume is de­ ned by the relation
" » (e + 12 v 2 ) =
1
½
o
t
½
d3 r » (e + 12 v 2 );
(3.29)
f
which furnishes the following explicit expression for the internal energy density e:
1
d ef
"» e =
½
o
t
d3 r » (e + 12 v v 0 ):
½
(3.30)
f
The latter integrand includes contributions from the microscale internal energy density e, as well as a kinetic energy contribution arising from the distinction between
1 2
v and 12 v 2 , with v 0 = v ¡ v.y Furthermore, for reasons identical to those eluci2
dated in the second footnote on p. 762, the di¬erence between the macroscale free
energy density and the macroscale internal energy density is expected to be negligible. Thereby, we obtain
1
d ef
"F =
½
o
t
½
d3 r » (e + 12 v v 0 );
(3.31)
f
wherein F denotes the macroscale volumetric free energy density.
(iii) Macroscopic ° uxes
The following generic de­ nition is adopted for super­ cial macroscale areal ®ux
densities.
y This phenomenon, whereby the clear distinction between `internal’ and `external’ energies at one
scale becomes blurred at a coarser scale, is analogous to comparable phenomena at the pre-continuum,
discrete molecular scale, whereby, say, the kinetic energy of discrete molecules becomes a part of the
internal energy of the system at the continuum scale (so that the sum of the kinetic energies of the
molecules is not identical with the kinetic energy of the continuum) (Kirkwood 1967).
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
767
De¯nition 3.2. At a point R of the porous medium, the macroscale ®ux-density
tensor M ² M (R)z of the comparable microscale ®ux density tensor M ² M (R; t)
is de­ ned generically as
d ef
ds M :
sk M =
t
(3.32)
s k fn g
Inasmuch as sk ² ds, the preceding expression constitutes a Eulerian de­ nition
of the macroscale ®ux provided that M can be shown to be independent of the
orientation sk of the cell face. This fact is demonstrated by Brenner (1972) (see
also Nitsche & Brenner (1989)) for spatially periodic ­ elds M , albeit without the
temporal integration. This Eulerian de­ nition eventually enables us to satisfy the
physical requirements demanded of any proposed de­ nition of a macroscale ®ux
density, e.g. the stress tensor, since in this case, the right-hand side of the above
equation respectively represents the force exerted on a mass ®owing across an area
which is `in­ nitesimal’ at the macroscale. The following results, proofs of which are
furnished in Appendix C, are consequences of the de­ nition (3.32).
Theorem 3.3. When the instantaneous gradient rM of the microscale °ux density M is spatially periodic, the macroscale °ux density M (R), de¯ned so as to
possess the areal property (3.32) (albeit only to O(klk) terms), is explicitly given by
the expression
M =
1
½
r ds M :
o
t
@½
(3.33)
o
Moreover, M de¯ned in this way possesses a macroscale gradient, expressible as
rM =
1
½
o
t
@½
o
t
@½
o
ds M :
(3.34)
ds M :
(3.35)
Upon scalar contraction the latter yields
1
r M =
½
o
The preceding properties possessed by M enable us to explicitly identify the following super­ cial areal ®ux densities, each respectively associated with a di¬erent
choice of physical property.
Mass ° ux-density vector:
1
d ef
JM =
½
o
» r ds v:
t
@½
(3.36)
o
Proof of the fact that the above de­ nition is consistent with one’s expectation that
J M should be expressible in terms of the macroscale variables » and v as J M = " » v
is presented in Appendix D.
z
Note that M ² M (R ; t) in the case of a ®ow varying temporally on the macroscopic time-scale.
Proc. R. Soc. Lond. A (2000)
768
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
Di® usion ° ux-density vector:
1
d ef
N1 =
½
r ds N1 ;
o
t
@½
(3.37)
o
where N1 = » 1 v + J1 . Here, N 1 denotes the net macroscale ®ux-density of species 1
(relative to the stationary porous medium) arising from both microscale convection and di¬usion processes. Analogous to the comparable microscale decomposition
appearing in the di¬use interface description, one can separate N 1 into a sum of
convective and di¬usive ®uxes:
N 1 = " » 1v + J 1;
(3.38)
where the macroscale di¬usion ®ux J 1 is de­ ned as
1
d ef
J1 =
½
o
r ds (» 1 v 0 + J1 );
t
@½
(3.39)
o
and v 0 is de­ ned following (3.30), being the deviation from the mean velocity. Similarly, one can de­ ne the macroscale momentum ° ux-density dyadic
d ef
PM = "» vv + P;
with the macroscale stress tensor P de­ ned as
d ef 1
r ds (» vv 0 + P );
P =
½ o t @½ o
(3.40)
(3.41)
in which P ² ¡ pI + ¿ denotes the microscale stress tensor.
Appendix E discusses the manner in which P as de­ ned above can be expressed
as the sum of a macroscale pressure p and a (generally asymmetric) macroscale
deviatoric stress tensor ¿ :
P = ¡ pI + ¿ ;
(3.42)
wherein
d ef
p = ¡
1
3½ o
r ds (» vv0 + P )
t
@½
(3.43)
o
and
1
d ef
¿ = ¡
½
o
(rI ¡
t
@½
1
Ir)
3
ds (» vv0 + P ):
(3.44)
o
As further discussed in Appendix E, it is the macroscopic pressure gradient that is
physically signi­ cant, rather than the macroscale pressure itself. Indeed, according
to (E 8),
rp =
½
1
o
ds p:
t
@½
(3.45)
o
Similarly, the net macroscale energy ° ux-density vector is de­ ned as
d ef
J e = " » (e + 12 v 2 )v ¡
Proc. R. Soc. Lond. A (2000)
P v + q;
(3.46)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
769
with
1
d ef
q =
½
o
t
@½
r ds [q + » (e + 12 v 2 )v 0 ¡
P v0 ]
(3.47)
o
identi­ ed as the macroscopic heat ®ux-density vector. This completes the de­ nitions
adopted for the di¬use Darcyscale ®uxes pertinent to our analysis.
(iv) Macroscopic source terms
Macroscale source terms are de­ ned as being the volume averages of the respective
microscale source terms. Again, this is consistent with measurements that would be
registered by a macroscale experimental probe; that is, a probe whose aperture is of
the order of the cell size.
Upon using equation (3.23) de­ ning " » i , the momentum source term is given by
2
" » i gi :
i= 1
The macroscopic energy dissipation due to the di¬usive ®ux is given by
2
iv
(" »
+ J i ) gi ;
i= 1
where the de­ nition (3.39) of J i has been used.
To e¬ect the averaging procedure enabling the transition from microscale to macroscale we require use of the following generic de­ nition of the macroscale gradient.
(v) Macroscopic gradients
De¯nition 3.4. For a function of volumetric density à whose macroscale value
is de­ ned as
1
d ef
"Ã =
½
o
d3 r Ã;
t
½
(3.48)
o
the corresponding macroscopic spatial gradient is de­ ned as
1
d ef
rà =
½
f
ds Ã;
t
@½
(3.49)
o
and the macroscale time derivative as
@ Ã d ef ¡1 1
= T
½ f
@t
½
d3 r Ã(t + 12 T ) ¡
o
½
1
f
d3 r Ã(t ¡
½
1
T)
2
:
(3.50)
o
The latter is consistent with our time-element identi­ cation in (3.17). Functions
à satisfying assumption 3.1 possess a macroscopic time derivative that is identically
zero:
@Ã
= 0:
@t
Proc. R. Soc. Lond. A (2000)
(3.51)
770
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
However, in the general case where the ®ow is macroscopically unsteady, (3.50) translates as
1
@Ã
=
½ f
@t
d3 r
t
½
o
@Ã
:
@t
(3.52)
The above identi­ cations permit derivation of the pertinent ­ eld equations governing transport at the macroscale starting from the corresponding microscale transport
equations, a topic addressed in the following section.
(e) Macroscale equations
Averaging the microscale equations over time and space (the former over one timeperiod T and the latter over the unit cell ½ o ), and using the preceding de­ nitions of
the macroscopic time derivative together with the expressions for the macroscopic
®ux densities and their divergences jointly with the boundary conditions (3.15),
(3.16), yields the following di¬use Darcyscale ­ eld equations.
Di® usion equations:
@(" » 1 )
+ r (" »
@t
@(" » 2 )
+ r (" »
@t
1
v) = ¡ r J 1 ;
(3.53 a)
2
v) = ¡ r J 2 :
(3.53 b)
Momentum equation:
@(" » v)
+ r (" » v v) = ¡ rp + r ¿ +
@t
2
" » i gi + F^ ;
(3.54)
i= 1
where
d ef
F^ =
ds P
t
(3.55)
sp
denotes the macroscopic volumetric force density arising from the immobility of the
bed particles.y
Energy equation:
@[" » (e + 12 v 2 )]
+ r [" » (e + 12 v 2 )v ¡
@t
2
P v + q] ¡
(" »
i
v + J i ) gi = 0: (3.56)
i= 1
y The genesis of F^ lies in the no-slip condition on the solid particle surfaces. Slip or other forms
of boundary conditions resulting from the presence of contact lines would presumably lead to di¬erent
boundary conditions at the solid surfaces, and hence would also change the constitutive form of the force
appearing in the ­ nal expression. However, there does not exist an analogous term in the macroscale
energy equation (3.56) even for varying temperature ­ elds. The underlying reason for this is that the
force F^ is of a physicochemical nature, arising from the stationarity of the particles, rather than being
externally imposed. Were the latter the case, this force would then have necessarily appeared in the
energy equation. This subtle feature plays a crucial role in determining the form of Darcy’s law from
considerations of irreversible thermodynamics, as is outlined in x 3 f.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
771
Table 2. Identi¯cation of the respective forces and ° uxes in the entropy source terms appearing
in the generalized entropy balance
physical entity
di® usion ° ux
body-force density
deviatoric stress (symmetrized)
heat ° ux
generalized force
generalized ° ux
Jk
F^
¿s
q
r· k ¡ g k
v ¡U
D
rT
All quantities appearing in these di¬use Darcyscale ­ eld equations have been rigorously de­ ned so as to maintain consistency with the physical measurements discernible to a macroscale observer. Moreover, they are operationally expressed entirely
in terms of spatial and temporal quadratures of well-de­ ned microscale ­ elds. As
the subsequent analysis will be primarily concerned with inertia-free and isothermal
microscale ®ows, e¬ects arising from the existence of the macroscopic stress tensor
as well as from energy considerations are consequently irrelevant and hence will not
be further pursued except insofar as they prove necessary in our generic discussion
of macroscale irreversible thermodynamic principles.
(f ) Irreversible thermodynamics and macroscopic constitutive equations
Irreversible thermodynamics provides a convenient formalism for establishing the
functional forms of the macroscale constitutive equations required in the macroscale
­ eld equations, at least in circumstances where these constitutive relations are linear. The methodology employed involves the identi­ cation of forces and ®uxes in
the entropy source term of the generalized entropy balance, subsequent to which
one proposes a general linear relationship between these forces and ®uxes (deGroot
& Mazur 1984). The fundamental basis for constructing such an entropy balance
lies in the assumption of local thermodynamic equilibrium, which posits that thermodynamic identities valid at equilibrium continue to remain locally valid in nonequilibrium circumstances. Although one can thereby identify the constitutive forms
of the pertinent linear force{®ux relationships within the framework of linear irreversible thermodynamics, the formalism itself provides no insight into the magnitudes
of the phenomenological coe¯ cients appearing therein (with the possible exception
of their algebraic signs), including whether coupling occurs between the di¬erent
®uxes owing to the presence of non-zero phenomenological `cross-coupling’ coe¯ cients. Furthermore, it is not possible to predict the functional dependence, if any, of
the phenomenological coe¯ cients upon the (macroscale) physical properties. Despite
such limitations, the success of linear irreversible thermodynamics in providing valuable insights into the possible forms of constitutive equations cannot be denied.
In this section we apply such concepts to the dependent macroscale ­ eld variables
quantifying our `two-phase’ system.
In view of the relative simplicity of the preceding system of di¬use Darcy macroscale equations, the generalized entropy balance equation is easily obtained along
the lines outlined by deGroot & Mazur (1984). For brevity we sketch only essential
details.
Proc. R. Soc. Lond. A (2000)
772
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
The kinetic energy balance is obtained from the momentum equation (3.54) as
@( 12 " » v 2 )
+ r ( 12 " » v 2 v) = (r P ) v + F^ v +
@t
2
"»
i
v gi :
(3.57)
i= 1
Upon performing an analysis predicated along the same lines as outlined by deGroot
& Mazur (1984), the following entropy transport equation is obtained (where the
macroscopic volumetric entropy density s is de­ ned as usual in terms of macroscopic
thermodynamic quantities as T ds = de + p dv;y with v = 1=» ):
d(" » s)
=¡ r
T
dt
2
q¡
·
1
q rT
T
¡
kJ k
k= 1
2
J k ( r·
¡
k
¡
gk ) ¡
F^ v ¡
¿ s : D;
(3.58)
k= 1
wherein d=dt ² @=@ t + v r represents the macroscale material derivative, and · k
the macroscale chemical potential of species k. Furthermore, D denotes the macroscopic version of the deformation tensor (cf. (2.19)), and ¿ s the symmetric part of
the macroscopic deviatoric stress tensor. As the distinction between the microscale
temperature T and the macroscale temperature T is peripheral to the dynamical
issues of interest, it is not further addressed here, whence we denote the macroscale
temperature also by T .z Consistent with the assumptions of linear response theory, the respective forces and ®uxes can be identi­ ed from the entropy source term
(table 2). The above list constitutes only the irreversible contributions to the ®uxes.
To remind readers of the reference frame in which v is measured, these ®uxes include
the translational velocity U of the porous medium as a whole (measured relative to
an inertial reference frame). This also renders the constitutive equations objective
(Truesdell 1966).
Based on the above identi­ cations, together with the assumed absence of coupling
between ®uxes of tensorial orders di¬ering by an odd integer, which holds true for a
porous medium composed of centrally symmetric particles (Nye 1985), the following
macroscopic constitutive equations are obtained:
J k = ¤Ik (r·
k
2
F^ = ¤III v +
¡
gk ) + ¤II
k v
(r·
¤IV
k
k
¡
gk )
(k = 1; 2);
(k = 1; 2);
(3.59)
(3.60)
k= 1
¿ = ¤V : rv:
(3.61)
y Strictly speaking, the mean pressure p is not the same as the thermodynamic pressure, whence
one should de­ ne P = ¡p T I + ½ 0 I + ¿ , where p T now refers to the thermodynamic pressure, with
the di¬erence p ¡ p T denoted as ½ 0 . The subsequent analysis requires knowledge of only that portion
of P excluding the thermodynamic pressure pT . Since the diagonal tensor ½ 0 I contributes only to the
symmetric part of P , we combine ½ 0 I with ¿ s and denote the sum also by ¿ s . However, the distinction
between the thermodynamic and the mean pressures vanishes when r v = 0 (Brenner & Nadim 1996),
which condition applies to the class of ®ows eventually considered here.
z At issue here is the question of whether or not a continuum that is isothermal at the macroscale is
also necessarily isothermal at the microscale, and, if so, is T = T ?
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
773
(Here, in contrast to the tabulation in table 2, we have reverted to our earlier convention of measuring velocities relative to the ­ xed porous medium, whence U has been
set equal to zero.) The individual di¬usive ®uxes J 1 and J 2 are not independent variII
ables since J 1 + J 2 = 0, whence only two of the four coe¯ cients ¤I1 ; ¤I2 ; ¤ II
1 ; ¤2 are
independent. Furthermore, invocation of the Onsager reciprocity principle, assumed
applicable at the macroscale (deGennes 1983; Kalaydjian 1993), leads to the observaIII and ¤IV (k =
tion that only three of the phenomenological coe¯ cients, ¤Ik ; ¤II
k;¤
k
1; 2), are independent.
A major limitation of irreversible thermodynamics lies in its inability to provide
any indication of the magnitudes of these phenomenological coe¯ cients, except possibly their algebraic signs. Moreover, the above relations are purely phenomenological,
and hence do not constitute a self-contained theory. As a consequence, their applicability to the phenomena at hand needs to be rati­ ed either through comprehensive
experiments or a detailed microscale theory, or both. While the former veri­ cation
mode is left to experimentalists, x 3 m provides an illustrative analysis (albeit in the
somewhat restricted context of quasi-static microscale behaviour) of the manner in
which the second mode of con­ rmation may be implemented.
Under the dual assumptions of slow ®ows and weak macroscale gradients in the
mean velocity v, we may neglect e¬ects arising in the macroscale equations from
(i) the presence of ®uid inertia; (ii) the existence of a deviatoric stress tensor. Furthermore, neglect of coupling between the ®uxes J k and the force F^ leads, upon
substituting the above relation for F^ into (3.54), to the relation
2
rp = ¤III v +
»
i gi ;
(3.62)
i= 1
which constitutes a di¬use, macroscale, single-phase, multispecies Darcy’s law. In
the case where g1 = g2 ² g, say, where g is the acceleration of gravity, we ­ nd upon
inverting equation (3.62) that
v = (¤III )¡1 (rp ¡
» g):
(3.63)
It is not surprising that through linearized irreversible thermodynamic arguments
we have obtained a macroscale form of Darcy’s law, which is itself a linear constitutive relation. A more rigorous derivation of (3.63) from ­ rst principles, is presented
in x 3 m. Previous studies by others (Mokadam 1961; Kalaydjian 1993) have also indicated how Darcy’s law may be obtained by using constitutive relationships derived
from irreversible thermodynamics or some other basis. At a later stage (see equations
(3.73), (3.74)) we address the more general case where coupling between the ®uxes
is included.
Compared with earlier work by others, our main accomplishments thus far have
been (i) rigorous de­ nitions of the di¬use Darcyscale ­ eld quantities, each of which
was chosen to accord with its physically agreed-upon interpretation; and (ii) subsequent use of these de­ nitions to derive Darcyscale conservation equations. Such
rigorous de­ nitions were lacking in a number of previous studies also concerned with
the (albeit, singular) Darcyscale viewpoint. In addition, we have outlined a hitherto
unproposed scheme for viewing the ®ows from a di¬use rather than singular interface
viewpoint. This naturally led to a `two-component mixture’ macroscale formulation,
and eventually to the above single-phase Darcy’s law relating the total pressure gradient rp to the mixture velocity v. As the resulting system of equations is closed
Proc. R. Soc. Lond. A (2000)
774
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
with respect to the numbers of variables and equations, classical closure concepts
like capillary pressure do not explicitly arise at this level of description.
Sections 3 g{m, which follow, are devoted to clarifying the concept of phase-speci­ c
Darcy’s laws in an e¬ort to relate our work to more conventional singular, Darcyscale
formulations of two-phase ®ow phenomena. By so doing, our two-component mixture
framework renders rational the previously empirical two-phase Darcy ­ eld equations.
Concomitantly, the exercise furnishes valuable insights into the concept of capillary
pressure as a macroscale variable.
(g) Phase-speci¯c quantities
This section provides an interpretation of the conventional empirical singular
Darcyscale phase-speci­ c quantities appearing in the ubiquitous two-phase Darcy
®ow equation (1.2). This is achieved by correlating these phase-speci­ c ­ elds with
our two-component microscale quantities. The basis for accomplishing this lies in the
di¬use two-component mixture formulation that we have adopted at both the microand macroscales. The approach used in this section involves relating phase-speci­ c
(or, more properly, species-speci­ c) macroscale ­ elds appearing in the singular Darcyscale description (1.2) to their di¬use Darcyscale counterparts. Thus, in contrast
with many prior analyses (Marle 1977; Whitaker 1986; Kalaydjian 1993), singular
Darcyscale phase-speci­ c quantities will not be de­ ned as averages of their (singular
interface) microscale counterparts|so that, for instance,
1
p1 6=
½
o
d3 r p 1 :
t
½
(3.64)
o
Only inertialess (slow) isothermal ®ows are addressed in the ensuing discussion.
(i) Phase velocities
The di¬use Darcyscale formulation that we have adopted serves to de­ ne v and
J i in a rigorous manner. A macroscale observer then discerns the individual `phase
velocities’ v 1 and v 2 only indirectly, namely from measurements of J 1 and J 2 , by
pursuing the following argument. By de­ nition,
J 1 = » 1 (v 1 ¡
J 2 = » 2 (v 2 ¡
v);
v):
(3.65)
Consequently,
v1 =
J1
+ v;
» 1
v2 =
J2
+ v:
» 2
(3.66)
It is readily established that the above de­ nitions imply that
" » 1 v1 =
½
1
o
d3 r » 1 v 1 ;
t
½
(3.67)
f
together with a comparable relation for species 2. However, our interpretation of v 1
and v 2 is independent of (3.67) et seq., since at the di¬use interfacial microscale only
the mixture velocity v is physically relevant in the context of our model.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
775
(ii) Phase-speci¯c pressures
Before dealing with the concept of phase-speci­ c pressures (in the context of the
singular Darcy description), and as a preliminary to the subsequent development, we
de­ ne the following di¬use Darcyscale quantities.
Macroscale chemical potential. Consistent with the assumption stated in the Introduction, we hypothesize that conventional thermodynamic identities among microscale quantities hold equally well when macroscale thermodynamic quantities (such
as free energy) are involved. A similar assumption has been made by others (Marle
1977; Pavone 1989; Kalaydjian 1993), who likewise employed conventional thermodynamic identities at the macroscale. Thus, we de­ ne the macroscale chemical potentials:
d ef
·
1
=
@F
;
@» 1
·
@F
:
@» 2
d ef
2
=
(3.68)
The indicated di¬erentiations are to be performed while keeping constant all other
physical variables other than species density.
Phase-speci¯c pressures. With the above considerations in mind we propose the following de­ nitions for the singular Darcyscale phase-speci­ c pressures:
rp 1 =
»
1
»
r·
1;
»
rp 2 =
2
»
r·
2:
(3.69)
These relations are exact for homogeneous (and weakly inhomogeneous) systems at
thermodynamic equilibrium. Their extension to non-equilibrium situations, as in our
case, is linked to the assumption of local equilibrium|discussed in x 3 f. In fact, Kirkwood & Bearman (1958) in their discussion of transport processes in multicomponent
mixtures propose identical relations, albeit for species-speci­ c pressures. The above
de­ nitions do not incorporate e¬ects arising from incompressibility. This restriction is incorporated into our work by a novel method ­ rst proposed by deGennes
(1980), although in a di¬erent context. Its implementation involves assuming that
the incompressibility condition can be viewed as resulting from the action of ­ ctitious physicochemical forces, represented here by the gradient of some macroscale
^ . Explicitly,
potential U
r p1 =
»
1
»
r(·
1
+ U^ );
r p2 =
»
2
»
r(·
2
+ U^ );
(3.70)
^ is determined by imposing the incompressibility condition:
where U
rp1 + rp2 = rp:
(3.71)
^ in the above equations and substituting the resulting expression
Upon solving for U
into (3.70) one obtains
rp 1 =
»
»
1
rp +
» 1» 2
r(·
» 2
1
¡
·
2 );
r p2 =
»
»
2
rp +
» 1» 2
r(·
» 2
2
¡
·
1 ):
(3.72)
The latter pair of equations constitute the operational de­ nitions of the singular
Darcyscale phase-speci­ c pressures. It is appropriate here to review the evolution
Proc. R. Soc. Lond. A (2000)
776
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
that led to these de­ nitions. In particular, these pressures have not been obtained
by averaging respective phase-speci­ c microscale pressures p1 and p2 arising in the
singular interface view, as has been done in the past by other researchers. Indeed,
such individual pressure ­ elds have no existence in our di¬use, strongly inhomogeneous, single-phase view. Rather, in our di¬use microscale view, the pressure p of the
mixture (`solution’) as a whole was identi­ ed by its contribution to the linear momentum transport process. Rigorous de­ nitions of the di¬use Darcyscale ®uxes (in terms
of the microscale ®uxes) subsequently yielded the di¬use Darcyscale pressure of the
mixture. Within this mixture, individual species (or `phase-speci­ c’) macroscale pressures pi were identi­ ed with the corresponding gradients r· i in chemical potential,
consistent with the local equilibrium assumption. Incompressibility was incorporated
through ­ ctitious physicochemical forces. This conceptual foundation underlying rp1
and rp2 is subsequently used in the next section to explain the phase-equilibriumtype of relationship postulated (Marle 1981; Hassanizadeh & Gray 1993) to exist
between rp1 and rp2 .
Observe that we have not concerned ourselves here with the concept of phasespeci­ c stress tensors. Discussion of such quantities can be found, for example, in
the work of Kirkwood & Bearman (1958).
(h) Phase-speci¯c Darcy’s laws
In this section we use the above de­ nitions, jointly with the constitutive equations obtained in x 3 f, to derive phase-speci­ c Darcy laws underlying the singular
Darcyscale description of multiphase ®ows. As such an exercise entails manifold algebraic manipulations, pertinent details are relegated to Appendix F. The ­ nal results
obtained are
v 1 = (¤ 011 )¡1 (rp1 ¡
v2 =
(¤ 021 )¡1
(r p 1 ¡
" » 1 g 1 ) + (¤012 )¡1 (rp2 ¡
" » 1g1 ) +
(¤022 )¡1
( rp 2 ¡
" » 2 g 2 );
(3.73)
" » 2 g 2 );
(3.74)
in which the dyadics (or, equivalently, 3 £ 3 matrices) ¤0ij are tensorial phenomenological coe¯ cients. From our analysis it was established that only three of the above
four tensor coe¯ cients ¤0ij (i; j = 1; 2) are independent. However, we were unable
establish any symmetries of the cross-coupling tensors (deGennes 1983; Kalaydjian
1993).
The forms of the pertinent constitutive equations obtained from irreversible thermodynamic considerations, supplemented by our de­ nitions of phase-speci­ c quantities, thus yield the generalized two-phase form of Darcy’s law, namely equations
(3.73){(3.74), hypothesized by a number of other researchers (deGennes 1983; Pavone
1989; Kalaydjian 1993). Fundamental to these singular Darcyscale laws is, however,
the precursor multicomponent mixture, single-phase, Darcy’s law (3.63) that we
derived. This approach contrasts with that of other investigators, who treat the
phase-speci­ c Darcy’s laws as fundamental entities, necessitating a closure relationship for the capillary pressure Pc [c] (Auriault & Sanchez-Palencia 1986; Auriault
1987). Since our preceding phase-speci­ c relationships were developed from mixturebased laws, we are as a consequence able to provide rational and explicit expressions
for the capillary pressure as well as other elements appearing as ad hoc continuummechanical identities in singular Darcyscale theories (cf. x 3 i ).
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
777
Note. In this and subsequent sections we will be concerned exclusively with slow,
isothermal ®ows characterized by weak spatial gradients. The temperature dependence of the various physical quantities, especially the free energies, will not be
explicitly displayed unless essential to the development.
(i) Capillary pressure
As noted in x 1, the singular Darcyscale description of multiphase ®ows, being
based upon phase-speci­ c Darcy’s laws, requires knowledge of the so-called capillary
pressure Pc [c] to e¬ect closure of the resulting system of equations (Marle 1981; Hassanizadeh & Gray 1993). However, prior analyses have been unable to convincingly
identify the fundamental source of the hypothesized capillary pressure{saturation
relationship underlying such theories. While it has always been conjectured that the
capillary pressure{saturation relationship represents a macroscale manifestation of
the microscale singular-surface Laplace interfacial boundary condition (cf. (2.69))
(Whitaker 1986), attempts to relate quantities like p1 and p2 |resulting from volume averaging|to microscale matching conditions, which are valid only at interfaces
(Hassanizadeh & Gray 1993), are not wholly convincing. Moreover, the Laplace equation derives from equilibrium considerations, whereas the Darcy-®ow case necessarily
entails non-equilibrium processes! As such, prior theoretical arguments proposed for
the existence of a capillary-saturation relationship are at best conjectural and incomplete.
In contrast, the present work adopts an innovative formalism for describing microscale two-phase ®ows, where one no longer deals with intrinsically macroscopic terms
like `interfaces’. As a consequence, the Laplace boundary condition no longer manifests itself in an explicit manner, thereby enabling rigorous operational de­ nitions
to be set forth for the di¬use Darcyscale pressure ­ eld p in terms of quadratures
of microscale ­ elds. In turn, knowledge of this macroscale pressure, together with
use of other equally well-de­ ned macroscale ­ elds, has enabled us to derive singular
Darcyscale phase-speci­ c pressures p1 ; p2 through rigorous operational de­ nitions of
the latter. This led, in turn, to the expressions (3.72) for rp1 and rp2 .
Subtraction of the latter pair of equations yields
»
»
»
rp 1 ¡
»
1
rp2 = r(·
¡
1
·
2 ):
At equilibrium, when both r» 1 and r» 2 vanish identically, so that » 1 and »
hence » ) are each independent of position, it follows from the above that
»
»
(3.75)
2
p1 ¡
1
»
p2 = ·
»
1
¡
·
2:
2
(and
(3.76)
2
The left-hand side of the above equation constitutes our proposed de­ nition of the
heretofore elusive capillary pressure, namelyy
»
d ef
Pc =
»
p1 ¡
1
»
»
p2 :
(3.77)
2
y This pressure is de­ ned in a manner di¬erent from that advocated by others, who employ a volumeaveraging procedure (cf. the review by Adler & Brenner (1988)).
Proc. R. Soc. Lond. A (2000)
778
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
The de­ nition (3.77) constitutes the logical outcome of our designation of species(or phase-) speci­ c pressures, starting from the concept of the total pressure p of the
mixture. The above equation also embodies a de­ nitive scheme for calculating the
capillary pressure from the microscale data, since · 1 and · 2 are eventually de­ ned in
terms of their relationship to microscale quantities. We believe that our investigation
is the ­ rst to propose a rigorous, physically based, operational de­ nition of the
capillary pressure as a true continuum ­ eld concept. In a deeply intuitive study,
Hassanizadeh & Gray (1993) hint at similar considerations; but in view of the absence
of a suitable micro{macro homogenization scheme in their proposal, the physical
basis for such a hypothesis is lacking.
Until now our discussion was completely independent of the choice of the microscale constitutive equations adopted for J1 ; J2 and ¿c in equations (3.11 a) and (3.12).
In the following sections, namely xx 3 j {l, we adopt the constitutive equations outlined in x 2 c, enabling us thereby to provide a number of signi­ cant insights into
macroscale concepts which have heretofore de­ ed completely rational elucidation in
previous investigations. In view of our rigorous de­ nitions linking microscale and
macroscale physical quantities, we are able to furnish exact statements regarding the
functional dependencies of the macroscale physical quantities. These dependencies
constitute the appropriate form of an equation of state for macroscale-level modelling. We also display the complete set of equations required for macroscale modelling which embody the variables consistent with such an equation of state. An
experimentally determined equation of state, in conjunction with the complete set
of macroscale equations, thereby provides a framework for a macroscale modelling
scheme, one which is devoid of the need for actually computing the microscale solutions appropriate to the porous-medium geometry. Such a framework is developed
in the next three subsections.
(j ) Macroscale free energy density
Preliminary to seeking insights into the physical basis of capillary pressure, a
topic which is discussed in the next section, we ­ rst analyse the functional form of
the macroscale free energy density.
Calculation of the capillary pressure via (3.84) rests upon knowledge of the macroscale chemical potentials, · 1 and · 2 . However, the macroscale chemical potential is
calculated via (3.68) through the macroscale free energy density F , which is de­ ned
in equation (3.31):
1
d ef
"F =
½
o
t
½
d3 r » (e + 12 v v 0 ):
(3.78)
f
We consider the constitutive form proposed in footnotes on pp. 745 and 762 for the
de­ nition of the internal energy density e required above, thereby allowing us to
write (3.78) as
1
d ef
"F =
½
o
t
½
d3 r [f (c) + fid + 12 K(rc)2 + 12 v v 0 ]:
(3.79)
f
Note that the energy density F possesses contributions arising from three sources:
(i) microscale bulk internal energy;
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
779
(ii) microscale kinetic energy; and
(iii) gradient energy.
Contributions from sources (i) and (ii) are present in the macroscale energy density in
single-phase, binary mixture ®ows. However, the third contribution, representing gradient e¬ects, arises exclusively as a consequence of the presence of the second phase.
To facilitate the analysis of the respective functional dependence of the macroscale
free energy upon each of these contributions, we divide the contributions into two
groups, denoting by F b u lk the sum of the contributions (i) and (ii), and by F grad
that resulting from the gradient terms. Consequently,
"F b
u lk
=
1
½
o
t
½
d3 r (f + fid + 12 v v 0 ):
(3.80)
f
The functional dependence of each of the terms appearing above may be analysed
as follows.
(i) In view of the fact that f (c) = 0 for c = 0 and 1, f (c) is non-zero only
in the interfacial region. Furthermore, since f (c) is of O(1) everywhere (even
within the interfacial region, since there exists no explicit dependence of f
upon gradients) the contribution of f (c) to F b u lk is of O(¯ 1=3 ) (¯ denoting the
width of the interfacial region). As such, the contribution of this term proves
negligible compared with that of the other terms.
(ii) The contribution of fid ² cf1 + (1 ¡ c)f2 to F b
only upon c (in addition to temperature).
u lk
is functionally dependent
(iii) The contribution arising from the kinetic energy term can in®uence the macroscale free energy through the dependence of the latter upon the macroscale
temperature, in a manner not unlike that of the contribution of molecular
motion to the continuum-scale internal energy density. However, since subtle
di¬erences existing between the macroscale and microscale temperatures are
peripheral to the focus of this article we eschew this complication by assuming that the contribution arising from the microscale kinetic energy term is
negligible for the slow, low Reynolds number ®ows considered here.
The suggested functional dependencies of the individual contributions lead us to
conclude that F b u lk ² F b u lk (c; T ).
Consider the gradient contribution to the macroscale free energy:
" F grad =
½
1
o
d3 r » egrad =
t
½
f
½
1
o
t
½
d3 r 12 K(rc)2 :
(3.81)
f
By de­ nition, egrad is non-zero only in di¬use interfacial regions. Using (2.70), relating the interfacial surface tension to the concentration gradient, it is easily seen that
egrad is identically equal to ¼ a, where ¼ and a respectively refer to the interfacial
Proc. R. Soc. Lond. A (2000)
780
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
tension and the speci­ c surface areay (interfacial area per unit volume) between the
phases.z Based upon this identi­ cation, we discern that F grad ² F (a).
In combination, the above results lead to the following functional dependence of
the macroscale free energy density:
F ² F (c; a):
(3.82)
(The dependence of the latter upon the macroscale temperature T is not displayed
explicitly, as in all prior cases.) The identi­ cation of a as a relevant variable in the
macroscale description of two-phase ®ows is particularly signi­ cant. Hassanizadeh &
Gray (1993) also claim that F should depend upon a. However, they o¬er neither a
concrete proof of this fact (as we have done) nor is the basis for F (as required by the
distinction between species and phases) provided in their work. Other researchers
who have undertaken to model directly at the macroscale (with no concession to
microscale considerations) have also intuitively recognized the fact that saturation
alone might not constitute the only variable required for a complete macroscale
description of the phenomena. In a sense, the above analysis constitutes a proof,
albeit one based upon our constitutive assumptions, that the macroscale free energy
F is functionally dependent not only upon saturation, but upon speci­ c surface area
as well.
Upon using (3.68) we further obtain
·
i
² · i (c; a)
(i = 1; 2);
(3.83)
the signi­ cance of which is analysed in the next section.
(k ) Capillary pressure{saturation relationship
In this section we use results obtained in x 3 i, concerning the de­ nition of the
capillary pressure, jointly with the results obtained in the previous section, relating
to the functional dependencies of the macroscale free-energy density and chemical
potentials, to draw pertinent conclusions regarding the existence of a relationship
between capillary pressure and saturation, as well as to the likely source of the
hysteresis repeatedly observed in such experiments.
Upon combining equations (3.76) and (3.77) we obtain
Pc = ·
1
¡
·
2:
(3.84)
This relation is responsible for the phase-equilibrium type of relationship existing
between Pc and the saturation c, since both · 1 and · 2 are functions of c. In our
di¬use Darcyscale model, where the total pressure p and the macroscopic saturation c constitute the appropriate two-phase entities, the existence of a relationship
between the macroscale pressures p1 ; p2 and c is not surprising. In contrast, however,
y The concept of interfacial area is absent when the di¬use interface microscale model is adopted.
Since ¼ and a always appears as the product ¼ a, a microscale researcher can calculate this quantity as
1
½
o
d3 r
t
½ f
1
2
K (rc) 2
using his microscale solutions.
z Incorporation of wetting e¬ects would necessarily involve the solid{®uid surface tensions and areas,
neither of which is present in our expression since we eschewed such e¬ects in our analysis.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
781
other (singular) Darcyscale studies treat the microscale pressures p1 and p2 as being
independent microscale variables. As such, their macroscale counterparts, namely p1
and p2 , necessarily manifest themselves as being independent quantities. As a result,
the existence of a functional relationship between p1 and p2 , as embodied in the
capillary pressure concept, appears to be not only ad hoc, but indeed mystifying.
In combination, equations (3.84) and (3.83), yield Pc ² Pc (a; c). Thus, according to our analysis, the phenomenon of capillary pressure hysteresis arises from
the non-uniqueness of Pc (a; c) under speci¯cation of only one of the two governing variables|namely, the saturation. Speci­ cation of c alone does not su¯ ce to
render the value of Pc unique. Since the variables a and c are wholly independent
of one other, each being ­ xed appropriately by the solution of microscale problem,
situations can arise wherein identical macroscale saturations involving di¬erent as
become possible.y Such circumstances would necessarily lead to di¬erent values of
Pc , thereby explaining the apparent hysteresis loop as observed in Pc versus c curves.
In sum, our hypothesis is that this hysteresis is simply a result of the projection of
the three-dimensional Pc {c{a surface onto the two-dimensional Pc {c plane!
(l ) Dynamical equations for macroscale description
In the previous section we rationalized the nature of the hysteresis observed in
experimental measurements of capillary pressure versus saturation experiments as
being a direct consequence of the fact that the macroscale free energy density, and
hence the capillary pressure is, in addition to the saturation, also a function of the
speci­ c surface area. In this section we complement the above result by expounding
the argument that the mathematical modelling of multiphase ®ows needs to address
the issue of macroscale energy transport in order to complete the set of equations
governing two-phase ®ow phenomena.
As observed in the preceding section, a complete macroscale equation of state
would necessarily involve the relationship of capillary pressure to both speci­ c surface area and saturation. Such an equation of state complements Darcy’s laws (3.73),
(3.74), and the di¬usion equations (3.53). However, these equations do not describe
the evolution of the speci­ c surface area a, the latter constituting a variable wholly
independent of the saturation. Consequently, a complete the set of equations purporting to describe the temporal evolution of the independent variables requires that
we need to consider the macroscale equation governing energy transport.
Based upon the equivalence between the macroscale free energy density and the
macroscale internal energy density, and using the results of x 3 j we ­ nd that
e(c; a) = eb
u lk (c)
+ ¼ a:
(3.85)
Under the assumption of slow, isothermal ®ows with weak gradients, such has pervaded the entirety of our discussion, the macroscale conservation equation for e,
namely (3.56), acquires the following form:
@(" eb u
@(" ¼ a)
+ r (" v¼ a) = ¡
@t
@t
lk )
2
J i gi ¡
+
F^ v;
(3.86)
i= 1
y One possible mechanism by which such a situation might occur would correspond to the scenario
wherein the wetting properties of the ®uids are taken into account.
Proc. R. Soc. Lond. A (2000)
782
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
where J i and F^ obey the uncoupled version of the constitutive equations elucidated
in x 3 f. The above equation serves as the governing equation for the evolution of the
speci­ c surface area a. This evolution equation is non-conservative owing to the fact
that interfacial area can be either created or destroyed by the respective breakup or
coalescence of the droplets comprising the discontinuous phase.
The above equation complements the phase-speci­ c Darcy’s law, equations (3.73){
(3.74), as well as the di¬usion equations, (3.53). Along with a capillary pressure
equation of state of the functional form Pc (a; c), the resulting system of equations
provides a complete description of slow, isothermaly biphasic ° ows through porous
media.
(m) Illustrative example
This ­ nal section explicitly illustrates how our microscale model can be used to
rigorously derive and hence formally justify linear macroscale constitutive relations
that accord with those previously postulated axiomatically on the basis of irreversible
thermodynamic arguments. Additionally, it is shown explicitly how one may, in principle, calculate the Darcyscale phenomenological coe¯ cients appearing therein from
the speci­ ed microscale data. Section 3 f illustrated the form of the constitutive
equations, without however specifying the manner in which the phenomenological
coe¯ cients ¤0 could be theoretically calculated. The analysis of the present section,
designed to rectify this omission, is directed to the elementary case of two-phase ®ow
through a spatially periodic model of a porous medium under the additional constraints of small capillary and Reynolds numbers. At the macroscale this description
eventually leads to a homogeneous multiphase ®ow con­ guration, wherein macroscale
quantities like saturation, velocity, and pressure gradient are all uniform throughout
the porous medium (Whitaker 1986). Complete analytical treatment of the more general case presents insuperable di¯ culties owing to the severe nonlinearities present
in the problem formulation, most notably in ¿c and F (c). However, at least in principle, numerical evaluation of the requisite macroscale phenomenological coe¯ cients
is made possible by their rigorous de­ nitions in terms of the prescribed microscale
data. In turn, this enables use of solutions of the microscale equations to con­ rm the
predictions of linear irreversible thermodynamics, and concomitantly to determine
the numerical values and parametric dependence of the phenomenological coe¯ cients
embodied in the various ¤0 s.
(i) Physical description
In this section we elucidate the basic physics underlying the illustrative example
accompanying the mathematical details elaborated in subsequent sections.
Consider the generic case of uniform two-phase ®ow through a porous medium
animated by an external pressure gradient. This applied macroscopic pressure gradient generates a ®ow at the microscale, resulting in motion of the two species. In the
following sections we analyse the governing microscale equations arising from this
y For non-isothermal ®ows the above equation needs to be replaced by two species-speci­ c energy
equations together with the requirement of local thermal equilibrium in order to satisfy the mandated
equivalence between the required number of equations and the number of independent variables. However,
such issues take us far from the scope of the present work. As such, they will not be addressed here.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
783
pressure gradient. An analytical solution of this problem would be extremely complicated owing to the underlying nonlinearities. These, in turn, are related to the
classical nonlinearities appearing in the singular microscale approach, which result
from having to satisfy interfacial boundary conditions on a surface whose unknown
shape needs to be determined as a part of the solution of the problem.
However, in view of our prior arguments based on irreversible thermodynamic
concepts, it is to be expected that the two-phase Darcy’s law will be applicable in
some linear limit of the problem. Anticipating such a simpli­ cation, in what follows
we analyse the microscale problem in the linear limit. This limit corresponds to
the leading-order, small perturbation approach devised by Taylor (1934), wherein
the interfacial boundary conditions are initially enforced on a surface whose shape
is assumed known a priori. Departures from that assumed shape are then derived
sequentially by employing the solutions of the governing equations (sans the normal
stress interfacial boundary condition) at successively higher orders of approximation
in the perturbation parameter.
We consider the microscale equations and boundary conditions outlined in x 3 c
together with the constitutive equations outlined in x 2 c, subject to the constraint
that the macroscale pressure gradient be equal to the applied pressure gradient.
Equation (3.43) suggests that this macroscale pressure gradient constraint is equivalent to a comparable constraint imposed on the microscale pressure gradient. The
applied macroscopic pressure gradient thereby acts as an animating force driving
the microscale ®ow, which can, in principle, be determined by solving the governing microscale equations. The microscale velocity ­ eld generated in this manner is
manifested as a uniform macroscale ®ow (at least for the case of homogeneous ®ows,
under the assumed spatial and temporal ®ow characteristics outlined previously)|
the dependence of which on the macroscale pressure gradient is subsequently proved
to be linear, in accordance with the Darcy’s laws. This section thereby constitutes
an ab initio proof of the validity of the heretofore empirical two-phase Darcy’s laws,
albeit for the limited case wherein only the linear limit of the microscale equations
is considered.
(ii) Microscale problem
The microscale equations appropriate to a quasi-steady, zero Reynolds number,
small capillary number description of multiphase ®ow is obtained by setting Ca = ¯
and ½ D U=L = ¯ (¯ ½ 1). To simplify the analysis we ignore the e¬ects of gravity,
thereby obtaining the trio of equations:
@»
+ r (v» ) = 0;
@t
¡ ¯ rp + ¯ r ¿v + ¯ r [(rc)(rc) ¡
¯
2
@f
@c
¡
+ r (vc) = ¤r2
@c
@t
(3.87)
jrcj2 I] = 0;
¯ 2 r2 c :
(3.88)
(3.89)
These equations are to be solved subject to the periodicity conditions (3.6) together
with the no-slip and no-®ux boundary conditions (3.15) and (3.16), respectively,
on the bed particle surfaces. Furthermore, subject to a posteriori veri­ cation the
following expansions are proposed for the relevant physical variables (corresponding
Proc. R. Soc. Lond. A (2000)
784
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
expansions for the inner variables being obtained by imposing a tilde ornation):
p=¯
¡1
p¡1 (q1 ; q2 ; n; t) + p0 (q1 ; q2 ; n; t) + O(¯ );
fv; cg = fv; cg0 (q1 ; q2 ; n; t) + ¯ fv; cg1 (q1 ; q2 ; n; t) + O(¯ 2 ):
(3.90)
(3.91)
The ®uid motion is assumed to be driven by an externally imposed, time-independent, macroscale pressure gradient rp, taken to be uniform throughout the porous
medium. This macroscopic pressure gradient, de­ ned through (3.45), serves as a constraint on the solution of the above system of equations. While the solution of the
complete dynamical problem subject to this constraint possesses non-trivial uniqueness properties, the fact that the dynamics in this section is assumed to involve
a quasi-steady response of the pressure ­ eld requires that an instantaneous version of (3.45) (that is, excluding the time averaging) be speci­ ed in order to assure
uniqueness of the resulting solutions of the microscale equations. Thus, for purely
mathematical reasons pertaining to uniqueness, we invoke a quasi-steady analogue
of (3.45), namely that the instantaneous spatial average of the microscale pressure
­ eld constitutes the macroscale con­ guration-speci­ c pressure gradient:
rp =
1
½
o
ds p:
@½
(3.92)
o
(iii) Inner and outer equations
We employ the same notation as in x 2. To leading order in ¯ , a quasi-steady
analysis requires a priori speci­ cation of the location of `interface’ (cf. Nitsche (1989)
and the references cited therein). However, the shape of the interface also constitutes
an unknown of the problem, which has to be determined concurrently with the
solution of the equations of motion. The assumption of small Ca makes it possible for
one to assign, at the leading order, a macroscopic interfacial con­ guration, one whose
curvatures µ1 ; µ2 are independent of (q1 ; q2 ). This requires that a similar condition
be satis­ ed by the normal component of the zeroth-order velocity ­ eld v0 at the
interfacial surface. Continuation of the perturbation expansion to determine higherorder corrections to the interfacial shape is then carried out in a systematic manner
similar to that of x 2, with the concentration ­ eld c used to establish corrections to
the shape (say, by arbitrarily choosing the macroscale interface to coincide with the
isoconcentration surface c = 0:5). Details of such an analysis are similar in spirit to
those outlined in x 2 e, but quickly become algebraically unwieldy at higher order.
Inasmuch as such extended calculations are not wholly relevant in this initial foray
into the ­ eld, they are omitted here in the interests of brevity. Accordingly, we assume
that the leading-order terms in the perturbation expansion su¯ ce to adequately
represent the required microscale ­ elds. The resulting equations are summarized
below.
Outer equations
r v0§ = 0;
rp§
0
Proc. R. Soc. Lond. A (2000)
=·
§
(3.93)
r2 v0§ :
(3.94)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
785
These are to be solved subject to the following jump conditions imposed across the
interface, such conditions being derived in a manner similar to those outlined in x 2 d :
[[v0 ]] = 0;
[[n ¿v0 ]] Is = 0;
(3.95)
(3.96)
ds [[n ¿v0 ]] = 0;
(3.97)
n v0 = u:
(3.98)
I12
Each of the above conditions is to be imposed at the parent surface, the unit normal to
which is denoted by n (not to be confused with the cell number). Further, the integral
condition (3.97) imposed on the stress is not a result of the singular perturbation
analysis, but is rather a result of physical considerations accompanying the small
Ca approximation (Taylor 1934). The singular perturbation analysis results in the
boundary conditions previously derived in x 2 e (wherein equation (3.95) ² (2.59)
and (3.96) ² (2.66) + (2.67)). However, the normal stress condition (2.69) cannot be
satis­ ed on the surface of the drop with an a priori shape assumed. Hence, based on
physical considerations we impose the conditions (3.97) and (3.98) (Nitsche 1989) to
establish the higher-order corrections to the droplet shape.
Inner equations. This section is concerned primarily with the normal component of
the momentum equation in the inner region. Its purpose is to illustrate the fact that
the microscale pressure gradient possesses a purely normal component in the interfacial region, whence its contribution to the macroscopic pressure gradient vanishes
upon integration over the interfacial region contained within the unit cell. Perturbation analysis of both the concentration equation and the tangential component
of the momentum equation are performed in a manner similar to that employed to
establish (2.52), (2.66) and (2.67) in x 2 e.
For the normal component, at the leading order it is found that
¡
@~
c0
@ p~¡1
+ (µ1 + µ2 )
@n
~
@~
n
2
= 0;
(3.99)
representing the Laplace condition imposed on the leading-order outer pressure ­ elds,
p§
¡1 (Nitsche 1989). At the next order in ¯ we ­ nd that
¡
@~
v31
@
@ p~0
2·
+
@n
~
@n
~
@~
n
¡
c1
@~
c0 @~
2 @(· À )
+ 2(µ1 + µ2 )
~
@n
~ @n
~
3 @n
= 0;
(3.100)
where c~1 is the O(¯ ) inner concentration term. Integration yields
¡ p~0 + 2·
@~
v31
@n
~
¡
2
·
3
À +2
(µ1 + µ2 )
c1
@~
c0 @~
~
@~
n @n
d~
n = ¬ (q1 ; q2 );
(3.101)
where, as already noted,
À = h1 h2
@
@q1
v~10
h2
Proc. R. Soc. Lond. A (2000)
+
@
@q2
v~20
h1
+
1
@
v31
1 @~
+ v~30
~ 2 ~h1
@~
n h
~
h2 h1 @ n
:
(3.102)
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
786
Upon taking the limit of the above expression as n
~ ! §1, and using (2.54), (2.55)
and (2.65), it is found that
lim À ² h1 h2
n
~ ¡!
1
v~10
h2
@
@q1
+
v~20
h1
@
@q2
+
@
1 @v30
+ v30 (µ1 + µ2 ) ! 0:
@n
~
h2 h1 @n
(3.103)
The assumption that v~30 ; µ1 ; µ2 are independent of (q1 ; q2 ), requires that
h1 h2
@
@q1
v~10
h2
+
v~20
h1
@
@q2
² independent of (q1 ; q2 );
(3.104)
whence
À ² independent of (q1 ; q2 ):
(3.105)
These intermediate results are required in the next section.
(iv) Composite solutions
Solutions of the inner and the outer equations may be used to obtain a composite
solution, de­ ned generically as (Nayfeh 1973)
Ãc0 = Ã0 + Ã~0 ¡
Ãm c ;
(3.106)
where à m c denotes the matching condition. However, in view of the fact that the
boundary layer is of the intermediate layer kind, it is necessary to de­ ne two (rather
than one) composite solutions:
§
Ãc0
= Ã0§ + Ã~0 ¡
Ãm§c :
(3.107)
Note that these composite ­ elds are continuous across the `interface’. Thus, at the
leading order we take the microscale ­ elds in equations (3.21){(3.29) to be repre§
sented by Ãc0
. Further, we de­ ne the continuous composite ­ eld
¡
+
Ãc0 = Ãc0
+ Ãc0
;
(3.108)
¡
+
with Ãc0
; Ãc0
respectively de­ ned on either side of the interface, as appropriate, and
set equal to zero on the other side.
Use of the above de­ nitions together with the results (2.57) and (2.47) yields
§
vc0
= v0§ ;
(3.109)
vc0 = v0+ + v0¡ :
(3.110)
c§
~0 ;
c0 = c
(3.111)
whence
We also have that
»
c0
=»
0
1
+ (»
0
1
¡
»
0
c0 :
2 )~
(3.112)
d~
n + ­ (~
n) + ± (~
n);
(3.113)
Moreover,
§1
§
p§
c0 = p0 ¡ 2
(µ1 + µ2 )
n
~
Proc. R. Soc. Lond. A (2000)
c1
@~
c0 @~
~
@~
n @n
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
787
where (3.105) has been used to denote
@ v~31
@n
~
2·
¡
+
@v30
@n
+
2·
2
·
3
¡
q30
À ²­:
(3.114)
The ­ elds p~¡1 (p§
~ , are denoted by ± (~
n)
¡1 are constants), which are functions only of n
in the above. The development of the prior section is to be construed as preliminary
to the observation that ­ is functionally dependent only upon n
~ . In addition, we
have that
pc0 = p+c0 + p¡
c0 :
(3.115)
(v) Macroscopic quantities
With equations (3.110){(3.115) de­ ning the continuous microscale ­ elds, equation
(3.26) furnishes the macroscopic velocity,
"» v =
1
½
d3 r [»
o
t
½
0
1
¡
(»
0
1
¡
» 02 )~
c0 ] vc0 :
(3.116)
f
The macroscopic pressure gradient is obtained as
1
rp =
½
o
@½
ds p+1 +
o
1
½
o
@½
ds p¡
1:
(3.117)
o
(The extra terms in (3.113) result in a contribution to (3.92) of the form ¹ (~
n)i3 in the
interfacial region. Integration of this purely normal component over the interfacial
transition region results in no net contribution.)
The macroscopic di¬usion ®ux J i is obtained from (3.39) as
J1 =
½
1
o
r ds [» 01 c~(vc0 ¡
t
@½
v)];
(3.118)
o
where we have used the fact that J~1 = 0 at leading order; J 2 is obtained from J 1 as
J 2 = ¡ J 1:
(3.119)
Computations of the macroscale ­ elds are seen to involve only the solutions of the
outer equations, jointly with the solution for c~0 . The solution for c~0 is completely
determined once the location of the interface is speci­ ed via an appropriate constitutive equation for F (cf. (2.13) and (2.52)). This solution is independent of the
externally applied macroscopic pressure gradient, rp. However, the time-dependence
of the microscale quantities in the above equations arises from the dynamical evolution of the interface, which depends nonlinearly on rp through the relation
dq30
= u(rp):
dt
(3.120)
In the general case the latter equation implies that the macroscale quantities depend
nonlinearly upon rp. However, for the special dynamical evolution considered, namely, time periodic or almost time periodic motions, it is possible to replace the time
Proc. R. Soc. Lond. A (2000)
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
788
(a)
(b)
Figure 5. Space-averaging domain L for (a) periodic motions (the domain L is a straight line);
and (b) almost-periodic motions (due to the densely ¯lling ergodic motion, the domain L constitutes the entire cell).
averaging by a space-averaging procedure (Adler et al . 1985). Consequently, we will
henceforth use the equivalence (refer to ­ gure 5 for the de­ nition of L )
²
:
(3.121)
L
t
As such, it is unnecessary to be concerned explicitly with the relationship between
the dynamical evolution of the system and the macroscale pressure gradient rp.
(vi) Linearity and Darcy’s laws
The outer equations are linear, enabling solutions thereof to be obtained in terms
of modi­ ed tensor ­ elds (V , P , T ) de­ ned as follows:
§
fv § ; p§ ; ¿ § g = fV § ; P § ; T
g rp:
(3.122)
In terms of these modi­ ed ­ elds the outer equations adopt the respective forms
r V0§ = 0;
rP1§
=·
(3.123)
§
r
2
V0§ ;
(3.124)
subject to the following jump conditions imposed across the interface:
I=
[[V0 ]] = 0;
Is [[n T v0 ]] = 0;
(3.125)
(3.126)
ds [[n T
(3.127)
1
½
o
@½
v0 ]]
ds P0+ +
o
= 0;
1
½
o
@½
ds P0¡ :
(3.128)
o
In the terminology of homogenization theories (Bensoussan et al. 1978; Auriault
1987), the above problem constitutes the `unit cell’ problem. Observe that the solution of the above quasi-steady problem is independent of the macroscopic pressure
gradient rp.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
789
With the identi­ cations made in this and previous sections the Darcy multiphase
®ow laws can be obtained as follows. From (3.116),
v = K rp;
(3.129)
in which
K=
1
1
"»
½
L
o
d3 r [»
½
0
1
¡
(»
0
1
¡
» 02 )~
c0 ] Vc0 ;
(3.130)
f
with Vc0 de­ ned in a manner similar to vc0 . Since Vc0 and c~0 are independent of
rp, the above equation constitutes the macroscopic mixture Darcy law. Moreover,
equation (3.130) gives the prescription for calculating the phenomenological tensor
coe¯ cient ¤III appearing in (3.63) (namely, ¤III = K ¡1 ), which was the objective
of this section.
Furthermore, from (3.118),
J 1 = K 0 rp;
(3.131)
wherein
1
K0 =
L
½
o
r ds [» 01 c~(Vc0 ¡
@½
K)]:
(3.132)
o
Using the above relations together with the de­ nitions of the phase-speci­ c macroscale quantities, one obtains the following phase-speci­ c Darcy’s laws:
v1 =
K+
K0
» 1
rp 1 ;
v2 =
K¡
K0
» 2
In the simpli­ ed spatially periodic model analysed above, r·
whence the capillary pressure adopts the simple form
Pc = const:
rp 2 :
k
(3.133)
= 0 and rc = 0,
(3.134)
throughout the system. This spatial uniformity is consistent with the homogeneity
of the macroscale ­ elds.
(vii) Discussion
Calculations presented in this section have served to illustrate the formalism
involved in computing the macroscale ­ elds. The desire for analytical tractability
forced us to choose an extremely simple example. Nevertheless, the analytical scheme
leading to the multiphase Darcy law description is clearly non-trivial despite the
apparent simplicity of the model. Numerical solutions of the unit cell problem also
promise to be formidable. It needs to be emphasized here that despite the resemblance
of the system of outer equations to the conventional singular interface equations, our
formalism de­ nes the pertinent macroscale ­ elds in terms of the complete, continuous ­ elds, rather than in terms of the outer solutions alone. As such, it is likely that
a schism will develop in more comprehensive studies between conventional singular
interface analyses, where one deals with discontinuous ­ elds, and our di¬use interface
analysis, involving continuous ­ elds. Moreover, in our simplistic case, the interfacial
Proc. R. Soc. Lond. A (2000)
790
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
conditions played a subsidiary role in determining the solutions of the unit cell problem. In particular, it was unnecessary to deal with averaging these conditions at any
stage in the analysis (in contrast with the treatments of Whitaker (1986) and Gray
& Hassanizadeh (1990)). A useful exercise would be to illustrate the above formalism
for other similar problems so as to highlight the conceptual advantages and the formalism embodied in our approach. However, such an e¬ort entails manifold details,
a number of which prove to be repetitive in terms of the examples addressed in this
section and in x 2 e. As such, we defer such an exercise to a future publication.
Outstanding issues regarding uniqueness of solutions of the outer equations are
easily resolved by pursuing conventional uniqueness arguments, e.g. employing the
energy dissipation theorem (Keller et al. 1967). As such, further discussion of the
uniqueness issue is omitted here. Questions pertaining to the symmetry properties of
the phenomenological coe¯ cients K; K 0 are also omitted in the interests of brevity,
except for the observation that these dyadics are not generally symmetric unless u
in (3.98) is set to zero. This conclusion is consistent with that of other researchers
(Whitaker 1986; Auriault 1987).
4. Conclusions and summary
The preceding analysis outlines a novel general framework for modelling multiphase
®uid ®ow through porous media. Our perspective involves treating pore-scale ®ow
phenomena from a di® use interface point of view, in contrast to the more usual singular interface viewpoint. In our scheme, two-phase ®ows are likened to the Navier{
Stokes ®ow of a single-phase, two-component mixture possessing steep gradients in
physicochemical properties across the interfacial transition zone. Initially, to keep
the arguments as simple as possible, we examined the admissibility of immiscible
two-phase mixtures under the di¬use interface model in the absence of the porous
medium. This involved a rigorous singular perturbation analysis to derive the conventional bulk-phase equations, as well as the kinematical and dynamical matching
conditions|the latter yielding the usual Laplace interfacial tension boundary condition.
Using such a framework we shifted the focus to comparable ®ows in porous media,
where we de­ ned all macroscale (di¬use Darcyscale) quantities in a rigorous and
physically plausible manner in terms of their microscale, Navier{Stokes counterparts.
Special emphasis was placed on oft-neglected considerations regarding the time evolution of the pore-level ­ elds (Whitaker 1986). Though the analysis was performed
for a special geometric model of porous medium, namely one possessing a spatially
periodic skeletal structure, we also indicated the existence of other cases to which we
believe our scheme could be applied. Using physically based de­ nitions we derived
the di¬use Darcyscale transport equations from their microscale counterparts. Arguments from irreversible thermodynamics were also invoked to suggest linear forms
for the macroscale constitutive equations connecting the dynamical and kinematical
macroscale ­ elds.
We believe that one of the main achievements of this work lies in clarifying the
distinction between phase-speci­ c and species-speci­ c quantities. This identi­ cation
was based upon a combination of physical and thermodynamic arguments. Such an
interrelationship was shown to result from correlating the respective di¬use and singular Darcyscale approaches. This led us phenomenologically to the underlying basis
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
791
for the existence of classical phase-speci­ c Darcy’s laws. In addition, it furnished a
number of insights into both the de­ nition and interpretation of capillary pressure
as a continuum Darcyscale ­ eld concept. Eventually, we used a simple illustrative
example, derived formally from our microscale equations by singular perturbation
methods, to provide some justi­ cation for the linear macroscopic constitutive equations originally proposed purely phenomenologically on the basis of irreversible thermodynamics.
Our development invoked a number of simplifying assumptions at various stages of
the analysis. It is useful to recall these, and hence the potential limitations thereby
imposed, so as to concurrently point towards possible future research directions
involving more realistic models of multiphase ®ows, albeit built upon the generic
framework embodied in our di¬use interface approach. One of our main assumptions
involved the neglect of wetting properties. A wide variety of experimental observations (Stokes et al. 1986) have demonstrated the relevance of wettability and contact
line phenomena to multiphase ®ow in porous media, neither of which phenomena fall
under the purview of this article. Furthermore, we have restricted our attention
exclusively to isothermal ®ows. In non-isothermal situations the averaging procedure as well as the de­ nition of macroscale temperature can be expected to raise
fundamental questions (although we believe these to be relatively unimportant in
most applications). Moreover, no numerical results were cited to corroborate our
macroscale theory. Model porous systems which facilitate numerical computations
of the pertinent microscale problems, simultaneously suggesting comparable experimental observations, represent potential future research directions towards ratifying
the theoretical framework presented here.
Our analysis deliberately straddled the boundary between physical concepts and
mathematical rigor, albeit without straining too much into the latter, as is usual in
such cases. We believe that the di¬use interface, mixture-theory model outlined in
this paper has the potential to clarify a number of fundamental issues prerequisite to
a deeper understanding of two-phase ®ow phenomena in porous media. At a practical
level the model outlined herein can also be used to e¬ect explicit numerical microcontinuum computations based upon representative geometric models of porous media.
Such calculations are greatly facilitated by the fact that the major computational
issue of temporally evolving `interfaces’, together with their concomitant spatial discontinuities, is no longer explicitly present in our di¬use interface model (although
it is, of course, implicitly present). Furthermore, our model can presumably be used
to acquire useful insights into the forms of the pertinent constitutive equations as
well as into the magnitudes of the phenomenological coe¯ cients appearing therein.
This work was supported by the O± ce of Basic Energy Sciences of the US Department of Energy
as well as by a `Grand Challenge Award’ from the Mathematical, Information and Computational
Sciences Division within the O± ce of Computational and Technology Research of the O± ce of
Energy Research of the DOE. V.G. acknowledges a number of useful conversations with Dr
Ehud Gavze and R. Radhakrishnan during the course of this work.
Appendix A. Notation
(Numbers in parentheses indicate the equations in which the symbol is de­ ned or
makes its ­ rst appearance.)
Proc. R. Soc. Lond. A (2000)
792
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
Nomenclature
a
Ca
c
c
D
e
e
eb u lk
F
F
Fb
u lk
F grad
Speci­ c surface area of the interface between the phases (3.81)
Microscale capillary number (2.29)
Microscale composition variable (volume fraction of species 1) (2.5)
Macroscale composition variable of species 1 ( ² s1 ) (3.24)
Deformation gradient tensor (2.20)
Microscale internal energy density (3.14)
Macroscale internal energy density (3.29)
Macroscale internal energy density of a miscible mixture of
composition c (3.85)
Cahn{Hilliard free energy at the microscale (2.13)
Macroscale free energy (3.31)
Macroscale free energy density of a miscible mixture of composition c
(3.80)
The contribution to the macroscale free energy density arising from the
presence of concentration gradients (3.81)
F^
Force manifest at the macroscale due to the stationarity of the
particles (3.55)
Fp
Physicochemical capillary forces at the microscale ² r ¿c (2.18)
F¬ e
External forces acting on the species ¬ at the microscale (2.4)
f (c)
Excess free energy of the mixture at the microscale (2.24)
f 0 (c)
Total free energy of the mixture at the microscale (footnote on p. 745)
fid (c)
The microscale free energy density of an ideal mixture of composition c
(footnote on p. 745)
g1
Gravitational force acting on species 1 (3.63)
g2
Gravitational force acting on species 2 (3.63)
g
Gravitational force per unit mass (3.63)
(h1 ; h2 ; h3 ) Metrical coe¯ cients for the parent surface-­ xed coordinate system
(2.41)
I
Three-dimensional idemfactor
Is
Surface idemfactor (2.2)
(i1 ; i2 ; i3 ) Triad of orthogonal unit vectors characterizing the parent surface-­ xed
coordinate system (2.53)
J¬
Di¬usive ®ux of species ¬ (¬ = 1; 2) at the microscale (2.8)
J¬
Di¬usive ®ux of the ¬ phase (¬ = 1; 2) at the macroscale (3.39)
K
Phenomenological coe¯ cient appearing in the microscale free energy F
(2.13)
K
Permeability tensor (1.1)
k
Boltzmann constant
l
Characteristic dimension of the unit cell of the spatially periodic
porous medium
l
Intrinsic microscale angular momentum at the microscale (3.28)
lk
Basic lattice vector characterizing a unit cell (3.1)
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
M
M
n
n
~
P
Pc
P
rp
rp
rp ¬
(q1 ; q2 ; q3 )
q
q
Re
Rn
R
dR
sp
sk
ds
T
u
v
v¬
v
v¬
¯
"
µi
¤
¤
·
¬
¹
»
»
¬
»
0
¬
»
»
¬
¼
½
o
½
f
½
D
¿v
793
Microscale ®uxes (3.32)
Macroscale ®uxes (3.32)
Normal coordinate in the outer region (2.31)
Normal coordinate in the inner region (2.33)
Microscale stress tensor ² ¡ pI + ¿ (3.41)
Capillary pressure (1.2)
Macroscale stress tensor (3.41)
Microscale pressure gradient
Macroscale pressure gradient (E 7)
Macroscale pressure gradient of the ¬ phase (¬ = 1; 2) (3.69)
Interfacial surface-­ xed coordinate system (2.31)
Microscale heat ®ux vector (3.14)
Macroscale heat ®ux vector (3.47)
Microscale Reynolds number (2.28)
Position vector of a lattice point (3.1)
Macroscale position vector (3.2)
Displacement vector between two adjacent lattice points (3.3)
Surface of the particles within a unit cell (3.15)
Directed surface area of the kth face of the unit cell (3.4)
Macroscale di¬erential surface element ² sk (3.4)
Time period of averaging (3.10)
Translational velocity of a drop (3.98)
Velocity of the mixture (solution) at the microscale (2.7)
Velocity of the species ¬ (¬ = 1; 2) at the microscale (2.3)
Velocity of the mixture at the macroscale (3.26)
Velocity of the ¬ phase (¬ = 1; 2) at the macroscale (3.66)
Dimensionless microscale interfacial thickness (2.29)
Porosity or void fraction (3.21)
Principal curvature of the interface (i = 1; 2) (B 4)
Mobility coe¯ cient in the microscale di¬usion equation (2.11)
Phenomenological coe¯ cients appearing in the generalized force{®ux
relationships and in the derivation of the Darcy law (3.59)
Macroscale chemical potential of the ¬ phase (¬ = 1; 2) (3.68)
Microscale interfacial thickness (2.14)
Density of the mixture at the microscale (2.5)
Mass per unit volume of the species ¬ (¬ = 1; 2) at the microscale (2.3)
Density of pure component ¬ (2.5)
Density of the mixture at the macroscale (3.22)
Density of the ¬ phase (¬ = 1; 2) at the macroscale (3.23)
Interfacial tension
Super­ cial volume of a unit cell (3.21)
Fluid domain within a unit cell (3.20)
Characteristic di¬usion time-scale at the microscale (2.29)
Microscale viscous stress tensor (2.19)
Proc. R. Soc. Lond. A (2000)
794
¿c
¿
¿
Ã
Ã
Á~i
Áci
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
Microscale capillary stress tensor (2.4)
Microscale deviatoric stress tensor ² ¿v + ¿c (3.13)
Macroscale deviatoric stress tensor (3.44)
Generic microscale volumetric density of physical properties (3.5)
Generic macroscale physical properties (3.20)
Physical property volumetric density in the inner region at the ith
order of the perturbation expansion (2.34)
Composite ­ eld (inner + outer - matching conditions) at the ith order
of the perturbation expansion ² Ác+ + Ác¡ (3.106)
Ái§
Outer ­ eld on either side of the interface at the ith order of
perturbation expansion (2.32)
§
Áci
Composite solution ­ eld on either side of the interface at the ith order
of perturbation expansion (3.107)
Jump operator across the cell faces(3.5)
Jump operator across the interface (2.2)
Macroscopic gradient operator (3.49)
Macroscale time-derivative (3.50)
Transposition operator
k:::k
[[: : : ]]
r
@=@ t
y
Appendix B. Scaling in the inner region
The following identities are easily obtained (Happel & Brenner (1983); the truncated versions appearing below retain only those terms necessary in the subsequent
analysis):
~ ² i1 @ + i2 @ + i3 @ ;
(B 1)
r
~ 2 @q2
~
n
¯ @~
h
h1 @q1
~
~
¿ ~3
~²~
~2 @ ¿ 1 + @ ¿ 2 + 1 @
~ Á
;
(B 2)
h1 h
r
~2
~ ~
¯ @n
@q2 ~
@q1 ~
h1 h
h1
h2
1 @¿ ~
1 @
h2 @ ¿ ~
@ ~
h1 @ ¿ ~
@ ~
~1~
~ 2¿ ~ ² h
:
(B 3)
+ 2
+
r
h2
~1~
~
~ h
¯ @n
@q2 ~
@q1 ~
h2 @ n
h1 @q2
h2 @q1
Further, as a consequence of our requirement that the curvatures be macroscopic,
the principal curvatures (Edwards et al. 1991)
h1
1 @~
;
(B 4)
µ1 = ¡
~ n~ = 0
¯ h10 @ n
h2
1 @~
;
(B 5)
µ2 = ¡
~ n~ = 0
¯ h20 @ n
are assumed to be O(1). (Refer to x 2 e for the de­ nitions of h1 ; h2 .)
The following expressions furnish the leading-order terms for some of the quantities
appearing in equations (2.28). Using (B 2) we obtain
~ 1 ~h2
I:D=h
Proc. R. Soc. Lond. A (2000)
@
@q1
v~1
~2
h
+
@
@q2
v~2
~1
h
+
1 @
~
¯ @n
v~3
~
~1
h2 h
;
(B 6)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
¿v = ¡
795
1 @ v~3
v3
1 @~
¡ 23 · (I : D)
+ · i3 i3 2
n
¯ @~
~
¯ @n
v1
1 @~
v3
~ 1 @~
¡ v~1 µ1 +
+ · (i1 i3 + i3 i1 ) h
~
¯ @n
@q1
v2
1 @~
@~
v3
;
¡ v~2 µ2 +
h2
+ · (i2 i3 + i3 i2 ) ~
~
¯ @n
@q2
2
·
3
[i1 i1 + i2 i2 ]
¿c = ¡ Ca¡1
co
(i1 i1 + i2 i2 ) @~
@~
n
¯
(B 7)
2
:
(B 8)
~ ¿~v is
At O(1=¯ 2 ) the contribution to r
i1
@~
v10
@
·
@n
~
@n
~
+ i2
@~
v20
@
·
@n
~
@n
~
+ i3
@
@n
~
@~
v30
;
@n
~
4
·
3
(B 9)
whereas the corresponding contribution at O(1=¯ ) is
@v30 1
@n
~ h20
@~
v10
@~
v10
@~
v30
@~
v11
µ1
(µ1 + µ2 ) + ·
+·
¡ v~10 µ1
+ h10
@n
~
@~
n
@q1
@n
~
@v30 1
@
¡ 23 ·
h10 h20 h30
@n
~ h10
@q2
@~
v20
@~
v20
@~
v30
@~
v21
@
µ2
(µ1 + µ2 ) + ·
¡ v~20 µ2
+·
+ h10
·
+
@n
~
@~
n
@q1
@n
~
@n
~
· @v20
@
· @v10
@
+ h10 h20 h30
h10 h20 h30
~
@q2 h20 @ n
~
@q1 h20 @ n
2 @(· À ) 2 @ v~30
@~
v31
@
(µ1 + µ2 ) :
(B 10)
+ 3·
¡
2·
+
@n
~
~
3 @n
@n
~
@n
~
@
¡
@q1
@
·
+
@n
~
i1 h10 h20 h30
+ i2
+ i3
2
·
3
In the above, À denotes the quantity
À = h10 h20
@
@q1
v~10
h20
+
@
@q2
v~20
h10
+
1
@
v31 )
1 @(~
+ v~30
~1
@n
~ ~
~
h20 h10 @ n
h2 h
:
(B 11)
Using the solution of the di¬usion equation (2.52) at leading order, we ­ nd that
~ ¿~c is identically zero, whereas at O(1=¯ ) we obtain
the O(1=¯ 2 ) contribution to r
c0
~ ¿~c = i3 Ca¡1 (µ1 + µ2 ) @~
r
@~
n
2
:
(B 12)
Summarized below are the respective forms of the various matching conditions
derived within the framework of the singular surface model for the speci­ c geometry chosen to represent the interface (Edwards et al . 1991). The velocity matching
condition for this material interface is
v0+ = v0¡
Proc. R. Soc. Lond. A (2000)
at
n = 0;
(B 13)
796
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
or, equivalently,
[[v0 ]] = 0:
(B 14)
The stress matching condition is
n [[P0 ]] = ¡ 2¼ Hn
at
n = 0;
(B 15)
where ¼ denotes the interfacial tension, n the unit normal to the interface (pointing
in the direction from ¡ to +), and
H ² 12 (µ1 + µ2 )
(B 16)
the mean curvature of the interface. In terms of the semi-orthogonal curvilinear
coordinates (q1 ; q2 ; n) the above conditions can be written as
@v30
@n
¡ p0 + 2·
@v
@v30
+ 10 ¡ v10 µ1
@n
@q1
@v20
@v30
¡ v20 µ2
+
h20
@n
@q2
·
h10
·
= ¡ 2¼ H;
(B 17)
= 0;
(B 18)
= 0:
(B 19)
Consistent with our earlier assumptions, surface rheological e¬ects have been assumed absent.
Appendix C. Proof of theorem 3.3
This section is largely based on Brenner (1972). We use the following identities,
proofs of which are left as an exercise for the interested reader.
Identity C 1. A time periodic (or almost time periodic) function à possessing a
spatially periodic gradient can be decomposed as follows:
0
¶
à = Ã(R;
t) + R F + R F (t);
(C 1)
¶
where Ã(R;
t) is a spatially and time-periodic function, and
F =
1
½
o
t
0
@½
ds Ã;
(C 2)
F:
(C 3)
o
ds à ¡
F =
@½
o
It follows easily that
0
F = 0:
(C 4)
t
Identity C 2. For ö a spatially periodic function, the following identity holds:
t
Proc. R. Soc. Lond. A (2000)
1
ds ö = sk
½ o
sk
¶
r ds Ã:
t
@½
o
(C 5)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
797
In what follows, we outline the proof of theorem 3.3.
Proof . The macroscopic ®ux density tensor M is de­ ned as
ds M :
sk M =
t
(C 6)
sk
Using identity C 1 we have that
0
¶ (R; t) + R G + R G (t)];
ds [M
sk M =
t
(C 7)
sk
0
where G and G are de­ ned as in (C 2) and (C 3). Upon using (C 4) and neglecting
O(klk) terms (i.e. writing R = R + O(klk)], we obtain
¶ (R; t) + sk R G:
ds M
sk M =
t
(C 8)
sk
Identity C 2 thereby yields
sk M =
½
1
o
¶ + sk R G:
r ds M
sk
t
@½
(C 9)
o
¶ in terms of M , G and G0 (t), and upon further
Using identity C 1 to write M
neglecting terms of order O(klk), one thereby obtains
sk M =
1
½
o
r ds M +
sk
t
@½
o
½
1
o
0
r ds G :
sk
t
@½
(C 10)
o
Subsequent use of (C 4) jointly with the fact that sk can be of arbitrary magnitude
and direction thereby furnishes a proof of the ­ rst part of theorem 3.3.
To prove the second part of the theorem we observe that the spatial periodicity of
the porous medium requires that sk fng = sk fn + dng, where n + dn denotes the
face of an adjacent cell. Consequently, (3.32) implies that
sk [M (Rn+
d n)
¡
ds [M (sk fn + dng) ¡
M (Rn )] =
M (sk fng)]:
s k fn g
t
(C 11)
However, since Rn ² R and Rn+
[M (Rn+
d n
² R + dR, it follows that
d n)
¡
M (Rn )] = lk rM :
Application of identity C 1 to the right-hand side of (C 11) yields
ds [M (sk fn + dng) ¡
t
M (sk fng)] = lk G sk :
sk
Upon combining the above two results, one veri­ es the second part of theorem 3.3.
Proc. R. Soc. Lond. A (2000)
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
798
Appendix D. On the equivalence of mass ° ux and momentum
density de¯nitions of the macroscale velocity
Consider (3.36), representing the de­ nition of the macroscopic velocity based upon
the mass ®ux:
1
(D 1)
r ds ½v:
"» v =
½ o t @½ o
Gauss’s divergence theorem yields
1
"» v =
½
r (½vr) d3 r;
o
t
½
(D 2)
o
in which we have used the fact that in our di¬use interface model all ­ elds are
continuous throughout the unit cell, except at the surface of the solid particle, where
v = 0. The right-hand side can be simpli­ ed further as
"» v =
1
½
rr (½v) d3 r +
o
t
½
o
1
½
o
½v d3 r;
t
½
(D 3)
o
which upon using the microscale equation of continuity (2.9) for the mixture can be
written alternatively as
"» v =
1
½
r ¡
o
t
½
o
@»
@t
d3 r +
1
½
o
½v d3 r:
t
½
(D 4)
o
Interchange of the space and time integrations in the ­ rst integral gives
"» v =
1
½
rd3 r
o
½
¡
t
o
@»
@t
+
½
1
o
½v d3 r:
t
½
(D 5)
o
The ­ rst integral vanishes in view of the assumed steadiness of the macroscale ­ elds,
whence (3.26) is obtained.
Appendix E. On the de¯nitions of the macroscopic pressure
and the stress tensor
The macroscopic stress tensor is de­ ned in (3.41). Separation of the latter into respective pressure and deviatoric stress contributions is arbitrary. Usually, the quantity of
signi­ cance is the pressure gradient rather than the pressure itself (Brenner 1972).
Here, we adopt the convention of assigning the mean stress to the pressure term.
Thus, we de­ ne
d ef
p = ¡
1
I
3
: P;
(E 1)
whence
¿ =P ¡
1
I(I
3
: P ):
(E 2)
This decomposition furnishes expressions (3.43) and (3.44). Upon applying identity C 1 to (3.44) one obtains, upon neglecting O(klk) terms,
¿ =¡
Proc. R. Soc. Lond. A (2000)
½
1
o
(rI ¡
t
@½
o
1
Ir)
3
ds (» vv0 + P¶ ):
(E 3)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
799
(The value of the integral multiplying R G is identically zero.) Since the above
quadrature involves averaging spatial- and time-periodic functions respectively over
a unit cell and one time period, we ­ nd that
r ¿ = 0:
(E 4)
r P = ¡ rp:
(E 5)
Thus,
Using theorem 3.3 it follows that
1
rp = ¡
½
o
ds P ;
(E 6)
ds (¡ pI + ¿ ):
(E 7)
t
@½
o
which requires that
1
rp = ¡
½
o
t
@½
o
Explicitly,
rp =
1
½
o
1
ds p ¡
t
@½
o
½
ds ¿ :
o
t
@½
(E 8)
o
The second integral vanishes in view of the spatial periodicity of ¿ , whereupon (3.45)
is recovered.
Appendix F. Generalized Darcy’s law
We assume a general form for the constitutive equations, allowing for coupling
between the ®uxes. However, recognition that J 1 + J 2 = 0 requires that only two of
the phenomenological tensors ¤1k , ¤2k (k = 1; 2) be independent quantities. Consequently,
J 1 = ¡ J 2 = ¤01a r(·
1
¡
·
2)
+ ¤01b v
(F 1)
¡
·
(F 2)
and
F = ¤01c v + ¤01b r(·
1
2 );
where the symmetry of the coupling tensors has been invoked. Substituting (F 2)
into the linear momentum equation (3.54) gives
2
» i g i + ¤01b r(·
¡ rp + ¤01c v + "
1
¡
·
2)
= 0;
(F 3)
i= 1
equivalently,
2
v = (¤01c )¡1
rp ¡
» ig i
"
i= 1
Proc. R. Soc. Lond. A (2000)
¡
(¤01c )¡1 ¤01b r(·
1
¡
·
2 ):
(F 4)
800
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
(When coupling between the ®uxes is neglected, the above reduces to the Darcy
mixture law (3.63).) Furthermore, use of the relations
J 1 = " » 1 (v 1 ¡
v);
J2 = "» 2 (v 2 ¡
v)
(F 5)
yields
2
"» 1 v 1 = ¤00a r(·
1¡
·
00
2 ) + ¤b
rp ¡
"
» ig i
(F 6)
» i gi ;
(F 7)
i= 1
and
2
"» 2 v 2 = ¤00d r(·
2¡
·
00
1 ) + ¤c
rp ¡
"
i= 1
where the ¤00a , ¤00b , ¤00d and ¤00c are combinations of ¤01a , ¤01b and ¤ 01c . Knowledge of
the explicit forms adopted by ¤00a , ¤ 00b , ¤00d ; and ¤00c will prove unnecessary; however,
no symmetry in the coupling tensors could be identi­ ed despite the fact that only
three of the above four phenomenological coe¯ cients are independent. Use of (3.72)
gives
r(·
1
¡
·
2)
=
» r p2
;
» 2
» rp 1
¡
» 1
(F 8)
which may be written alternatively as
r(·
1
¡
·
2)
=
r p1
¡
» 1
"g1
¡
r p2
¡
» 2
"g2 :
(F 9)
Substitution into in (F 6) and (F 7) thereby yields
v 1 = (¤011 )¡1 (rp1 ¡
v2 =
(¤021 )¡1
( rp 1 ¡
"» 1 g 1 ) + (¤012 )¡1 (rp2 ¡
"» 1 g 1 ) +
(¤022 )¡1
( rp 2 ¡
"» 2 g 2 );
(F 10)
"» 2 g 2 );
(F 11)
where, again, the exact forms of the ¤0ij (i; j = 1 : : : 4) are unimportant.
References
Adler, P. M. 1992 Porous media: geometry and transports. Boston: Butterworth-Heinemann.
Adler, P. M. & Brenner, H. 1988 Multiphase ° ow in porous media. A. Rev. Fluid Mech. 20,
35{59.
Adler, P. M., Zuzovsky, M. & Brenner, H. 1985 Spatially-periodic suspensions of convex particles
in linear shear ° ows. 1. Description and kinematics. Int. J. Multiphase Flow. 11, 361{385.
Antanovskii, L. K. 1995 Phase ¯eld model of capillarity. Phys. Fluids A 7, 747{753.
Aris, R. 1962 Vectors, tensors, and the basic equations of ° uid mechanics. New Jersey: PrenticeHall.
Auriault, J. L. 1987 Nonsaturated deformable porous media: quasistatics. Transp. Porous Media.
2, 45{64.
Auriault, J. L. & Sanchez-Palencia, E. 1986 Remarques sur la loi de Darcy pour les ¶ecoulements
biphasiques en milieu poreux. J. de M¶ech. Th¶eorique et Appl. (Special Issue), pp. 141{156.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
801
Bear, J. & Bachmat, Y. 1984 Transport phenomena in porous media|basic equations. In Fundamentals of transport phenomena in porous media (ed. J. Bear & M. Y. Corapciglu), pp. 3{61.
NATO ASI Series E, Applied Science. Boston: Martinus Nijho® .
Bedford, A. & Drumheller, D. S 1983 Recent advances: theories of immiscible and structured
mixtures. Int. J. Engng Sci. 21, 863{960.
Bensoussan, A., Lions, J. L. & Papanicolaou, G. 1978 Asymptotic analysis for periodic structures .
New York: North-Holland.
Bohr, H. 1947 Almost periodic functions. New York: Chelsea.
Brenner, H. 1970 Rheology of two phase systems. A. Rev. Fluid Mech. 2, 137{176.
Brenner, H. 1972 Elements of transport processes in porous media. (Nine folios of unpublished
notes.)
Brenner, H. 1980 Dispersion resulting from ° ow through spatially periodic porous media. Phil.
Trans. R. Soc. Lond. A 297, 81{133.
Brenner, H. & Edwards, D. A. 1993 Macrotransport processes. Boston: Butterworth-Heinemann.
Brenner, H. & Nadim, A. 1996 The Lorentz reciprocal theorem for micropolar ° uids. J. Engng
Math. 30, 169{176.
Brillouin, L. 1946 Wave propogation in periodic structures . New York: McGraw-Hill.
Broadbent, S. R. & Hammersly, J. M. 1957 Percolation processes. 1. Crystals and mazes. Proc.
Camb. Phil. Soc. 53, 629{641.
Buckingham, E. 1907 Studies on the movement of soil mixtures. Bulletin, US Dept of Agriculture,
vol. 38. Washington, DC.
Bu® , F. P. & Saltsburg, H. 1957 Curved ° uid interfaces. 2. The generalized Neumann formula.
J. Chem. Phys. 26, 23{31.
Caginalp, G. 1986 An analysis of a phase ¯eld model of a free boundary. Arch. Ration. Mech.
Analysis 92, 206{245.
Caginalp, G. & Fife, P. C. 1988 Dynamics of layered interfaces arising from phase boundaries.
SIAM J. Appl. Math. 48, 506{518.
Cahn, J. W. 1965 Phase separation by spinodal decomposition in isotropic systems. J. Chem.
Phys. 42, 93{99.
Cahn, J. W. 1977 Critical point wetting. J. Chem. Phys. 66, 3667{3672.
Cahn, J. W. & Hilliard, J. E. 1958 Free energy of a nonuniform system. 1. Interfacial free energy.
J. Chem. Phys. 28, 258{267.
Cahn, J. W. & Hilliard, J. E. 1959 Free energy of a nonuniform system. 3. Nucleation in a
two-component incompressible ° uid. J. Chem. Phys. 31, 688{699.
Dahler, J. S. & Scriven, L. E. 1963 Theory of structured continua. Proc. R. Soc. Lond. A 275,
504{527.
Davis, H. T. 1995 Statistical mechanics of phases, interfaces, and thin ¯lms. New York: VCH.
Davis, H. T. & Scriven, L. E. 1982 Stress and structure in ° uid interfaces. Adv. Chem. Phys.
49, 357{454.
deGennes, P. G. 1980 Dynamics of ° uctuations and spinodal decomposition in polymer blends.
J. Chem. Phys. 72, 4756{4763.
deGennes, P. G. 1983 Theory of slow biphasic ° ows in porous media. PhysicoChem. Hydrodyn.
4, 175{185.
deGroot, S. R. & Mazur, P. 1984 Non-equilibrium thermodynamics. New York: Dover.
Dobran, F. 1984 Constitutive equations for multiphase mixtures of ° uids. Int. J. Multiphase
Flow 10, 273{305.
Dussan, V. E. B. 1979 On the spreading of liquids on solid surfaces: static and dynamic contact
lines. A. Rev. Fluid Mech. 11, 371{400.
Edwards, D. A., Brenner, H. & Wasan, D. T. 1991 Interfacial transport processes and rheology .
Boston: Butterworth-Heinemann.
Proc. R. Soc. Lond. A (2000)
802
Downloaded from http://rspa.royalsocietypublishing.org/
V. Ganesan and H. Brenneron December 3, 2014
Eliassen, J. D. 1963 Interfacial mechanics. PhD thesis, University of Minnesota, Minneapolis.
Fixman, M. 1962 Viscosity of critical mixtures. J. Chem. Phys. 36, 310{318.
Fleming, P. D., Yang, A. J. M. & Gibbs, J. H. 1976 A molecular theory of interfacial phenomena
in multicomponent systems. J. Chem. Phys. 65, 7{17.
Gibbs, J. W. 1957 The collected works of J. Willard Gibbs. New Haven: Yale University Press.
Goodrich, F. C. 1981 The theory of capillary excess viscosities. Proc. R. Soc. Lond. A 374,
341{370.
Gray, W. G. & Hassanizadeh, S. M. 1990 Mechanics and thermodynamics of multiphase ° ow in
porous media including interphase boundaries. Adv. Water Resour. 13, 169{186.
Gunton, J. D., San Miguel, M. & Sahni, P. S. 1983 The dynamics of ¯rst-order phase transitions.
In Phase transitions and critical phenomena (ed. M. S. Green & C. Domb), vol. 8, pp. 267{466.
New York: Academic.
Happel, J. & Brenner, H. 1983 Low Reynolds number hydrodynamics. The Netherlands: Kluwer.
Hassanizadeh, S. M. & Gray, W. G. 1993 Thermodynamic basis of capillary pressure in porous
media. Water Resour. Res. 29, 3389{3405.
Hohenberg, P. C. & Halperin, B. I. 1977 Theory of dynamic critical phenomena. Rev. Mod.
Phys. 49, 435{480.
Ishii, M. 1975 Thermo-° uid dynamic theory of two-phase ° ows . Paris: Eyrolles.
Kalaydjian, F. 1993 A macroscopic description of multiphase ° ow in porous media involving
spacetime evolution of ° uid/° uid interface. Transp. Porous Media. 2, 537{552.
Keller, J. B. 1985 Darcy’ s law for ° ow in porous media and the two-space method. In Physical
mathematics and nonlinear partial di® erential equations (ed. J. H. Lightbourne & S. M.
Rankin), pp. 429{443. Lecture Notes in Pure and Applied Mathemtics, vol. 30. NATO ASI
Series E, Applied Science. Boston, MA: M. Decker.
Keller, J. B., Rubenfeld, L. A. & Molyneux, J. E. 1967 Extremum principles for slow viscous
° ows with applications to suspensions. J. Fluid Mech. 30, 97{125.
Kirkwood, J. G. 1967 Selected topics in statistical mechanics. New York: Gordon & Breach.
Kirkwood, J. G. & Bearman, R. J. 1958 Statistical mechanics of transport processes. 11. Equations of transport in multicomponent systems. J. Chem. Phys. 28, 136{145.
Kirkwood, J. G. & Bu® , F. P. 1949 The statistical mechanical theory of surface tension. J.
Chem. Phys. 17, 338{343.
Langer, J. S. 1974 Metastable states. Physica 73, 61{72.
Larson, R. G., Scriven, L. E. & Davis, H.T. 1977 Percolation theory of residual phases in porous
media. Nature 268, 409{413.
Leverett, M. C. 1941 Capillary behavior in porous media. Trans. AIME 142, 152{169.
Marle, C. M. 1977 On macroscopic equations governing multiphase ° ow with di® usion and
chemical reactions in porous media. Int. J. Engng Sci. 20, 643{662.
Marle, C. M. 1981 Multiphase ° ow in porous media. Houston: Gulf.
Mavrovouniotis, G. M. & Brenner, H. 1993 A micromechanical investigation of interfacial transport processes. 1. Interfacial conservation equations. Phil. Trans. R. Soc. Lond. A 345, 165{
207.
Mavrovouniotis, G. M., Brenner, H., Edwards, D. A. & Ting, L. 1993 A Micromechanical investigation of interfacial transport processes. 2. Interfacial constitutive equations. Phil. Trans.
R. Soc. Lond. A 345, 209{228.
Mokadam, R. G. 1961 Thermodynamic analysis of the Darcy law. Trans. Am. Soc. Mech. Engrs
83, 208{213.
Nayfeh, A. H. 1973 Perturbation methods. New York: Wiley.
Neumann, S. P. 1977 Theoretical derivation of Darcy’ s law. Acta Mechanica 25, 153{170.
Nigmatulin, R. I. 1979 Spatial averaging in the mechanics of heterogeneous and dispersed systems. Int. J. Multiphase Flow 5, 353{385.
Proc. R. Soc. Lond. A (2000)
Downloaded
from http://rspa.royalsocietypublishing.org/
3, 2014
A di® use interface
model of two-phase ° ow on
inDecember
porous media
803
Nitsche, L. C. 1989 Multiphase ° ows through spatially-periodic models of porous media. PhD
thesis, Massachusetts Institute of Technology, Cambridge, USA.
Nitsche, L. C. & Brenner, H. 1989 Eulerian kinematics of ° ow through spatially-periodic models
of porous media. Arch. Ration. Mech. Analysis 107, 225{292.
Nye, J. F. 1985 Physical properties of crystals: their representation by tensors and matrices.
Oxford: Clarendon.
Ono, S. & Kondo, S. 1960 Molecular theory of surface tension in liquids. In Encyclopedia of
physics (ed. S. Flugge), vol. 10, pp. 134{280. Springer.
Ottino, J. M. 1989 The kinematics of mixing: stretching, chaos and transport . Cambridge University Press.
Pavone, D. 1989 Macroscopic equations derived from space averaging for immiscible two-phase
° ow in porous media. Rev. I Fr. Petrol. 44, 29{41.
Pego, R. L. 1989 Front migration in the nonlinear Cahn{Hilliard equation. Proc. R. Soc. Lond.
A 422, 261{278.
Philip, J. R. 1970 Flow in porous media. A. Rev. Fluid Mech. 2, 177{204.
Quintard, M. & Whitaker, S. 1988 Two-phase ° ow in heterogeneous porous media: the method
of large scale averaging. Transp. Porous Media. 3, 357{413.
Scriven, L. E. 1960 Dynamics of a ° uid interface. Equation of motion for Newtonian surface
° uids. Chem. Engng Sci. 12, 98{108.
Slattery, J. C. 1964 Kinematics of di® usion in a heterogeneous surface. Chem. Engng Sci. 19,
453{455.
Starovitov, V. N. 1994 A model of motion of a two-component liquid with e® ect of capillary
forces. Zh. Prikl. Mekh. Tekhn. Fiz. 35, 86{92.
Stokes, J. P., Weitz, D. A., Gollub, J. P., Dougherty, A., Robbins, M. O., Chaikin, P. M. &
Lindsay, H. M. 1986 Interfacial stability of immiscible displacement in a porous medium.
Phys. Rev. Lett. 57, 1718{1721.
Taylor, G. I. 1934 The formation of emulsions in de¯nable ¯elds of ° ow. Proc. R. Soc. Lond.
A 146, 41{48.
Teltezke, G. F., Scriven, L. E. & Davis, H. T. 1982 Gradient theory of wetting transitions. J.
Colloid Interface Sci. 87, 550{571.
Tolman, R. C. 1949 The super¯cial density of matter at a liquid{vapor boundary. J. Chem.
Phys. 17, 118{127.
Truesdell, C. 1966 Continuum mechanics. 1. The mechanical foundations of elasticity and ° uid
dynamics. New York: Gordon & Breach.
Turski, L. A & Langer, J. S. 1980 Dynamics of a di® use liquid{vapor interface. Phys. Rev. B 22,
2189{2195.
Whitaker, S. 1986 Flow in porous media. 2. The governing equations for immiscible two-phase
° ow. Transp. Porous Media. 1, 105{125.
Whitaker, S. 1999 The method of volume averaging . Boston: Kluwer.
Proc. R. Soc. Lond. A (2000)