HAL Id: hal-03042108
https://hal.archives-ouvertes.fr/hal-03042108
Submitted on 5 Dec 2020
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.
Comparison between direct numerical simulations and
eective models for uid-porous ows using penalization
Charles-Henri Bruneau, Didier Lasseux, Francisco J Valdés-Parada
To cite this version:
Charles-Henri Bruneau, Didier Lasseux, Francisco J Valdés-Parada. Comparison between direct nu-
merical simulations and eective models for uid-porous ows using penalization. Meccanica, Springer
Verlag, 2020, 55 (5), pp.1061-1077. �10.1007/s11012-020-01149-7�. �hal-03042108�
Comparison between direct numerical simulations and
effective models for fluid-porous flows using penalization
Charles-Henri Bruneau
1
, Didier Lasseux
2
, Francisco J. Valdés-Parada
3
Corresponding author: bruneau@math.u-bordeaux.fr
1
Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33405, Talence, France
INRIA, IMB, UMR 5251, F-33405, Talence, France
2
CNRS, I2M, UMR 5295, Esplanade des Arts et Métiers, 33405 Talence, France
3
Departamento de Ingeniería de Procesos e Hidráulica, División de Ciencias Básicas e
Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa,
Av. San Rafael Atlixco 186, 09340, CDMX, Mexico
Abstract: This work is devoted to a numerical study of two-dimensional incompressible flows in a
fluid-porous medium system placed on an impermeable wall. Flow in the whole system is studied
either at the pore scale or at the Darcy scale (using a one-domain approach) and the models at
both scales are solved with a penalization method using the same formal Navier-Stokes equations
modified by a Darcy-like term. Several effective medium penalized models are considered for the
simulations. Various flow regimes are investigated ranging from laminar to turbulent. The velocity
profiles inside and outside the porous medium show significant discrepancies between the different
penalization models compared to the direct numerical simulations. This work motivates further
studies about the spatial variations of the penalization coefficient to be introduced in the effective
medium models in order to better reproduce the physics near the fluid-porous medium boundaries.
Keywords: Fluid-porous media, Direct numerical simulation, Penalization method, Effective
medium models.
1 Introduction
Momentum transport between a porous medium and a fluid has been a topic of intense research interest
over, at least, half a century, starting with the pioneering theoretical and experimental study by Beavers
and Joseph (BJ) [6], who analyzed creeping flow in a channel partially filled with a porous medium. These
authors proposed to couple Darcy’s law with the Stokes equation by means of a Robin-type matching
boundary condition. Shortly after, Saffman [36], used a statistical approach to extend Darcy’s law outside
of the porous medium bulk. Then, the boundary condition is recovered from a boundary-layer analysis
imposing a step function for the spatial variations of the permeability and porosity. More recently, this
boundary condition has been derived using the homogenization technique [20] and the volume-averaging
method [43]. However, numerical simulations performed at the pore-scale have shown that the slip velocity
in the BJ model is extremely sensitive to the position of the fluid-porous medium boundary [26, 27].
One alternative to the use of the BJ model is to match the Darcy-Brinkman and Stokes equations by
means of two boundary conditions. In this regard, Ochoa-Tapia and Whitaker [34] proposed that the velocity
should be continuous (which contrasts with the BJ statement) and that there should be a jump in the viscous
stress at the boundary. These ideas may be viewed as an extension of those by Neale and Nader [32] who used
continuity conditions of both the velocity and stress in order to couple the Darcy-Brinkman model to Stokes
equations. Other alternatives of matching conditions between these two equations have been proposed. For
example, Cieszko and Kubik [16] used two linearly-independent BJ-type equations to carry out the matching,
which may be interpreted as a discontinuity on both the velocity and the stress. However, a point against
these works is that the validity of the Darcy-Brinkman model is questionable, in general, as pointed out by
Auriault [5].
1
The above works can be classified as two-domain approaches because they rely on the matching between
the governing equations in the porous medium with those in the fluid by means of convenient boundary
conditions. However, from the above references, it appears that no definitive answer has been proposed in
this modeling approach. Conversely, the one-domain approach makes use of a single average equation, with
position-dependent coefficients, which can be used to model transport in the porous medium, the fluid and in
the transition zone [4]. This has the advantage that no boundary conditions are required at the fluid-porous
interface, at the expense of requiring knowledge of the spatial variations of the effective-medium coefficients
such as the permeability (in the case of creeping flow) or the apparent-permeability (in the case of inertial
or even turbulent flow). The prediction of these spatial variations is still an active research field in which
sharp [21] or smooth [42] variations of the permeability have been proposed for creeping flow conditions.
Furthermore, the one-domain approach can be used in conjunction with the differential equations in the
two-domain approach to derive the corresponding boundary conditions [34] and to predict their coefficients
[42]. For this reason, the one-domain approach is used in this work.
Literature on flow between a porous medium and a fluid has been mainly dedicated so far to creeping
flow in a channel partially filled with a porous medium (BJ configuration) and extensions to the inertial (or
turbulent) regime are still scarce. A complete one-domain averaged (but unclosed) model valid everywhere
and including inertia was derived in [33] and recently revisited in [35] (see Eqs. (17) and (18) therein).
In the former reference, a closure scheme was proposed for a two-domain approach while ignoring inertia
in the homogeneous part of the porous region. The unclosed one-domain model could be used here as
this is the most complete description available so far. Nevertheless, the solution of this model is still an
important challenge as it requires to predict the porosity gradient in the interfacial region and the pore-scale
pressure and velocity fields to close the model. This is obviously not a computational procedure which can
be envisaged as a routine solution to the problem. An alternative to circumvent this difficulty is to use an
empirical closure scheme as proposed in [8, 9, 37] for the BJ configuration. In these works, an empirical
form of the spatial variation of the porosity was assumed. It was used in empirical correlations between the
porosity and the permeability and inertial resistance coefficient usually employed for the homogeneous part
of the porous medium under consideration.
Our purpose in the present work is to investigate whether an effective-medium model, derived for a
homogeneous porous medium (i.e. far from the interfaces with a fluid region), which is much simpler than
the complete one-domain averaged model, remains accurate when employed in the whole domain using a
penalization approach. This effective-medium model can be formally viewed as a Navier-Stokes equation
written in terms of spatial averages of the velocity and pressure with addition of a Darcy-like penalization
term. This is performed on a system different from the BJ configuration with the idea that, in many practical
situations, the mean flow is not fully developed along the fluid-porous medium interface and that this interface
is not necessarily parallel to the mean flow. The configuration under consideration here is similar to those
used to study bluff bodies [1]. It is made of a two-dimensional rectangular porous slab consisting of a square
pattern of parallel cylinders immersed in a semi infinite fluid region limited by a solid wall. The aim is
first to accurately compute the incompressible flow orthogonal to the cylinders’ axes inside and outside the
porous medium. This is achieved by a Direct Numerical Simulation (DNS) of the penalized Navier-Stokes
equations. The DNS serves as a reference for further comparisons with the effective homogeneous model
and has archival value for subsequent studies in this particular configuration. In a second step, the porous
zone is considered as a homogeneous medium taking into account its properties, namely, the permeability
coefficient and porosity. Three different forms of the Darcy-like penalization term in the effective-medium
model are investigated. Ultimately, the purpose is to show that the flow solution for a fluid-porous system
obtained with the effective model using a penalization approach on a coarse grid is much faster than with
DNS.
Both DNS and the solution of the effective medium models are performed using a penalization approach
as it will be explained in further details below. An appealing feature of this approach is that the formulation
of the momentum equations is the same (although pressure and velocities do not have the same scale for
the DNS and effective medium simulations). The only difference lies in the way penalization is applied to
the Darcy-like term. The penalized equations are approximated by an accurate finite differences scheme and
solved using a multigrid procedure involving several grid levels. The code is highly parallelized with Message
Passing Interface (MPI) directives.
Models and results are carefully analyzed to determine which effective-medium model yields results
2
100`
c
500`
c
600`
c
˜
˜
˜
˜
12,000 obstacles: Ω
s
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . . . . .
ε = 0.913
Fluid region: Ω
f
600`
c
30`
c
x
∗
y
∗
`
c
a
a)
100`
c
500`
c
600`
c
˜
˜
˜
˜
Equivalent homogeneous medium: Ω
p
ε = 0.913
Fluid region: Ω
η
600`
c
30`
c
x
∗
y
∗
b)
Figure 1: Domain configuration for a) the DNS and b) the effective medium approach.
closest to those obtained from DNS. In particular the flow fields inside the porous medium are scrutinized
to investigate the relevance of the different effective medium models close to and far from the porous zone.
The following sections are organized as follows. In section 2, the pore scale and heuristic effective
medium models are reported. The use of the penalization procedure for both models is explained. Section
3 is dedicated to the description of the numerical method. The numerical results corresponding to DNS
are compared with those obtained with the effective medium approach in section 4. In Appendix B, it is
shown that some improvement to the effective medium approach can be obtained empirically. Conclusions
are provided in section 5.
2 Modeling
The aim of this work is to compare as accurately as possible the DNS of an unsteady two-dimensional
Newtonian one-phase flow inside and around a porous medium to heuristic versions of effective medium
models using the penalization approach. To do so, consider a domain Ω partially filled with a regular
pattern of in-line parallel cylinders of square cross-section as the one depicted in Fig. 1. The pattern is
obtained from 400 x
∗
- and 30 y
∗
-periodic repetitions of a unit cell of size `
c
× `
c
(see Fig. 1a).
The solid cylinders Ω
s
, have an edge size a = 0.295`
c
, yielding a porosity of the unit cell = 0.913.
The entire domain, saturated by the fluid occupying the region Ω
f
, has a size 600`
c
× 600`
c
and the porous
zone, corresponding to the portion occupied by the cylinder pattern, is placed at 100 6 x
∗
/`
c
6 500,
0 6 y
∗
/`
c
6 30. The lower wall is impermeable, while the vertical size of the domain is large enough to avoid
the influence of the upper boundary condition due to the blockage as shown in [15]. The flow approaches the
3
system from the left boundary with a constant velocity v = u
∞
e
x
, while at the exit, on the right boundary,
the traction is specified following a formulation provided in [12].
The aim is first to compute a reference flow using the DNS around the 12, 000 squares by solving the
penalized Navier-Stokes equations for the genuine unknowns velocity u, and pressure p, in the fluid domain
Ω
f
. For the sake of simplicity and in order to use a formulation that is uniform for both the DNS and
the computation of the effective model, the flow is solved in the whole domain Ω using the penalization
method [2, 25, 10]. This allows to easily handle the 12, 000 solid squares without the need of an adapted
mesh. The advantages and disadvantages of this method are discussed in detail in[31]. It should be pointed
out that the volume (or Brinkman) penalization method has been extensively used over the past twenty
years to solve problems in complex geometries [10]. It is particularly suited for all kind of spectral or finite
differences approximations [24, 23] but is also used with other approximations [39, 30, 40]. The method is
applied directly to simulate incompressible flows or is extended to other situations including moving solid
obstacles [22] and heat transfer [3] for instance. The aim is always to easily handle fixed or moving complex
geometries [19, 38, 18, 7]. Authors found agreement with laboratory experiments or with benchmark results
obtained with unstructured meshes.
Under these circumstances, using `
c
, u
∞
, `
c
/u
∞
and ρu
2
∞
as the references for the lengths, velocity, time
and pressure, respectively, the model can be written in the dimensionless form as
∂
t
u + u · ∇u = −∇p +
1
Re
∇
2
u − Γ · u in Ω
T
(1a)
∇ · u = 0 in Ω
T
(1b)
B.C.1 u = 0 at 0 6 x 6 600, y = 0 (1c)
B.C.2 ∇u · n = 0 at 0 6 x 6 600, y = 600 (1d)
B.C.3 u = e
x
at x = 0, 0 6 y 6 600 (1e)
B.C.4 Σ · n +
1
2
(u · n)
−
(u − u
ref
) = Σ
r ef
· n at x = 600, 0 6 y 6 600 (1f)
I.C. u = 0 in Ω at t = 0 (1g)
In these equations, u, p, t, x and y represent the dimensionless velocity, pressure, time and spatial
coordinates, respectively. The Reynolds number of the flow is given by Re = ρu
∞
`
c
/µ and Ω
T
= Ω × (0, T ),
with T being the simulation time. The traction boundary condition on the exit section involves the stress
tensor Σ and uses the notation a = a
+
− a
−
. It yields a well-posed problem [12] and conveys properly the
vortices through the artificial boundary. It requires a reference flow that is selected as the calculated flow
on the front of the last column of cells [11]. Finally, Γ · u is the penalization term, which for the DNS, is
fixed to be Γ = 10
−16
I so that it vanishes in Ω
f
. On the opposite, Γ = 10
8
I in Ω
s
so that the velocity
becomes of the order 10
−8
inside each square. It is worth mentioning that the penalized model in Eqs. (1)
is also used to solve the effective-medium equations that are taken to be formally the same as the pore scale
model. This represents a heuristic approach since a rigorous derivation of the effective equations is still an
active research field. A recent work proposed an up-scaled model that is non-local in time together with the
associated closure problems [29]. The up-scaled momentum equation does not correspond to Eq. (1a), but,
for the sake of simplicity, the heuristic model of Eqs. (1) is used here as a first approach. For this model,
u and p must be understood as the intrinsic average of the point-wise velocity and pressure. The intrinsic
average of any physical quantity ψ is defined by hψi
β
= 1/V
β
R
V
β
ψdV , V
β
being the region of the fluid phase
(of volume V
β
) within a periodic unit cell of the porous medium. Consequently, the entire porous medium
is conceived as a pseudo-continuum (Ω
p
) as sketched in Fig. 1b. Following the same trend of ideas used for
the DNS, Γ = 10
−16
I in the fluid region, which is now denoted as Ω
η
.
The simplest effective-medium model is obtained by fixing Γ =
Re
k
−1
I, where k is the dimensionless
intrinsic permeability of the porous medium, which is obtained from the solution of an intrinsic closure
problem as detailed below. In the following, this model is referred to as the K model. A second-order model,
denoted as the K2 model, introduced in [41] following the experiments in [17], is also employed using an
Ergun formulation given by Γ =
Re
k
−1
+ 0.1429
3/2
k
−1/2
kuk
I. In this formulation, k can be estimated
from Ergun’s formula or from the closure procedure. Both have been tested and, since the two values of k
4
