1 Introduction
Flows in porous media are common in nature and in technical applications, with velocities through the pores which range from very low values (e.g. flow in biological tissues, ground water infiltration) to very large (e.g. canopy flow). When a fluid flows in a porous medium, there is usually a strong separation of scales between the macroscopic scale defined by the global size of the problem and the microscopic length typical of the solid inclusions of the porous medium (cf. figure 1). Two approaches are generally used to treat this type of phenomenon: the first one is to implement a pore-scale numerical simulation of the flow in the medium, reproducing closely the geometry of the solid skeleton. The second consists in an effective macroscopic simulation where the microscopic structure is lost, but some auxiliary problem is introduced to characterize the pore-scale behaviour. The first path has been followed, for instance, by Breugem & Boersma (Reference Breugem and Boersma2005), Kuttanikkad (Reference Kuttanikkad2009) and Matsumura & Jackson (Reference Matsumura and Jackson2014). One of the difficulties of this kind of analysis is the computational cost of the simulation: the smaller is the ratio between the microscopic and the macroscopic length scale, the larger is the resolution which must be used to discretize the geometry of the skeleton. Such a big effort presents the additional disadvantage of providing a solution to only a particular configuration of the porous medium. With the second choice we have a concise and rapid description of the fluid behaviour, satisfactory from a macroscopic point of view. A precursory idea of this method consists in assuming that simulating a fluid flowing through a porous medium is the same as thinking of the coupled fluid–solid medium as a continuum for which an effective conductivity, permeability or viscosity, different from that of the fluid, can be defined. These effective medium approaches were first postulated from an empirical point of view by Darcy (Reference Darcy1856) and then modified using analytical considerations (Brinkman Reference Brinkman1949). Even if the equations were born as empirical laws, several strategies have been developed in time in order to derive these equations analytically, deducing the macroscopic behaviour from the local description (see Davit et al.
Reference Davit, Quintard, Bell, Byrne, Chapman, Kimpton, Lang, Oliver, Pearson and Waters2013, for a review). This procedure is called upscaling. Different techniques are available for upscaling: they start from a representative elementary volume (REV) and generate an equivalent macroscopic continuous model called the homogenized model. One of the most used approaches is the volume averaging method thoroughly described by Whitaker (Reference Whitaker1998). It consists in considering an elementary cell representative of the porous structure (
$V$
in figure 1); the size of the REV is of the order of the pore size. A fundamental assumption is that of periodicity, over
$V$
, for the unknown quantities. Starting from the Navier–Stokes equations (NSE), it is assumed that the solution (velocity and pressure fields) can be decomposed into a mean part plus a fluctuation (Gray Reference Gray1975). The fluctuation represents how much the solution is far from its mean part and it must satisfy a zero average condition over the REV, that is fundamental to deduce the effective equations. After substituting the decomposition of the unknown fields inside the NSE, and considering the spatial average of the equations over the REV, one obtains the momentum equation for the mean flow, with forcing terms related to the fluctuations. In this way, the Darcy’s law (Whitaker Reference Whitaker1986), its Brinkman’s correction (Quintard & Whitaker Reference Quintard and Whitaker1994) and the Forchheimer’s equation (Whitaker Reference Whitaker1996) have been deduced theoretically.
Figure 1. View of a generic porous medium within an elementary cell
$V$
.
$V_{f}$
is the volume occupied by the fluid and
$V_{s}$
is that occupied by the solid, so that
$V=V_{f}+V_{s}$
.
${\it\Gamma}$
is the fluid–solid microscopic interface.
Another strategy widely used is described by Mei & Vernescu (Reference Mei and Vernescu2010) and consists in implementing a homogenization technique based on a multiple scale analysis. The starting point of this technique is the same as the volume averaging method: the NSE are taken to be valid over
$V_{f}$
in figure 1 and a proper expansion of the unknown fields is assumed. If the convective term is sufficiently small (of the order of the ratio between length scales) the resulting equations at leading order are self-contained, i.e. there is no need to introduce a closure relation. Darcy’s, Brinkman’s and Forchheimer’s equations can be deduced analytically also in this case. Both approaches briefly outlined above require closure relations when nonlinear terms are present. Much work is present in the literature: in particular, Mei & Auriault (Reference Mei and Auriault1991) have examined the effect of weak fluid inertia in the porous medium, finding that Forchheimer’s correction to Darcy’s law should be at least cubic (instead of quadratic) for isotropic media. This has been confirmed by Firdaouss, Guermond & Le Quéré (Reference Firdaouss, Guermond and Le Quéré1997). The same occurs for orthotropic media (Skjetne & Auriault Reference Skjetne and Auriault1999). Auriault (Reference Auriault2009) investigated the domain of validity of Brinkman’s equation, finding that it is valid for swarms of fixed particles or a fixed bed of fibres at very low concentration and under precise conditions which depend on the separation of the scale parameter
${\it\epsilon}=l/L$
, with
$l$
and
$L$
two different representative length scales.
Homogenization provides a point of contact between the microscopic and the macroscopic worlds, and permits to transfer information from one point of view to the other. In this, it differs from purely microscopic approaches such as those by Jackson & James (Reference Jackson and James1986), van der Westhuizen & du Plessis (Reference van der Westhuizen and du Plessis1996), Tamayol & Bahrami (Reference Tamayol and Bahrami2009) and Yazdchi et al. (Reference Yazdchi, Srivastava and Luding2011), and it is also different from purely macroscopic points of view, such as that by Battiato (Reference Battiato2012) in which analytical forms of the permeability are taken from the literature.
All the works and the methods cited up to now refer to the solution of the flow deeply inside a porous medium, far from the boundaries. The problem of the interface conditions between the pure fluid and the porous region is crucial and amply discussed in the literature. Jäger & Mikelić (Reference Jäger and Mikelić1996, Reference Jäger and Mikelić2000) spent much effort on this problem, concluding, through functional analysis and introducing an auxiliary problem based on homogenization, that the widely used condition of Beavers & Joseph (Reference Beavers and Joseph1967) and its modification by Saffman (Reference Saffman1971) have a mathematical justification. Beavers & Joseph (Reference Beavers and Joseph1967) observed that the penetration of the velocity into the porous bed extends over a length proportional to
$\sqrt{\mathscr{K}}$
, with
$\mathscr{K}$
the bed’s permeability. This is equivalent to specifying a jump in the average velocity at the interface. More elaborate conditions have also been proposed, simulating the flow inside the porous medium with the Darcy–Brinkman equation, introducing an effective viscosity
${\it\mu}_{e}$
(which can be evaluated experimentally, see Givler & Altobelli Reference Givler and Altobelli1994), and imposing continuity of the normal and tangential components of the stress tensor at the interface (Hill & Straughan Reference Hill and Straughan2008). Another interface condition has been developed by Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995) via the definition of the excess surface and bulk stress tensors to be obtained from the governing equations holding in the porous and pure fluid regions. This condition becomes, under certain hypotheses, the continuity of the effective velocity and pressure over the macroscopic interface. Finally, a strategy often used to couple two different media is the penalization method. It consists in solving the NSE in the whole domain, with a forcing term added in the porous region to take into account the presence of the structure. This forcing term is multiplied by the characteristic equation of the porous domain and it goes to zero smoothly outside of it. This approach has been pursued, among others, by Angot et al. (Reference Angot, Bruneau and Fabrie1999) and Bruneau & Mortazavi (Reference Bruneau and Mortazavi2008). Cimolin & Discacciati (Reference Cimolin and Discacciati2013) have shown that even if the penalization method is stable and easy to implement, the comparison of the solution against experimental data near the interface is not satisfactory.
Summing up, the interface conditions can be essentially classified into three classes: the first one involves pressure and velocities which are linked over the interface directly (continuity or jump); the second one involves them indirectly, linking the normal to the interface components of the stress tensor. The second class cannot be employed if the equation to be solved in the porous medium is Darcy’s law. The third class includes all those methods which use a filter to go from the porous region to the fluid region; a drawback of this latter approach is that there is no general physical justification for the choice of the filter. Some authors (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2004; Jamet & Chandesris Reference Jamet and Chandesris2009) have proven that the variation of the permeability near the interface is not an intrinsic property of the porous medium, but depends on the properties of the flow.
In the first part of the present paper, using a homogenization technique, a strategy to develop models for flow through a bed of dense, rigid fibres is explained (figure 2). Particular attention is paid to the flow regime, so that two models are deduced: a classical one for Stokes’ flows and an original one which accounts for inertial effects. For each regime two sets of equations apply: one microscopic and one macroscopic. Solving both sets of equations, a complete description of the flow inside the porous medium is obtained. In the second part of the work the models are validated on the basis of numerical simulations which solve for the fluid flow through a real, porous geometry (constituted by fibres), and against experiments of turbulent flows through a canopy (Ghisalberti & Nepf Reference Ghisalberti and Nepf2004, Reference Ghisalberti and Nepf2006, Reference Ghisalberti and Nepf2009).
Figure 2. View of the particular porous medium studied in this work. The solution is assumed
$V$
-periodic, on the elementary cell defined by dashed lines.
2 Multiple scale analysis
Let us focus, for the moment, only on the flow inside a porous medium, neglecting the interface with the pure fluid domain. We consider an unbounded, rigid, porous medium saturated by an incompressible Newtonian fluid of constant density
${\it\rho}$
and dynamic viscosity
${\it\mu}$
. The velocity and pressure fields in the fluid domain
$V_{f}$
are ruled by the incompressible Navier–Stokes equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}}{\partial x_{i}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Re\frac{\partial u_{i}}{\partial t}+Reu_{j}\frac{\partial u_{i}}{\partial x_{j}}=-\frac{1}{{\it\epsilon}}\frac{\partial p}{\partial x_{i}}+{\rm\nabla}^{2}u_{i}, & \displaystyle\end{eqnarray}$$
with
$u_{i}=0$
on the solid–fluid interface
${\it\Gamma}$
. The equations above have been rendered dimensionless assuming that the global pressure gradient is balanced by the local viscous term, i.e.
$$\begin{eqnarray}O\left(\frac{{\rm\Delta}P}{L}\right)\sim O\left(\frac{{\it\mu}U}{l^{2}}\right),\end{eqnarray}$$
which defines the pressure scale, and time is normalized by
$l/U$
, where
$l$
is a microscale length which represents a characteristic dimension of the solid inclusions and
$L$
is a macroscale associated to the global pressure gradient
${\rm\Delta}P$
. In this way the ordering symbol
${\it\epsilon}$
is defined as
$$\begin{eqnarray}{\it\epsilon}=\frac{l}{L}\ll 1,\end{eqnarray}$$
and the Reynolds number, based on microscale length
$l$
and reference velocity
$U$
, is
$Re=Ul/{\it\nu}$
, with
${\it\nu}={\it\mu}/{\it\rho}$
the kinematic viscosity.
2.1 The expansion in
${\it\epsilon}$
Since the structure of the medium defines two characteristic scales, we need to introduce the fast (microscopic) and slow (macroscopic) variables
$\boldsymbol{x}$
and
$\boldsymbol{x}^{\prime }={\it\epsilon}\boldsymbol{x}$
with the subsequent change in the derivative:
$$\begin{eqnarray}\frac{\partial }{\partial x_{i}}\rightarrow \frac{\partial }{\partial x_{i}}+{\it\epsilon}\frac{\partial }{\partial x_{i}^{\prime }},\end{eqnarray}$$
plus the expansion:
$$\begin{eqnarray}f=f^{(0)}+{\it\epsilon}f^{(1)}+\cdots \,,\end{eqnarray}$$
where
$f$
represents a generic variable of the unknown field, function of
$(\boldsymbol{x},\boldsymbol{x}^{\prime },t)$
. This choice highlights the dual nature of the problem, i.e. the partition between microscopic and macroscopic aspects of the problem connected by the fundamental concept of averaging (cf. the next section). Once we have substituted the expansion (2.6) into the governing equations, grouping each term in orders of
${\it\epsilon}$
and considering only those of order
$0$
and
$1$
, we obtain the following equations defined over
$V_{f}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}^{(0)}}{\partial x_{i}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}^{(1)}}{\partial x_{i}}+\frac{\partial u_{i}^{(0)}}{\partial x_{i}^{\prime }}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\partial p^{(0)}}{\partial x_{i}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Re\left(\frac{\partial u_{i}^{(0)}}{\partial t}+u_{j}^{(0)}\frac{\partial u_{i}^{(0)}}{\partial x_{j}}\right)=-\frac{\partial p^{(1)}}{\partial x_{i}}-\frac{\partial p^{(0)}}{\partial x_{i}^{\prime }}+\frac{\partial ^{2}u_{i}^{(0)}}{\partial x_{j}\partial x_{j}}. & \displaystyle\end{eqnarray}$$
Equation (2.9) yields
$p^{(0)}=p^{(0)}(\boldsymbol{x}^{\prime },t)$
, i.e. the leading-order term of the pressure varies only over the macroscale (and time). Equation (2.10) contains the inertial terms which may, or may not, be negligible depending on the order of magnitude of the microscopic Reynolds number,
$Re$
.
2.2
$Re=O({\it\epsilon})$
If the Reynolds number is of order
${\it\epsilon}$
, the left-hand side of (2.10) goes to higher order, so that, following Mei & Vernescu (Reference Mei and Vernescu2010), the solution can be written formally as
$$\begin{eqnarray}u_{i}^{(0)}=-\unicode[STIX]{x1D612}_{ij}\frac{\partial p^{(0)}}{\partial x_{j}^{\prime }},\quad p^{(1)}=-A_{j}\frac{\partial p^{(0)}}{\partial x_{j}^{\prime }}+p_{0}^{(1)}(\boldsymbol{x}^{\prime },t),\end{eqnarray}$$
where
$\unicode[STIX]{x1D612}_{ij}$
is a tensor,
$A_{j}$
is a vector and
$p_{0}^{(1)}$
is an integration constant (with respect to the integration variable
$\boldsymbol{x}$
). Substituting (2.11) into (2.7) and (2.10) and simplifying, it follows that the coefficients should satisfy the relations:
$$\begin{eqnarray}\frac{\partial \unicode[STIX]{x1D612}_{ij}}{\partial x_{i}}=0,\quad -\frac{\partial A_{j}}{\partial x_{i}}+{\rm\nabla}^{2}\unicode[STIX]{x1D612}_{ij}=-{\it\delta}_{ij};\end{eqnarray}$$
this is a forced Stokes problem with the unknowns
$\unicode[STIX]{x1D612}_{ij}$
and
$A_{j}$
for
$i,j=1,2,3$
(in principle a system of 12 equations and 12 unknowns) with boundary condition
$$\begin{eqnarray}\unicode[STIX]{x1D612}_{ij}=0\quad \text{on }{\it\Gamma},\end{eqnarray}$$
plus periodicity over
$V$
for
$\unicode[STIX]{x1D612}_{ij}$
and
$A_{j}$
. We note that the solution of (2.12), (2.13) depends only on the microscopic spatial scale. Moreover, to ensure uniqueness of the solution we can impose
$$\begin{eqnarray}\langle A_{j}\rangle =0,\end{eqnarray}$$
for each
$j$
, where the volume average over a unit cell
$\langle \cdot \rangle$
is defined by:
$$\begin{eqnarray}\langle f\rangle :=\frac{1}{V}\int _{V_{f}}f\,\text{d}V.\end{eqnarray}$$
By taking the volume average of (2.11), noting that
$p^{(0)}$
and
$p_{0}^{(1)}$
do not depend on
$\boldsymbol{x}$
and using (2.14), we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \langle u_{i}^{(0)}\rangle =-\mathscr{K}_{ij}\frac{\partial p^{(0)}}{\partial x_{j}^{\prime }},\quad \text{with}~\mathscr{K}_{ij}=\langle \unicode[STIX]{x1D612}_{ij}\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \langle \,p^{(1)}\rangle ={\it\vartheta}p_{0}^{(1)},\quad \text{with}~{\it\vartheta}=\frac{V_{f}}{V}, & \displaystyle\end{eqnarray}$$
where
$\mathscr{K}_{ij}$
is the dimensionless permeability tensor, the dimensional counterpart of which scales like
$l^{2}$
, and
${\it\vartheta}$
the porosity. Equation (2.16) is Darcy’s law, which is a first-order approximation of the Navier–Stokes equations in
${\it\epsilon}$
. We thus have a macroscopic problem formed by the volume average of (2.8) which, after application of the spatial averaging theorem (Mei & Vernescu Reference Mei and Vernescu2010), leads to
$$\begin{eqnarray}\frac{\partial \langle u_{i}^{(0)}\rangle }{\partial x_{i}^{\prime }}=0,\end{eqnarray}$$
plus (2.16), in which the tensor
$\mathscr{K}_{ij}$
, calculated through (2.12), appears. Darcy’s law is based on the assumption that
$Re=O({\it\epsilon})$
. In view of the practical developments of this work, we now want to understand how to treat the problem when
$Re$
is not infinitesimal. One obvious extension consists in taking higher-order terms in the expansion (2.6), as done in Mei & Vernescu (Reference Mei and Vernescu2010). This leads to Darcy–Forchheimer’s equation, valid in the range
${\it\epsilon}<Re<1$
. In § 2.3 an alternative is proposed.
2.2.1 Microscopic simulations
In order to solve (2.16), we first need to calculate the permeability tensor over the normalized domain shown in figure 3. For both cases
$Re=O({\it\epsilon})$
and
$Re=O(1)$
an opensource solver based on OpenFOAM (www.openfoam.com) is used. A time-marching solver is employed to reach a steady solution. The differential operators are discretized using the finite volume method and the Gauss scheme is employed to express the divergence, gradient and Laplacian operators; we have verified that the results reported here are grid converged.
Figure 3. Isocontours of
$\unicode[STIX]{x1D612}_{11}$
in the microscopic domain for the case of porous media made by packed spheres (a) and cylindrical fibres (b).
We have first computed the case of an isotropic porous medium formed by packed spheres; our results for
$\hat{\mathscr{K}}=\langle \unicode[STIX]{x1D612}\rangle$
match perfectly those obtained by Zick & Homsy (Reference Zick and Homsy1982) for varying values of the porosity
${\it\vartheta}$
(cf. figure 4). The values of Zick & Homsy (Reference Zick and Homsy1982) are Stokes flow results, obtained through the technique of integral equations, for different packings of uniform spheres. The much-quoted Kozeny–Carman empirical formula (Carman Reference Carman1939):
$$\begin{eqnarray}\hat{\mathscr{K}}=\frac{1}{5}\left(\frac{V_{s}}{|{\it\Gamma}|}\right)^{2}\frac{{\it\vartheta}^{3}}{(1-{\it\vartheta})^{2}},\end{eqnarray}$$
is close to the theoretical findings only for
${\it\vartheta}$
around 0.5; in the formula,
$|{\it\Gamma}|$
is the interface area and
$V_{s}$
the solid volume. A better fit through the data, to represent the isotropic permeability as a function of the porosity, is given by:
$$\begin{eqnarray}\hat{\mathscr{K}}=\frac{1}{5}\left(\frac{V_{s}}{|{\it\Gamma}|}\right)^{2}\frac{{\it\vartheta}^{5/2}}{(1-{\it\vartheta})^{47/30}}\end{eqnarray}$$
(dashed red line in figure 4), valid in the range
$0.476<{\it\vartheta}<1$
.
The second step is to find the permeability of our particular porous medium and its orthotropic structure reduces to two the components of the tensor to be computed:
$\unicode[STIX]{x1D612}_{11}$
and
$\unicode[STIX]{x1D612}_{33}$
represent the longitudinal and transversal components of the permeability, respectively. Once
$\unicode[STIX]{x1D612}_{11}$
and
$\unicode[STIX]{x1D612}_{33}$
are found, we have a pointwise solution inside the elementary cell
$V$
; as before, we carry out the weighted integral of
$\unicode[STIX]{x1D612}_{ij}$
over its microscopic domain, as by (2.15), to obtain
$\mathscr{K}_{11}$
and
$\mathscr{K}_{33}$
. The results are displayed in figure 5 against a set of theoretical and experimental data from the literature, with close agreement for
${\it\vartheta}$
ranging from 0.3 to 0.9. Thanks to figures 4 and 5, one can conclude that model (2.12) reproduces well the values of the permeability both in isotropic and anisotropic cases. The homogenization technique is thus suitable from a microscopic point of view and for low-
$Re$
flows, despite the fact that the permeability components vary over a wide range of values, from
$O(10^{-4})$
to
$O(10^{-1})$
depending on the porosity, an indication of the fact that
$l^{2}$
is the proper scale for the permeability only at large porosities
${\it\vartheta}$
’s. This is confirmed by numerical solutions of (2.12), which yield
$\mathscr{K}_{ij}\rightarrow {\it\delta}_{ij}$
for
${\it\vartheta}\rightarrow 1$
.
Figure 4. Permeability versus porosity for regular arrays of spheres: homogenization theory (▫ (red)), empirical law of Kozeny–Carman (—— (blue)), theoretical results by Zick & Homsy (Reference Zick and Homsy1982) (▹ (blue)).
Figure 5. Components of the permeability tensor versus porosity for regular arrays of cylinders. The present results are represented by bullets, red for
$\mathscr{K}_{11}=\mathscr{K}_{22}$
, black for
$\mathscr{K}_{33}$
.
2.3
$Re=O(1)$
When
$Re$
is not small, (2.10) remains unaltered. We consider a steady solution of (2.10) through a Oseen approximation of the advective term, i.e. we assume that
$$\begin{eqnarray}u_{j}^{(0)}\frac{\partial u_{i}^{(0)}}{\partial x_{j}}\approx \overline{U}_{j}\frac{\partial u_{i}^{(0)}}{\partial x_{j}},\end{eqnarray}$$
where
$\overline{U}_{j}$
is the mean fluid velocity vector through the medium defined via a spatial average:
$$\begin{eqnarray}\overline{U}_{j}:=\frac{1}{V_{Tot}}\int _{V_{Tot}}\langle u_{j}^{(0)}\rangle \,\text{d}V,\end{eqnarray}$$
with
$V_{Tot}$
the whole macroscopic volume of the porous medium (fluid plus solid). We can also define the mean fluid velocity vector through the pores as
$\boldsymbol{v}^{p}=\overline{\boldsymbol{U}}/{\it\vartheta}$
and the pore Reynolds number,
$Re_{p}=U_{p}l/{\it\nu}$
, where
$U_{p}$
is the scale of the pore velocity
$\boldsymbol{v}^{p}$
.
The Oseen approximation has been proposed for cases with
$Re\gg 1$
by Gustafsson & Protas (Reference Gustafsson and Protas2013), with satisfactory results. We adopt it here in order to maintain a linear problem yielding (2.16). The choice of a global estimate of the velocity
$\overline{U}$
, defined via (2.22), is important to establish a two-way link between the microscopic and the macroscopic set of equations. Hence, we can write the solution of the new linearized equations as in (2.11) and search for
$\unicode[STIX]{x1D612}_{ij}$
and
$A_{j}$
which satisfy
$$\begin{eqnarray}\frac{\partial \unicode[STIX]{x1D612}_{ij}}{\partial x_{i}}=0,\quad -{\it\vartheta}Re_{p}{v^{p}}_{l}\frac{\partial \unicode[STIX]{x1D612}_{ij}}{\partial x_{l}}=\frac{\partial A_{j}}{\partial x_{i}}-\frac{\partial ^{2}\unicode[STIX]{x1D612}_{ij}}{\partial x_{g}^{2}}-{\it\delta}_{ij},\end{eqnarray}$$
with the same boundary conditions as those used in problem (2.12). For the definition of
$\boldsymbol{v}^{p}$
to make sense, we need to implement an iterative method, as we will explain in § 4.1. Finally, Darcy’s law has the same form as the canonical one (2.16), but the permeability tensor arises from iterating between the macroscopic equation (2.16) and the system (2.23) valid at the microscopic level. Thus, the permeability in this case is no longer a simple material property, but depends on the Reynolds number and the orientation of the velocity vector; it has been termed apparent permeability by Edwards et al. (Reference Edwards, Shapiro, Bar-Yoseph and Shapira1990). In the simple case in which the fluid has a single main direction of motion
$x_{1}$
, with constant pore velocity, i.e.
$\boldsymbol{v}^{p}=(1,0,0)$
, (2.23) becomes:
$$\begin{eqnarray}\frac{\partial \unicode[STIX]{x1D612}_{ij}}{\partial x_{i}}=0,\quad -{\it\vartheta}Re_{p}\frac{\partial \unicode[STIX]{x1D612}_{ij}}{\partial x_{1}}=\frac{\partial A_{j}}{\partial x_{i}}-\frac{\partial ^{2}\unicode[STIX]{x1D612}_{ij}}{\partial x_{g}^{2}}-{\it\delta}_{ij}.\end{eqnarray}$$
2.3.1 Effect of
$Re_{p}$
on the apparent permeability
The procedure briefly outlined in § 2.3 for calculating the apparent permeability at varying values of the pore Reynolds number is here validated against simulations of the Navier–Stokes equations for the flow through regular arrays of cylinders arranged in a square lattice. The reference simulations by Edwards et al. (Reference Edwards, Shapiro, Bar-Yoseph and Shapira1990) are two-dimensional (as such they furnish only the value of
$\mathscr{K}_{11}$
), and are carried out with a finite element solver for
$Re_{p}$
up to approximately 400 and porosities
${\it\vartheta}$
in the range 0.4–0.8. The comparison between our solutions of (2.24) and the reference values of
$\mathscr{K}_{11}$
is shown in figure 6, as a function of the Reynolds number and the porosity of the medium, together with our solutions for
$\mathscr{K}_{33}$
. The Oseen-based approximations of the first component of the permeability tensor match satisfactorily the reference values, providing confidence in the model developed, particularly for larger porosities.
Figure 6. Components of the apparent permeability as a function of the pore Reynolds number. The arrows point in the direction of increasing values of
${\it\vartheta}$
, from
${\it\vartheta}=0.4$
to
$0.8$
in intervals of
$0.1$
. Symbols in (a) refer to the results by Edwards et al. (Reference Edwards, Shapiro, Bar-Yoseph and Shapira1990).
3 Interface conditions
The theory described thus far is based on the spatial homogeneity of the porous medium and is not valid, in general, when homogeneity is broken by, e.g., the presence of an interface with a different medium. We must thus pay particular attention to the conditions to be imposed to transfer information from the pure fluid region (denoted in the following by
$F$
with the corresponding unknown fields
$\left.\cdot \right|_{F}$
) to the porous medium (denoted by
$P$
with the corresponding unknown fields
$\left.\cdot \right|_{P}$
) and vice versa: near the actual physical interface a more or less thin buffer layer exists where a matching must be enforced.
Several interface conditions are present in the literature and a recent review can be found in Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013). We have tested many of them comparatively; all conditions contain a certain degree of arbitrariness and require the a posteriori calibration of one or more parameters. We should also highlight the fact that a ‘double’ set of conditions is needed, one that links the fluid equations to Darcy’s law, and the other which takes the results computed in
$P$
and transforms them into suitable interface conditions for the Navier–Stokes solver. Among the different strategies that we have tested, we outline and compare three of them.
A. The first approach consists in solving Darcy’s equation in the porous region and imposing
$$\begin{eqnarray}\left.p^{(0)}\right|_{P}=\left.p\right|_{F}+\text{constant},\end{eqnarray}$$
at a fictitious interface positioned a small distance
${\it\delta}$
, which we call the penetration depth, below the physical interface, in order to transfer information from the fluid region to the porous region. The presence of the constant term has no influence on the solution of Darcy’s equation, since the pressure is always specified (both in
$P$
and
$F$
) up to an arbitrary constant; the constant is however maintained, at least formally, in (3.1) to highlight the fact that pressure may be discontinuous when crossing the interface. More specifically, a pressure jump has been predicted by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) and Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013) when the porous medium is anisotropic or in the presence of flow inertia (see also Sahraoui & Kaviany Reference Sahraoui and Kaviany1992), and it varies along the interface direction. On the other hand, continuity of pressure across the interface, up to the order of the pore size, was advocated by Ene & Sánchez-Palencia (Reference Ene and Sánchez-Palencia1975). We will come back to the pressure jump condition in § 4.6.
To have adequate conditions for the Navier–Stokes equation, we use
$$\begin{eqnarray}\left.u_{i}\right|_{F}=\left.\langle u_{i}^{(0)}\rangle \right|_{P},\end{eqnarray}$$
applied at the same fictitious interface as above. Here
$\left.u_{i}\right|_{F}$
is not averaged because it is the effective velocity, computed over a grid in which each point corresponds to a microscopic cell. Enforcing the continuity of both pressure and velocity below the real interface has been proposed and tested by Le Bars & Worster (Reference Le Bars and Worster2006), with the scope of accounting for inertia within the porous layer in the immediate vicinity of the fluid domain. Le Bars & Worster (Reference Le Bars and Worster2006) have shown that this is essentially equivalent to the interface condition of Beavers & Joseph (Reference Beavers and Joseph1967), which is applied at the physical interface, and have estimated the distance
${\it\delta}$
to vary as
$$\begin{eqnarray}{\it\delta}=c\sqrt{\frac{\mathscr{K}}{{\it\vartheta}}},\end{eqnarray}$$
where
$\mathscr{K}$
represents the permeability (in their case a scalar since the medium is isotropic) and
$c$
is a constant to be determined (
$c=1$
in their work). By enforcing the conditions at
${\it\delta}$
, we are sufficiently far from the interface for the periodicity assumption of homogenization theory to remain tenable.
B. The second approach (e.g. Kaviany Reference Kaviany1995) consists in modifying Darcy’s equation adding a viscous term through which we can impose continuity of the velocity and the tangential and normal components of the stress tensor over the interface. We note that Brinkman’s equation can be obtained from homogenization theory (Mei & Vernescu Reference Mei and Vernescu2010), adopting slightly different normalizations of the starting equations on the assumption that the fibres are sufficiently sparse. The resulting equation, in terms of the scales used in the multiple scale theory, is
$$\begin{eqnarray}\langle u_{i}^{(0)}\rangle =-\mathscr{K}_{ij}\frac{\partial p^{(0)}}{\partial x_{j}^{\prime }}+\mathscr{K}_{ij}{\it\epsilon}^{2}\frac{{\it\mu}_{e}}{{\it\mu}}\frac{\partial ^{2}\langle u_{j}^{(0)}\rangle }{\partial x_{k}^{\prime 2}},\end{eqnarray}$$
where
${\it\mu}_{e}$
is the (a priori unknown) effective viscosity of the homogenized medium. In order to match the results at the interface with the fluid region
$F$
, we should scale the pressure in Brinkman’s equation with
${\it\rho}U^{2}$
, so that the newly normalized equation reads:
$$\begin{eqnarray}\langle u_{i}^{(0)}\rangle =-\mathscr{K}_{ij}{\it\epsilon}^{2}Re_{L}\frac{\partial p}{\partial x_{j}}+\mathscr{K}_{ij}{\it\epsilon}^{2}\frac{{\it\mu}_{e}}{{\it\mu}}\frac{\partial ^{2}\langle u_{j}^{(0)}\rangle }{\partial x_{k}^{2}},\end{eqnarray}$$
with
$Re_{L}$
the macroscopic Reynolds number,
$Re_{L}=Re/{\it\epsilon}$
, with the primes
$^{\prime }$
now omitted from the name of the macroscale variable,
$x_{i}^{\prime }$
. We also omit the superscript
$(0)$
from the variable
$p$
in the porous domain. If
${\it\mu}_{e}$
were available we could impose the following conditions at the fluid–porous interface:
$$\begin{eqnarray}\left.({\it\sigma}_{ij}n_{j})n_{i}\right|_{F}=\left.({\it\sigma}_{ij}n_{j})n_{i}\right|_{P}\quad \text{and}\quad \left.({\it\sigma}_{ij}n_{j})t_{i}\right|_{F}=\left.({\it\sigma}_{ij}n_{j})t_{i}\right|_{P}\end{eqnarray}$$
to transfer data from the fluid domain to the porous medium
$P$
ruled by Brinkman’s equation, plus (3.2) to go from the porous domain to the fluid domain. In (3.6),
${\it\sigma}_{ij}$
is the stress tensor,
$n_{i}$
and
$t_{i}$
are, respectively, the unit vectors normal and tangent to the interface. Such conditions can be written in the present case as
$$\begin{eqnarray}\left.\left(\frac{\partial u_{1}}{\partial x_{3}}+\frac{\partial u_{3}}{\partial x_{1}}\right)\right|_{F}=\left.\frac{{\it\mu}_{e}}{{\it\mu}}\left(\frac{\partial \langle u_{1}^{(0)}\rangle }{\partial x_{3}}+\frac{\partial \langle u_{3}^{(0)}\rangle }{\partial x_{1}}\right)\right|_{P}\end{eqnarray}$$
and
$$\begin{eqnarray}\left.\left(-p+\frac{2}{Re_{L}}\frac{\partial u_{3}}{\partial x_{3}}\right)\right|_{F}=\left.\left(-p+\frac{{\it\mu}_{e}}{{\it\mu}}\frac{2}{Re_{L}}\frac{\partial \langle u_{3}^{(0)}\rangle }{\partial x_{3}}\right)\right|_{P}.\end{eqnarray}$$
We note that the stress from the pure fluid region
$F$
is transmitted to the homogenized region
$P$
, and is distributed to both the inclusions (the fibres) and the fluid contained within them.
C. The third interface condition tested consists in applying the same strategy as in A with
${\it\delta}=0$
, i.e. the matching is enforced precisely at the physical interface, but the permeability components are driven smoothly to a large absolute value near the interface, through a filter which forces the permeability to vary rapidly over a distance of the order of
${\it\epsilon}$
. This strategy has been proposed by Chandesris & Jamet (Reference Chandesris and Jamet2008) and Jamet & Chandesris (Reference Jamet and Chandesris2009); in the present case the filter is taken directly from three-dimensional numerical simulation results with account of all individual fibres.
Figure 7. View of the macroscopic domain with the boundary conditions imposed when solving Darcy’s equation in the
$P$
-region.
4 First validation case: laminar flow in a cavity
We first solve for the flow in a composite domain, in which there is both a pure fluid region and a porous region formed by densely packed rigid fibres: the choice of considering a lid-driven cavity, filled up to a certain height (each fibre is in length
$0.33$
times the side length of the square cavity) by the porous medium, as sketched in figure 7, allows us to analyse the orthotropic character of the medium. In fact, because of the non-negligible vertical velocity inside the porous region, we are able to quantify how much
$K_{33}$
affects the fluid flow and we can impose and compare different interface conditions. On the other hand, this case is more difficult than the classical pressure-driven channel flow with porous walls studied by many authors before.
The method used is classified as a domain decomposition method; it is an iterative method based on splitting interface conditions, which play the role of transmission conditions (cf. Quarteroni & Valli Reference Quarteroni and Valli1999), between the
$F$
and
$P$
regions. Searching for a steady solution we formalize the numerical procedure with the following steps:
-
(i) We solve for the two-dimensional NSE in the
$F$
-region, expressed in dimensionless form as: (4.1)
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}}{\partial x_{i}}=0, & \displaystyle\end{eqnarray}$$
(4.2)with initial interface condition
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}}{\partial t}+u_{j}\frac{\partial u_{i}}{\partial x_{j}}=-\frac{\partial p}{\partial x_{i}}+\frac{1}{Re_{L}}{\rm\nabla}^{2}u_{i}, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}=0$
;
$Re_{L}$
is based on the characteristic dimension of the cavity (
$L$
) which can be identified with the macroscopic length of the homogenization technique, so that
$Re_{L}=Re/{\it\epsilon}$
. The solver developed and validated in Zampogna (Reference Zampogna2016) is based on the fractional step method, using fourth-order finite differences for the space discretization and a Crank–Nicolson scheme for the time advancement (Perot Reference Perot1993).
-
(ii) We take the divergence of (2.16) and use (2.18) to have in the
$P$
-region: (4.3)plus homogeneous Neumann conditions on the solid walls. At the interface the pressure is imposed (cf. (3.1)). If the Brinkman’s model is considered instead, conditions (3.7) and (3.8) are used at the interface, together with no slip along the walls, to solve for the equations in the porous domain and find the pressure
$$\begin{eqnarray}\displaystyle \mathscr{K}_{11}\frac{\partial ^{2}p}{\partial x_{1}^{2}}+\mathscr{K}_{33}\frac{\partial ^{2}p}{\partial x_{3}^{2}}=0, & & \displaystyle\end{eqnarray}$$
$\left.p\right|_{P}$
and the velocity components (
$\langle u_{1}^{(0)}\rangle ,\langle u_{3}^{(0)}\rangle$
).
-
(iii) When using Darcy’s model we still need to compute the velocity in the porous medium, using (3.5) with
${\it\mu}_{e}$
set to zero. For either models the final step of the iteration procedure consists in enforcing condition (3.2) at the interface for the subsequent solution in the
$F$
-region.
4.1 Double interactions for the
$Re=O(1)$
case
The scheme of the
$Re=O({\it\varepsilon})$
case remains unchanged also when
$Re=O(1)$
, except for one point. In this case the permeability tensor cannot be computed before starting the macroscopic simulation but must be computed each time that Darcy’s velocity is updated because
$\unicode[STIX]{x1D612}_{ij}$
depends directly on the magnitude and direction of the mean velocity inside the porous medium (cf. 2.23). A guessed value of
$\unicode[STIX]{x1D612}_{ij}$
is thus chosen to initialize the iterative procedure and a reasonable choice is to use the permeability of the Stokes flow case. After each solution in the
$P$
-region, the mean velocity components are computed with (2.22) and used to estimate the permeability for the successive time step.
4.2 Direct numerical simulations
In order to compare the different interface conditions and assess which of them fits the problem analysed in this work, three-dimensional direct numerical simulations (DNS) of the incompressible NSE in the real geometry, accounting for all fibres (50 for the simulations presented here, so that
${\it\epsilon}=0.02$
), are performed with OpenFOAM. Also in this case all the operators (divergence, gradient and Laplacian) are discretized using the finite volume method, with the Gauss linear scheme. A Crank–Nicolson scheme is used to march in time.
Figure 8. Spatial convergence study for the DNS in the case of
$Re_{L}=100$
. The different grids are built with the software snappyHexMesh (http://www.openfoam.org/version2.3.0/snappyHexMesh.php).
In the transverse direction (
$x_{2}$
) the size of the domain is taken equal to
${\it\epsilon}$
(in dimensionless terms) and periodic boundary conditions are enforced. The numerical mesh needed to have grid-converged solutions is formed by over 8 million cells; figure 8 provides an example of a convergence test for a measure (the PDF of the permeability components) introduced further in (§ 4.3).
An example of a DNS result (for
$Re=100$
) is shown in figure 9, where a slice of the steady converged solution at
$x_{3}=0.25$
(thus within the porous medium) is shown, through isocontours of the first velocity component, near the solid boundaries and around
$x_{1}=0.5$
. The figure also shows a cut through the numerical grid.
Figure 9. Zoom of the macroscopic domain of the DNS inside the porous medium with
${\it\vartheta}=0.8$
. Slice of the cavity at
$x_{3}=0.25$
. Only few cylinders (out of 50) are shown inside and at the extremities of the cavity. The colours and the white curves in (a) represent the
$x_{1}$
-component of the velocity field and its isocontours. In (b) a cut through the grid is shown.
With the results from the DNS we can:
-
(i) evaluate the components of the permeability tensor and verify the results of homogenization theory;
-
(ii) test the different interface conditions and, for example, calibrate the constant
$c$
in the definition of
${\it\delta}$
(3.3); -
(iii) validate the unknown fields at a macroscopic level, verifying the appropriateness of the homogenized equations.
To compare the DNS results with the two-dimensional simulations described in the previous section we must extract in a proper way the averaged values of the fields: the whole domain is thus decomposed into elementary cubic cells of size
${\it\epsilon}^{3}$
and all fields are averaged over each microcell using the definition (2.15); in this way two-dimensional fields are obtained from the DNS results and, in particular, we have a sampling of the solution in the porous medium (in the case of a layer of
$N$
filaments of height
$0.33$
we have a Cartesian grid, with velocity components and pressure centred in the cell, composed by
$0.33\,N^{2}$
elements in
$P$
). Using the grid built, the pressure gradient is computed with fourth-order finite differences. In the internal region of the porous medium, assuming that Darcy’s approximation is valid (i.e.
$\langle \boldsymbol{u}^{DNS}\rangle =\langle \boldsymbol{u}^{(0)}\rangle$
and
$\langle \,p^{DNS}\rangle =p^{(0)}$
),
$\mathscr{K}_{ij}$
is computed using (2.16), using the same scalings which have led to (3.5). Finally, we need to account for the influence of the boundary layers at the solid walls and at the interface, and for the presence of isolated regions in which the pressure gradient becomes negligibly small; thus, the values of
$\mathscr{K}_{ij}^{DNS}$
(superscript
$DNS$
is used to indicate that the value has been estimated on the basis of the full simulations) are evaluated on the basis of their probability density functions (PDF). To compute the PDF we divide the values of
$\mathscr{K}_{ii}$
in a number of intervals such that each one of them has an extent equal to 10 % of the value of
$\mathscr{K}_{ii}$
for which we find the maximum of the PDF.
4.3 The components of the permeability tensor
The complete simulations have been performed for Reynolds numbers
$Re_{L}$
varying from
$1$
to
$1000$
, for fixed value of
${\it\epsilon}=0.02$
and
${\it\vartheta}=0.80$
. The a posteriori treatment of the DNS results within the porous medium permits to extract the value of the permeability components over each elementary cubic cell, yielding the probability density functions displayed in figure 10.
Figure 10. Probability density function of
$\mathscr{K}_{11}$
and
$\mathscr{K}_{33}$
for different values of
$Re_{L}$
. The apparent permeability is computed a posteriori applying Darcy’s law to the averaged solution of the DNS. (a)
$Re_{L}=10$
, (b) 100, (c) 1000.
It is immediately apparent that the PDF pinpoints sharply the values of
$\mathscr{K}_{11}$
and
$\mathscr{K}_{33}$
for
$Re_{L}$
up to
$100$
, whereas for larger
$Re_{L}$
’s the two components of the permeability display a broader variability. The peak values are nonetheless well defined and are reported in table 1, against corresponding results from the homogenization theory. The latter have been computed including inertial terms in the equations, with the same
$U_{p}$
found in the direct simulations of the full geometry. The values of the permeability components from the DNS are approximately 40 % larger than those calculated from the system (2.12), (2.13); in absolute terms, the difference is of order
${\it\epsilon}$
, consistent with the theory. We will see (figure 21) that something similar occurs also in another flow case, under turbulent conditions. In the narrow range of pore Reynolds number considered here, the full simulations do not allow us to extract a clear trend of how the permeability varies, whereas the theory indicates a slow decrease of both components with
$Re_{p}$
.
Table 1. Values of the permeability components from the DNS and homogenization theory.
$U_{p}$
represents the dimensionless mean pore velocity in the whole porous medium. Because of continuity, the mean pore velocity is oriented only along
$x_{1}$
.
$d_{1}$
and
$d_{3}$
are the distances from the interface starting from which
$\mathscr{K}_{11}$
and
$\mathscr{K}_{33}$
can be assumed constant when condition C is applied (cf. figure 16).
The disagreement between theory and direct simulations appears to be relatively mild, when seen in the scale of figure 6, but nonetheless deserves to be clarified, since it occurs also for flow cases (
$Re_{L}=1$
) in which the motion of the fluid through the pores is clearly ruled by Stokes’ equation. We believe that, on the one hand, the flow geometry plays a role here, since the two sidewalls constrain the fluid, forcing it to recirculate; a much better agreement would have been obtained had we considered a pressure-driven channel flow, as done by many authors (Le Bars & Worster Reference Le Bars and Worster2006; Hill & Straughan Reference Hill and Straughan2008; Battiato Reference Battiato2012). On the other hand, we would expect a closer match between the DNS results and the theory had we pushed the asymptotic development to third order, modifying Darcy’s equation with a Forchheimer term, as described by Mei & Vernescu (Reference Mei and Vernescu2010). This is indirectly supported by one direct simulation at
$Re_{L}=100$
, conducted with 200 fibres (instead of 50) – for which the small parameter
${\it\epsilon}$
of the expansion is small (and equal to
$0.005$
), for the same porosity coefficient
${\it\vartheta}=0.80$
– which yields values of the apparent permeability closer to the theoretical ones in table 1.
In the following, the results from the three-dimensional numerical simulations will be compared to results from the coupled two-dimensional Navier–Stokes/Darcy equation, where in the latter model the values used for the apparent permeability are those of the direct numerical simulations.
4.4 The macroscopic fields
A representative comparison between the DNS results and those from the model system is displayed in figures 11 and 12 in terms of the streamlines of the flow and the pressure field, respectively. In the model system we have enforced interface conditions A (cf. § 4), with
${\it\delta}={\it\epsilon}/5$
, i.e. we have used the penetration depth
${\it\delta}$
that, as we will show, yields the best agreement between the model and the direct simulations. The results are very close to one another and, in particular, the complex behaviour of the flow within the porous domain appears to be correctly captured by Darcy’s equation (with a tensorial permeability). Also the pressure field of the homogenized model matches that of the DNS with the constant term of (3.1) appropriately chosen.
Figure 11. Streamlines in the composite domain. (a,b) Represent the solution for
$Re_{L}=100$
; (c,d) display results for
$Re_{L}=1000$
. In all cases the porosity is
${\it\vartheta}=0.80$
. (a,c) Show solutions of the homogenized model, (b,d) are the DNS results. Representative velocity profiles will be later displayed along the vertical line
$x_{1}=0.5$
.
Figure 12. Isocontours of the pressure field in the composite domain. Panels (a,b) represent the solution for
$Re_{L}=100$
; (c,d) display results for
$Re_{L}=1000$
. In all cases the porosity is
${\it\vartheta}=0.80$
. Panels (a,c) are the solutions of the homogenized model, (b,d) are the DNS results.
It is significant that the model problem requires a CPU time of the order of 30 minutes using one 2.0 GHz Intel®Xeon® processor and less than 2 Gb of RAM to reach convergence, whereas each DNS, on the same computer, needs a CPU time of over 500 h and more than 8 Gb of RAM, running in parallel on 16 processors. It is also noticeable that for different interface conditions the results are qualitatively similar, and all the model solutions look ‘correct’, with minor differences noticeable (corresponding images are not shown for reasons of space). This statement, however, justifies a closer look at the interface.
4.5 Focus on the interface
To assess the quality of the interface conditions we focus on
$x_{1}=0.5$
. The first condition (A of § 3) yields the results displayed in figure 13, for
$Re_{L}=100$
. A value of
${\it\delta}$
between
$0$
and
${\it\epsilon}/5$
provides a good match with the DNS results in the
$F$
-region and at the interface, whereas increasing
${\it\delta}$
above
${\it\epsilon}/2$
yields values of both velocity components near the interface which progressively overestimate the ‘true’ values. Within the porous domain (
$x_{3}\leqslant 0.25$
), the difference between the solutions computed for varying
${\it\delta}$
’s and the DNS is negligible.
Figure 13. View of the macroscopic profile of
$u_{1}$
and
$u_{3}$
(lower panels) for
$x_{1}=0.5$
(middle of the cavity) for the DNS and four choices of
${\it\delta}$
(
$Re_{L}=100$
). Two zooms are highlighted at the interface. Interface strategy: A.
Figure 14. View of the macroscopic profile of
$u_{1}$
and
$u_{3}$
(lower panels) for
$x_{1}=0.5$
for the DNS and four choices of
${\it\delta}$
(
$Re_{L}=1000$
). Interface strategy: A.
When
$Re_{L}$
is equal to 1000, an excellent agreement with the DNS results is achieved for
${\it\delta}={\it\epsilon}/5$
, as shown in figure 14. It is important to note that a poor choice of
${\it\delta}$
(e.g.
${\it\delta}\geqslant {\it\epsilon}$
) degrades the solution also in the fluid region
$F$
, as figure 14 clearly demonstrates. As far as interface condition A is concerned, we observe an increase of the penetration depth,
${\it\delta}$
, with
$Re_{p}$
, as one would intuitively expect.
We now turn to interface condition B. Brinkman’s equation has been discretized and directly matched to the Navier–Stokes results in the
$F$
-region; the velocity profiles obtained at
$x_{1}=0.5$
are compared to the reference profiles in figure 15, for
$Re_{L}=100$
and different values of the effective viscosity. As
${\it\mu}_{e}$
approaches zero, the Darcy limit is recovered; the upper value of
${\it\mu}_{e}/{\it\mu}$
considered is
$30$
, close to that determined experimentally by Givler & Altobelli (Reference Givler and Altobelli1994). Furthermore, even the theoretical value,
${\it\mu}_{e}/{\it\mu}=(1-{\it\vartheta})^{-1}=5$
, deduced on the basis of volume averaging theory by Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995), has been tested. The immediate observation is that a value of the effective viscosity, which is either too large (30) or too small (0.3), degrades significantly the solution in the
$F$
-region, as a direct consequence of the poor estimation of the velocity in the proximity of the interface. Conversely, when
${\it\mu}_{e}$
is a few times larger than the dynamic viscosity of the fluid, i.e. it is close to the theoretical approximation, the solution appears to behave better when compared to the DNS, although close inspection of the fields within the porous region (not shown) show differences between the model solution and the ‘exact’ result which are larger than those found with interface strategy A. This negative aspect of Brinkman’s model, coupled to the fact that the effective viscosity is not available a priori, is somehow tempered by the fact that Brinkman’s equation has been used here in a system for which the fibres are not sufficiently sparse, as prescribed by Mei & Vernescu (Reference Mei and Vernescu2010). A similar conclusion on the poor performance of Brinkman’s model in handling porous–fluid interfaces has also been reported by Sahraoui & Kaviany (Reference Sahraoui and Kaviany1992), Kaviany (Reference Kaviany1995) and James & Davis (Reference James and Davis2001).
Results from the last interface condition (C) are displayed in figures 16 and 17; here the permeability varies as the macroscopic interface (at
$x_{3}^{ITF}=0.33$
) is approached. Whereas this is physically reasonable, it is unknown how the permeability should be taken to vary. In this specific case we use the DNS results to define the shape of the filter to apply on
$\mathscr{K}$
. The permeability inside the porous medium is recovered using Darcy’s law, as explained in the previous section. Here, we need a value of permeability for each
$x_{3}$
in the porous zone. It is sampled for each
$x_{1}$
in the cavity and then a mean value along the
$x_{1}$
direction is extrapolated (solid and dotted lines of figure 16). It is clear that up to a distance
$d=O({\it\epsilon})$
from the interface, the value of
$\mathscr{K}_{ii}$
is constant; then it increases. In particular we observe that for the case
$Re_{L}=100$
shown in figure 16 the component
$\mathscr{K}_{11}$
has
$d_{1}=(5/2){\it\epsilon}$
, while for
$\mathscr{K}_{33}$
it is
$d_{3}=(3/2){\it\epsilon}$
. This seems to indicate that
$\mathscr{K}_{ii}^{-1}$
is proportional to
$d_{i}$
; in table 1 this trend is further explored for varying
$Re_{L}$
’s. The distance
$d_{i}$
(
$i=1,3$
) defines an interfacial layer whose thickness decreases slowly with the increase of the Reynolds number
$Re_{L}$
. In this it differs from the protrusion height
${\it\delta}$
of approach A. The agreement at the interface between the homogenized solution and the solution of the DNS is however not very good (cf. figure 17) and this is likely due to the fact that near the interface Darcy’s law is inappropriate and a model which accounts for inertia is needed.
Figure 15. View of the macroscopic profile of
$u_{1}$
and
$u_{3}$
(lower panels) for
$x_{1}=0.5$
and four choices of
${\it\mu}_{e}/{\it\mu}$
(
$Re_{L}=100$
). The images on the right highlight the degradation of the solution with increasing
${\it\mu}_{e}$
, in particular close to the interface. Interface strategy: B.
Figure 16. Inverse of the filtered permeability for the cases
$Re_{L}=100$
. The two lines (dotted and solid) arise from the DNS and are
$x_{1}$
-averaged values;
$d_{1}$
and
$d_{3}$
are the lengths over which the components of permeability vary near the interface positioned in
$x_{3}=0.33$
.
Figure 17. Macroscopic profile of
$u_{1}$
and
$u_{3}$
(lower panels) for
$x_{1}=0.5$
computed using the filtered permeability shown in figure 16 (
$Re_{L}=100$
) and compared to the case in which the permeability components are maintained constant throughout the porous zone. Two zooms over the interfacial zone are highlighted on the right. Interface strategy: C.
4.6 The pressure jump at the interface
The availability of the DNS renders it possible to assess the presence of, and quantify, the pressure jump across the physical interface. Arguments for the existence of a pressure jump have been put forward, using homogenization theory coupled with a boundary layer analysis near the interface, by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) for the Stokes flow over an anisotropic porous bed. Similarly, a pressure slip at the interface has been reported by Sahraoui & Kaviany (Reference Sahraoui and Kaviany1992) in the presence of flow inertia; analytical arguments to quantify the pressure slip are given by Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2015) in the case of forced infiltration. The condition of Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) in dimensionless form, reads:
$$\begin{eqnarray}-[p]=p_{P}(x_{1},ITF^{-})-p_{F}(x_{1},ITF^{+})=\frac{C}{Re_{L}}\frac{\partial u_{1}}{\partial x_{3}}(x_{1},ITF),\end{eqnarray}$$
with the first term above representing the pressure jump across the interface, positioned at
$x_{3}=ITF$
and
$C$
a constant which depends only on the porous bed geometry. We have evaluated precisely this jump on the basis of the DNS, by averaging the microscopic data, both in
$F$
and
$P$
, over cubic unit cells of size
${\it\epsilon}^{3}$
. Thus, the pressure in the pure fluid is evaluated at
$x_{3}=ITF^{+}=ITF+({\it\epsilon}/2)$
, and the average pressure in
$P$
is evaluated at
$x_{3}=ITF^{-}=ITF-({\it\epsilon}/2)$
. The left panel of figure 18 displays
$|Re_{L}[p]|$
along the interface at the three values of
$Re_{L}$
considered; the right panel of the same figure reports a shear parameter defined as:
$$\begin{eqnarray}S=\frac{\partial u_{1}}{\partial x_{3}}(x_{1},ITF)+\frac{\partial u_{3}}{\partial x_{1}}(x_{1},ITF).\end{eqnarray}$$
In our simulations
$(\partial u_{3}/\partial x_{1})(x_{1},ITF)$
at the interface is typically two orders of magnitude smaller than the
$(\partial u_{1}/\partial x_{3})(x_{1},ITF)$
term, which thus dominates the parameter
$S$
. Comparison of the two panels makes it clear that, in the presence of inertia, a simple constant
$C$
cannot be found to satisfy the pressure jump condition, which is thus unsuited for use in the present case. Henceforth, although the pressure jump is present, condition (4.4) cannot be readily enforced.
Figure 18. (a) Pressure jump across the interface at three values of the Reynolds number. (b) Shear parameters
$S$
at the interface. Both panels are based on DNS results.
5 Second validation case: turbulent canopy flow
In this section we want to demonstrate that the model formed by (2.16), (2.24) plus the equations in the pure fluid domain and the interface conditions are adequate also to simulate a turbulent canopy flow. This kind of flow has been widely studied in the past since it develops very often in nature (a recent review is provided by Nepf (Reference Nepf2012)). Much literature (theoretical and experimental) is present to understand the mechanisms which induce the creation of a quasi-parallel mixing layer and the formation of honami and monami waves over canopies. Raupach, Finnigan & Brunet (Reference Raupach, Finnigan and Brunet1996) and Finnigan (Reference Finnigan2000) have classified the different structures which characterize the flow, on the basis of the geometrical properties of the canopy (e.g. its density) and the parameters of the mixing layer.
Figure 19. Setting of the experiments of Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004). Water flows over a layer of rigid cylindrical fibres. The whole system is slightly inclined at an angle
${\it\alpha}$
so that the fluid motion is driven by gravity.
In order to understand if the closure adopted for the
$Re=O(1)$
case in (2.24) makes sense we have used the experimental data of Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004) where a fully developed turbulent flow of water over a submerged layer of vertical rigid fibres is studied. The setting of the problem is sketched in figure 19; the porosity
${\it\vartheta}$
in the experiments ranges from
$0.96$
to
$0.99$
and
${\it\varepsilon}$
is approximately 0.2 (defined as the ratio between the diameter of the fibres and the height of the water column). The measurements of the mean velocity inside the canopy and the driving term are available so that the permeability
$\mathscr{K}_{11}$
can be estimated a posteriori using Darcy’s law which, in this case, can be written in dimensional variables as
$$\begin{eqnarray}\langle u_{1}^{(0)}\rangle =\mathscr{K}_{11}\frac{g\sin {\it\alpha}}{{\it\nu}},\end{eqnarray}$$
where
${\it\nu}=0.94\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
is deduced using
$Re_{ml}$
in table 1 of Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004) and
${\it\alpha}$
is the slope of the inclined system. The use of (5.1) in the lower zone of the canopy, termed ‘wake zone’, is justified on physical grounds since, as noted by Ghisalberti & Nepf (Reference Ghisalberti and Nepf2009), such a region is ‘governed by a simple balance of drag and hydraulic gradient, much like classical porous medium flow’.
Figures 20 and 21 display results for the microscopic permeability field
$\unicode[STIX]{x1D612}_{11}$
within the elementary volume, and the averaged values
$\mathscr{K}_{11}$
for variable pore velocity. The grey band represents the range of permeability
$\mathscr{K}_{11}$
computed using the homogenized model (2.24) for the corresponding porosities of each experiment. The experimental values are in reasonably good agreement with the homogenized values; discrepancies are of the same order as those shown in figure 6 and table 1.
Figure 20. Horizontal component of the permeability
$\unicode[STIX]{x1D612}_{11}$
within a unit cell, for varying
${\it\vartheta}Re_{p}=0,50,500,1500$
. Increasing
$Re_{p}$
, the value of
$\unicode[STIX]{x1D612}_{11}$
is reduced and the distribution initially loses symmetry with respect to the vertical mid-plane through the cylinder, to eventually regain it because of periodicity for
${\it\vartheta}Re_{p}$
above
$500$
. Here
${\it\vartheta}=0.975$
.
Figure 21. Horizontal component of the permeability
$\mathscr{K}_{11}$
versus
${\it\vartheta}Re_{p}$
of (2.24) from the measurements by Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004). Since the permeability is deduced from experimental data, the vertical and horizontal bars represent the uncertainty of the measurements which is approximately 10 % for each value. The grey band represents the range of permeabilities calculated by the homogenized model for finite Reynolds numbers for
$0.96\leqslant {\it\vartheta}\leqslant 0.99$
.
We now consider the problem from a macroscopic point of view. Observing the sketch of the velocity profile in figure 19 it seems reasonable to impose continuity of the streamwise velocity at a distance
${\it\delta}$
(the penetration depth) below the top of the canopy. This is the fictitious interface and (5.1) is thus applicable for
$x_{3}\in [0,h-{\it\delta}]$
. Since we are dealing with a fully developed turbulent flow, in the region above the fictitious interface (
$x_{3}\in [h-{\it\delta},x_{3}^{\infty }]$
), after using the Reynolds decomposition for the velocity (i.e.
$u_{i}=U_{i}+u_{i}^{\prime }$
) plus temporal averaging (
$\overline{\cdot }$
), the streamwise momentum equation reads:
$$\begin{eqnarray}\frac{\partial \overline{u_{1}^{\prime }u_{3}^{\prime }}}{\partial x_{3}}=g\sin {\it\alpha}.\end{eqnarray}$$
Thus, the turbulent stress balances the hydraulic gradient in the upper zone (termed the ‘exchange zone’), whereas the viscous term can be neglected. Assuming a constant mixing length model, such as that proposed by Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004), the Reynolds stress can be written as
$$\begin{eqnarray}-\overline{u_{1}^{\prime }u_{3}^{\prime }}=\left(cl_{{\it\delta}}^{\infty }\frac{\partial U_{1}}{\partial x_{3}}\right)^{2},\end{eqnarray}$$
where
$l_{{\it\delta}}^{\infty }$
is shown in figure 19,
$c$
is a constant and their product is a ‘mixing length’. Substituting (5.3) into (5.2) and solving for the mean flow it is found:
$$\begin{eqnarray}U_{1}(x_{3})=U_{1}(h-{\it\delta})+\frac{2}{3}\frac{\sqrt{g\sin {\it\alpha}}}{cl_{{\it\delta}}^{\infty }}[(x_{3}^{\infty }-h+{\it\delta})^{3/2}-(x_{3}^{\infty }-x_{3})^{3/2}].\end{eqnarray}$$
The condition at the fictitious interface yields
$U_{1}(h-{\it\delta})=\langle u_{i}^{(0)}\rangle (h-{\it\delta})=\mathscr{K}_{11}(g\sin {\it\alpha}/{\it\nu})$
. As we can observe from figure 22 the proposed solution fits the experiments with the unique value of
$c=0.086$
. Summing up, we can make the following two points:
-
(i) The values of the permeability deduced from the experiments slightly overestimate those evaluated via homogenization theory (cf. figure 21), and this is possibly due to the fact that in the experiments by Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004) the parameter
${\it\epsilon}$
is approximately
$0.2$
, a value which is not much smaller than one. -
(ii) The agreement between the experimental and theoretical velocity profiles in figure 22 confirms that the interface condition A of Le Bars & Worster (Reference Le Bars and Worster2006) is suitable in this turbulent case. Moreover, for the four experiments analysed, the order of magnitude of
${\it\delta}$
estimated from (3.3) matches that used here. However, if
$c$
in (3.3) were a simple constant, considering the trend of
$\mathscr{K}_{11}$
in figure 21,
${\it\delta}$
would decrease for increasing values of
$Re_{p}$
, contrary to observations (cf. the inset in figure 22). Thus, the ‘constant’
$c$
in (3.3) should be rendered an increasing function of
$Re_{p}$
for
${\it\delta}$
to increase monotonically with the increase of the pore velocity (at each given value of
${\it\theta}$
).
Figure 22. Analytical velocity profiles, in dimensional variables, against experimental results by Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004, Reference Ghisalberti and Nepf2006, Reference Ghisalberti and Nepf2009) for turbulent flow over a layer of rigid fibres. The symbols represent the velocity measurements (with uncertainty of 5 %) for cases B, H, I and J of table
$1$
of Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004). The analytical solution, (5.4), is represented by the solid lines. Two horizontal grey lines, which appear overlapped in the scale of the figure, represent the physical interface
$h$
(in dimensional terms
$h=0.139~\text{m}$
for case B and
$h=0.138~\text{m}$
for the other cases). In the inset the values of
${\it\delta}$
, which define the fictitious interface at
$x_{3}=h-{\it\delta}$
, are shown for the corresponding cases. The dashed line is drawn to guide the eye. In the four experiments we have used the following values for
$l_{{\it\delta}}^{\infty }$
: ✩: 0.202 m ▫: 0.245 m ○: 0.279 m ▹: 0.286 m. Furthermore, the microscale
$l$
(used to normalize
${\it\delta}$
in the inset) is equal to 5.06 cm in experiment B (with
${\it\vartheta}=0.99$
) and 2.83 cm in the other cases (for which
${\it\vartheta}=0.96$
).
6 Conclusions
A homogenization approach has been used to study the flow over and through a orthotropic porous medium. The work is developed essentially from two points of view: the microscopic and the macroscopic ones. The microscopic equations provide a characterization of the particular geometry of the porous medium, resulting in an estimate for the components of the permeability tensor; only two diagonal components of this second-order tensor are present in the particular medium considered since the terms relative to the plane normal to the fibres’ axes are equal. This result can be easily confirmed analytically (Milton Reference Milton2004). Since we are interested in the macroscopic field, the permeability tensor is averaged over the elementary cell and the values
$\mathscr{K}_{11}$
and
$\mathscr{K}_{33}$
are compared with values from the literature, obtaining good agreement. A new correlation for the (scalar) permeability
$\mathscr{K}$
of an isotropic medium is also provided as a function of the porosity
${\it\vartheta}$
of the medium, correcting the classical relation by Kozeny–Carman.
In the presence of inertia within the porous region, an approach to estimate the permeability is proposed, based on the Oseen linearization of the equations in the unit cell, plus an iterative approach which alternates between the microscopic and the macroscopic regions. With this technique, the apparent permeability from the theory approaches the values found in experiments and DNS, with the correct trend as a function of the pore Reynolds number. The theoretical values do not precisely match those available from experiments or DNS because the fibres considered here are neither sufficiently slender nor packed densely enough for the first-order equations of the multiple scale approximation to be fully adequate. To improve this, one could possibly push the multiple scale theory up to higher order.
The coupling of the flow problem within the porous medium with that in a pure fluid region neighbouring it is particularly important and different interface conditions have been evaluated. All of them contain some degree of arbitrariness and parameters which must be tuned. The main conclusion here is that Darcy’s law provides a good approximation of the velocity inside the porous medium, also when inertia in the flow within the
$P$
-region is non-negligible. The interface condition of Le Bars & Worster (Reference Le Bars and Worster2006), which is an extension of the one by Beavers & Joseph (Reference Beavers and Joseph1967), is found to be the most suitable, both in the case of laminar and turbulent flow. The penetration depth
${\it\delta}$
which needs to be imposed satisfies the order estimate of Le Bars & Worster (Reference Le Bars and Worster2006), for fixed pore Reynolds number. When
${\it\vartheta}Re_{p}$
increases also
${\it\delta}$
increases, and this should be accounted for when selecting
${\it\delta}$
in the model.
Conversely, the use of Brinkman’s correction to Darcy’s law should be avoided when the porosity is not large since, on the one hand, the value of the effective viscosity
${\it\mu}_{e}$
is uncertain and, on the other, results with this model appear to systematically overestimate the velocity near the interface between the fluid and the porous domain.
Finally, also a condition which employs permeability components which vary as the interface is approached do not appear adequate. The distance over which the permeability approaches the infinite value of the pure fluid region has been estimated to be of order
${\it\epsilon}$
on the basis of DNS, and depends mildly on Reynolds number, for the range considered which spans three orders of magnitude. The use of permeability components which change smoothly (even according to DNS results) is, however, not sufficient in a coupled NSE–Darcy system to represent inertial effects within a small porous layer adjacent to the pure fluid region, and this suggests the inclusion of a Forchheimer term in a thin-layer interface model (Firdaouss et al.
Reference Firdaouss, Guermond and Le Quéré1997).
Future work will focus on the flow of fluid over and through a bundle of flexible fibres.
Acknowledgements
The authors acknowledge the partial financial support provided by the EU to Wolf Dynamics s.r.l. through the PelSkin project (grant no. ACP2-GA-2013-334954-PEL-SKIN) and a grant from the Italian Ministry of Education through the PRIN 2012 initiative, project no. D38C1300061000.
