HAL Id: hal-03150992
https://hal.archives-ouvertes.fr/hal-03150992
Submitted on 24 Feb 2021
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
A novel one-domain approach for modeling ow in a
uid-porous system including inertia and slip eects
F. Valdés-Parada, Didier Lasseux
To cite this version:
F. Valdés-Parada, Didier Lasseux. A novel one-domain approach for modeling ow in a uid-porous
system including inertia and slip eects. Physics of Fluids, American Institute of Physics, 2021, 33
(2), pp.022106. �10.1063/5.0036812�. �hal-03150992�
A novel one-domain approach for modeling
flow in a fluid-porous system including
inertia and slip effects
Phys. Fluids 33
, 022106 (2021); doi: 10.1063/5.0036812
F. J. Vald
es-Parada
1
and D. Lasseux
2,a)
AFFILIATIONS
1
Departamento de Ingenier
ıa de Procesos e Hidr
aulica, Universidad Aut
onoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186,
09340 Ciudad de M
exico, Mexico
2
I2M, UMR 5295, CNRS, Univ. Bordeaux, 351 Cours de la Lib
eration, 33405 Talence Cedex, France
a)
Author to whom correspondence should be addressed: didier.lasseux@u-bordeaux.fr
ABSTRACT
A
new one-domain approach is developed in this work yielding an operational average description of one-phase flow in the classical Beavers
and Joseph configuration including a porous medium topped by a fluid channel. The model is derived by considering three distinct regions:
the homogeneous part of the porous domain, the inter-region, and the free fluid region. The development is carried out including inertial
flow and slip effects at the solid–fluid interfaces. Applying an averaging procedure to the pore-scale equations, a unified macroscopic
momentum equation, applicable everywhere in the system and having a Darcy form, is derived. The position-dependent apparent permeabil-
ity tensor in this model is predicted from the solution of two coupled closure problems in the inter-region and in the homogeneous part of
the porous medium. The performance of the model is assessed through in silico validations in different flow situations showing excellent
agreement between the average flow fields obtained from direct numerical simulations of the pore-scale equations in the entire system and
the prediction of the one-domain approach. Furthermore, validation with experimental data is also presented for creeping flow under no-slip
conditions. In addition to the fact that the model is general from the point of view of the flow situations it encompasses, it is also simple and
novel, hence providing a practical and interesting alternative to models proposed so far using one- or two-domain approaches.
I. INTRODUCTION
Modeling one-phase flow in a coupled fluid-porous medium sy s-
tem has been the su bject of active research since the pioneering work
of Beavers and Joseph.
1,2
Applications of such a configuration are
numerous, ranging from hydrology to chemical engineering. The
main difficulty lies in a physically relevant way of reconciling a macro-
scopic description of the flow with a Darcy equation (or its variants)
applicable in the bulk of the porous medium to the Navi er–Stokes
equation in the bulk fluid. Indeed, such a description requires one to
account for the rapid variat ion of to pology (and of the flow) in the sep-
aration zone betwe en the two media, as evidenced, for inst ance, in th e
experimental works of Goharzadeh et al.,
3
Morad and Khali li,
4
or
more recently, by Terzis et al.
5
As a generic configuration envisaged by Beavers and Joseph,
1
consider the system sketched in Fig. 1 consisting of a channel partially
filled with a porous medium. Let a single incompressible and
Newtonian fluid (the b -phase) flow through and above the porous
medium under steady conditions. This system may be decomposed
into three distinct regions (see Fig. 1)asfollows.
A. The x-region
This region corresponds to the homogeneous part of the porous
medium. For simplicity in the analysis, the r-phase, representing the
solid phase constitutive of the porous matrix, is assumed to be rigid
and homogeneous. Consequently, in this region, the porosity, the seep-
age velocity, and the perm eabil ity are position-independent quantities.
Here, Darcy’s law, or its modifications, is applicable as a result of
upscaling the microsc ale equations by making use of an averaging
domain representative of the porous structure and physical mecha-
nisms at play. Typically, the averaging volume corresponds to one or
more unit cell(s) if the structure is assimilated to be a periodic one.
B. The g x inter-region
This region is the transition zone where the average (or macro-
scopic) velocity experiences abrupt changes. For creeping-flow prob-
lems, the extent of this region has been proposed to be, at most, on the
order of 10 ‘
c
,
6
‘
c
being the size of a geometrical unit cell representa-
tive of the porous material. Note that, even if the porous medium
structure is assumed to be homogeneous, the velocity, the porosity,
and the permeability exhibit spatial variations in this part of the sys-
tem. Modeling flow in this region has been addressed using a general-
ized transport equati on (or penalization approaches
7,8
)asexplained
below, or even by direct numerical simulations.
9,10
C. The g-region
This region corresponds to the port ion of the free fluid above the
porous medium where the flow is one-dimensional and parallel to
the channel axis. It is limited by th e upper wall (also cons idered as the
r-phase) at z ¼ h. Here, the volume fraction of the fluid phase in the
averaging domain is equal to 1, except near the walls, and flow is
described by the Navier–Stokes equations.
Traditionally, flow and transport in this system have been studied
using a two-domain approac h, in which the average equations in the
g-region are coupled to those in the x-region by me ans of suitable
boundary conditions. Another alternative is the one-domain approach,
in which a single average equation is used to describe transport in the
three reg ions of the system. Both approaches have strengths and limi-
tations, which are briefly discussed in the following paragraphs.
The two-domain approach has rece ived far more attention than
the one-domain approach. This may be expla ined by the fact that the
model invo lves balance equations with constant transport coefficients ,
thus making their solution fairly easy to achieve, which represents an
appealing strength. However, there is no consensus about which equa-
tions should be used in the porous medium and which boundary con-
ditions should be employed. Originally, Beavers and Joseph
1
proposed
to couple Darcy’s law with the Stokes equation by means of a jump
condition in the velocity, which has the st ructure of a Newt on’s
cooling law equation. Late r, Neale and Nader
11
proposed to couple the
Darcy–Brinkman and Stokes equations with conditions of continuity
of both the stress and ve locity. Later, Ochoa-Tapia and Whitaker
12
fol-
lowed a similar path, albeit they proposed to consider a discontinuity
in the viscous stress. The jump conditions proposed by both Beavers
and Joseph
1
and Ochoa-Tapia and Whita ker
12
were written in te rms
of adjustable coefficients that needed to be determined experimentally.
In addition, the position of the dividing surf ace where th ese boundary
conditions are applied cannot be arbitrarily fixed. These issues have
been discussed in many works. Among them is the one by Vald
es-
Parada et al.,
13
in which an iterative methodology was proposed to
determine the dividing surface position as well as the values of the
jump coefficients in these models. These authors also derived a more
general two-domain model for which there are, in general, discontinu-
ities of both the velocity and stress. Recently, the Beavers and Joseph
boundary condition has been applied to ine rtial flow conditions, and it
was found that the value of the adjustable coefficient increases with the
Reynolds number.
14
Nevertheless, this boundary condition has been
recently found to be unsuitable for arbitrary flow directions.
15
Amore
detailed review and discussion of the Beavers and Joseph boundary
condition are available from Nield
16
and Auriault.
17
In the one-domain approach, the system is regarded as a pseudo-
continuum, and a single equation is used to model the flow at the mac-
roscopic scale everywhere. A heur istic or penalization approach is to
postulate that such a model is the Navier–Stokes equation with a
Darcy term involving a position-dependent permeability. A more rig-
orous approach was proposed by Ochoa-Tapia and Whitaker,
12
who
derived a generalized tra nsport equation resulting from averaging and
not upscaling the Stokes equation. This equation is more complicated
than th e penalization appro ach, and it has archival more than practical
value. In fact, it was shown in this reference that the one-domain
approach is necessary in the derivation of the jump conditions using
the volume averaging method,
18
justifying its archival value.
Furthermore, the methodolog y proposed by Vald
es-Parada et al.
13
showed that it is necess ary to accoun t for the spatial variations of the
effective coefficients present in the one -domain approach in order to
determine the values of the jump coefficients and the dividing surface
position. This difficulty in predicting the spatial variations of the
permeability and porosity may explain why the one-domain
approach has been less-frequently used than the two-domain
approach. In addition, the corresponding average model is, in gen-
eral, more complicated than the Darcy–Brinkman equation as orig-
inally derived by Ochoa-Tapia and Whitaker.
12
The few
applications of this modeling alternative use ad hoc coefficient s and
average expressions that are not alwa ys in agreement with the
results obtained from pore-scale simulations (s ee for example Refs.
7, 8,and19), and only an approximate closure scheme has been
proposed to compute the spatial variation s of the permeability.
6
Moreover, existing analyses have been mainly focused on flow in
the creeping regime, in the absence of rarefaction effects .
FIG. 1. Schematic representation of the flow in a channel partially obstructed with a
porous medium, showing the three regions of the system, the volume averaging
domain of radius r
0
, and the characteristic lengths L for the macroscopic size of the
porous medium, ‘
b
in the fluid phase, b, and ‘
r
in the solid phase, r, within the
porous medium. The x, g x, and g regions, respectively, correspond to
z z
x
; z
x
z z
g
, and z
g
z h. The unit normal vector at the solid-
fluid interface, directed from the b-phase to the r-phase, is denoted by n.
The objective of this work is to address these issues by proposing a
closed one-domain approach in the three-region system for one-phase
flow including inertial and Knudsen effects in the slip regime. The
model derived with this approach has a simple structure facilitating its
application, nevertheless different from the modified Darcy–Brinkman
equation used in the literature. The analysis is focused on the Beavers
and Joseph configuration sketched in Fig. 1. The idea is to obtain
average models that are valid in the three regions of the system and
are coupled by continuity conditions using a simplified version of the
volume averaging method. More specifically, this is achieved by deriv-
ing macroscale equations in both the x-andg x regions in which
the effective permeability tensors are obtained from coupled closure
problems. The solution in the g-region is obtained analytically since
the flow is one-dimensional along the channel axis (e
x
)inthisregion.
It is also coupled with that in the g x inter-region. The structure of
the average models in the three regions is shown to involve a Darcy-
like momentum equation, which can be condensed into a single
one-equation model with a position-dependent apparent permeability
tensor. The solution of the closure scheme to obtain this tensor is
computationally much less demanding than performing direct
numerical simulations. The unified model is quite practical and
appealing due to its simplicity. It is developed in a rather general flow
context, in which inertial and/or slip effects at the solid/fluid interfa-
ces may be present, widely extending the scope of the derived model
with respect to the existing works reported so far. At this point, it is
worth mentioning that there are some relevant works including slip
effects in dual-porosity media;
20–22
however, these works do not
include inertial effects and do not account for the spatial variations of
the apparent permeability tensor as it is done here.
With this aim in mind, the article is organized as follows. In
Sec. II, the governing equations at the microscale are pr esented, and
the derivations of the macroscale equations in each region are
reported, yielding the unified one-domain momentum equation.
Section III is dedicated to a series of illustrative examples showing the
performance of the one -domain approach. This is achieved by com-
parisons of the flow fields obtained from direct numerical simulations
at the mi croscale with predictions of the flow from the one -domain
approach developed here. It is carried out for creeping flow in the
absence and in the presence of slip effects as well as in the presence of
inertia. In addition, validation is also presented by comp arison with
experimental data unde r creeping flow and no-slip conditions.
Conclusions are drawn in Sec. IV.
II. MICROSCALE MODEL AND DEVELOPMENT OF THE
ONE-DOMAIN APPROACH
The configur ation under con sideration is the one repre sented in
Fig. 1. Before developing the one-domain macroscale model, the
underlying boundary val ue flow problem must be formulated at the
microscale.ThisisreportedinSec.II A foll owed by a ra pid prese nta-
tion of the derivation of the average mass-equation applicable every-
where in the system (Sec. II B). Finally, the macroscopic momentum
transport equation is derived in each region (Secs. II C–II E).
A. Microscale model
The governing equations for mass and momentum transport at
the pore-scale (in the three regions) are formulated in a rather general
version in which inertia and/or slip effects may be present. Assuming
incompressible flow, the mass and momentum balance equations in
the j-region (j ¼ x; gx; g) are given by
rv
j
¼ 0; in the b-phase; (1a)
qv
j
rv
j
¼rp
j
þ lr
2
v
j
; in the b-phase: (1b)
In these two equations, q and l denote the fluid density and dynamic
viscosity, both considered as constant, whereas v
j
and p
j
represe nt the
fluid velocity and pressure in th e j-region. Without loss of generality,
volume forces are not considered in the momentum transport equa-
tion. Taking into account the possible exist en ce of rare faction effects , in
particular if the flowing fluid is a gas, due to the size of the pores and/or
channel in conjunction with the thermodynamic conditions, the inter-
facial boundary condit ion can be formulated as follows:
23–26
v
j
¼n
k
k I nn
ðÞ
n rv
j
þrv
T
j
; at A
brj
: (1c)
Here, I and n, respectively, denote the identity tensor and the unit nor-
mal vector to the solid-fluid interface, A
brj
,inthej-region, directed
from the b-phase to the r-phase. In addition, k ¼ 1forj ¼ x; gx
(i.e., at the solid–fluid interface inside the porous material) and k ¼ 2
when j ¼ g (i.e., at the top surface of the channel). The mean free
path of the fluid molecules is denoted by k while n
k
¼ð2 r
k
Þ=r
k
is
a coefficient taking into account the reflection process of the molecules
at the solid wall related to the tangential accommodation coefficient,
r
k
. At this point, it is worth mentioning that the structure of the
boundary condition given in Eq. (1c) is a Navier-type slip, and it is
expressed in its complete form, i.e., including the transpose of the
velocity gradient, which is different from previous works.
27,28
Furthermore, this type of boundary condition has been pr oposed to
study other physical mechanisms such as flow over rough surfaces.
29,30
In this way, th e above bound ary cond ition is quite rich as it allows
studying physical situations beyond only rarefaction effects. In addi-
tion, ma croscopic boundary cond itions at the inlet and outlet of the
system (in the x-direction) should be provided to complete the state-
ment of th e pore-scale problem. However, they will not be used for the
derivations that follow and are hence left unspecified here.
In order to carry out a macroscale description, an averaging
domain, V,ofmeasureV and characteristic size r
0
, is defined. Two
averaging operators are considered, namely, the superficial and
intrinsic averages. For a piece-wise smooth function, w,definedinthe
fluid phase within th e j-region, they are, respectively, given by
hwi
j
¼
1
V
ð
V
bj
w dV; (2a)
hwi
b
j
¼
1
V
bj
ð
V
bj
w dV; (2b)
with j ¼ x; gx; g.Intheaboveequations,V
bj
(of measure V
bj
)
denotes the portion o f the averaging domain occupied by the fluid
phase in the j-region. Both averaging operators are related according to
hwi
j
¼
V
bj
V
hwi
b
j
: (2c)
The ratio V
bj
=V denotes the volume fraction of the fluid phase within
the averaging domain and it is, in general, a function of position. In
the porous medium bulk, this ratio is a constant corresponding to the
porosity, e.
B. Average mass balance equation
The development starts with the derivation of the macroscopic
mass balance equation in each region. Application of the superficial
averaging operator to the continuity equation (1a) and use of the spa-
tial averagin g theorem
31
hr ai
j
¼rhai
j
þ
1
V
ð
A
brj
n a dA; j ¼ x; gx; g (3)
(a denoting a vector field defined in the b-phase), together with the
interfacial boundary cond ition [Eq. (1c)], leads to the macrosc opic
mass conservation equation which takes the same form in the three-
regions. It is given by
rhv
j
i
j
¼ 0; in the j region; j ¼ x; gx; g: (4)
The attention must now be focused on the average momentum equa-
tion in each region and this is the purpose of Secs. II C–II E.
C. Analysis in the homogeneous part of the porous
medium (x-region)
This portion of the syste m has received cons iderable attention in
the literature; hence, the derivations provided here are presented in a
brief manner. More details can be found in references dedicated to the
averaging of one-phase flow in a porous medium in the presenc e of
inertia
32,33
and to the de rivation of a macroscopic model for one-
phase flow when slip is present at the solid–fluid interfaces.
26,34
For
the sa ke of convenience, the str ucture of the porous medium is mod-
eled as a periodic array of solid inclusions as illustrated in Fig. 2.
Hence, the averaging domain, V, may be defined as a finite arra y of
unit cells (each of size ‘
c
), with the following constraint for its charac-
teristic size, r
0
:
‘ r
0
L: (5)
Here, ‘ represents the largest characteristic size at the pore-scale [i.e.,
‘ ¼ maxð‘
b
;‘
r
Þ], while L is the smallest characteristic size of the
porous medium at the macroscale (i.e., the smallest dimension of the
porous medium in the three directions of space).
In order to derive the average momentum transport equation, it
is convenient to follow the short-cut approach suggested by Barre`re
et al.,
35
in which Eqs. (1) are considered in th e avera gin g do main, con-
ceived as a Representative Elementary Volume (REV), ma de of one or
more periodic unit cell(s) as sketched in Fig. 2. In this doma in, it is re a-
sonable to decompose the fluid pressure gradient according to
36
rp
x
¼rhp
x
i
b
x
þr
~
p
x
; (6)
~
p
x
being the pressure deviations, and regard rhp
x
i
b
x
as a consta nt
within the REV. This assumption is justified on the basis of the
separation of length scales given in (5). In this way, the pressure devia-
tions and the fluid velocity can be assumed to be periodic at the inlets
and outlets of the REV, and the pore-scale model can be written as
follows:
rv
x
¼ 0; in V
bx
; (7a)
qv
x
rv
x
¼r
~
p
x
þ lr
2
v
x
rhp
x
i
b
x
; in V
bx
; (7b)
v
x
¼n
1
k I nn
ðÞ
n rv
x
þrv
T
x
; at A
brx
; (7c)
wðr þ l
i
Þ¼wðrÞ; i ¼ 1; 2; 3; w ¼ v
x
;
~
p
x
; (7d)
h
~
p
x
i
b
x
¼ 0: (7e)
The last equation is an average constraint bounding the field of the
pressure deviations that is compliant with the assumption that hp
x
i
b
x
is constant within the REV, and it is necessary for the flow problem to
be well-pose d. In Eq. (7d), l
i
(i ¼ 1, 2, 3) repr esents the periodic lattice
vectors of the REV.
The formal solution of th i s problem can be written as
v
x
¼
F
l
rhp
x
i
b
x
; in V
bx
; (8a)
~
p
x
¼f rhp
x
i
b
x
; in V
bx
; (8b)
where F and f are closure variables that resu lt from solving the follo w-
ing problem in the (perio dic) REV
rF ¼ 0; in V
bx
; (9a)
q
l
2
jjrhp
x
i
b
x
jje
p
F
T
rF ¼rf þr
2
F þ I; in V
bx
; (9b)
F ¼n
1
k I nn
ðÞ
n rF þrF
ðÞ
T1
; at A
brx
; (9c)
wðr þ l
i
Þ¼wðrÞ; i ¼ 1; 2; 3; w ¼ F; f ; (9d)
hfi
b
x
¼ 0: (9e)
In Eq. (9b), e
p
is a unit vector in the direction of the macroscopic pres-
sure gradient. It is defined as
e
p
¼
rhp
x
i
b
x
jjrhp
x
i
b
x
jj
: (10)
In addition, in Eq. (9c), the sup erscript T1 denotes the transpose oper-
ator that permutes the first and second indices of a third order tensor
T,namely,T
T1
ijk
¼ T
jik
.
FIG. 2. Sketch of the averaging domain of size r
0
and a periodic two-dimensional
unit cell (of size ‘
c
) for the analysis of the flow problem in the homogeneous part of
the porous medium. Note that, in general, r
0
‘
c
. The solid phase, r,is
represented as a random distribution of circular inclusions embedded in the fluid
phase, b.
