Discrimination between Pore and Throat Resistances against Single-Phase Flow in Porous Media

Abstract

This study investigates the critical agents that cause non-Darrian flow in porous media. Four porous media different in morphology but similar in topology were studied numerically. By varying the throat diameters, the distinct roles of pores and throats in total dissipation were investigated using direct numerical simulation. Forchheimer model was selected to analyze the non-Darcian flow. In our simplified geometry, the ratio $\frac{{\mathrm{K}}_{\mathrm{app}}}{{\mathrm{K}}_{\mathrm{D}}}$ can best be correlated by non-Darcy effect (E). Total dissipation is directly related to the porous medium resistance against fluid flow. The energy dissipated in pores and throats was calculated by summing the dissipation in each computational segment. Pores are more prone to disobey the Darcy model than throats due to irregularity in fluid flow, and they are introduced as the cause of Darcy-model cessation. By increasing the pore-to-throat ratio, the non-Darcian flow in the pores begins sooner. The results show that the energy dissipation due to eddies is negligible. The dissipation in pores and throats was simulated through separate power-law equations, and their exponents were also extracted. The exponent for the pore body is equal to two when the viscous forces are dominant, and it increases by increasing the inertia force. The dissipation due to pore bodies is more apparent when the size of pore and throats are of the same order of magnitude. The relative losses of pore body increase as the velocity increases, in contrast to throats.

1. Introduction

High-velocity fluid flow in porous media occurs in many situations, including near-wellbore gas flow [1], fractures [2], filtration, and tight screens of cryogenic propellant tanks. Among the models introduced in the past decades to simulate fluid flow in porous media at high velocities, the quadratic Forchheimer equation has attained much more attention [3]. The application of the equation is verified empirically by measuring the global pressure gradient and flow rate [4,5,6]. However, these experiments cannot capture microscopic details [7,8]. Therefore, the origin of non-Darcian flow remains unknown in such studies unless more sophisticated methods are implemented [9,10].

Several studies have asserted the central role of fluid patterns in microscale in the characterization of transport properties of porous media. These studies can be categorized into three main groups: field-scale observation methods that benefit primarily from fluid flow imaging in porous media [11,12,13], pure numerical studies [14,15], and machine-learning modeling [16,17,18]. Many publications have emphasized the role of eddies and inertial core flow in the changing flow from Darcy to non-Darcy conditions [19,20,21,22,23,24]. Numerical investigations at pore-scale level approved that the eddy formation can rise to non-Darcian regime and increases the pressure drop [19]. Microscopic heterogeneity raises inertial effects, disturbs the flow, and distorts velocity and pressure fields [23].

Porous media commonly can be divided into the void and solid phases [25]. Recognizing the void space as a uniform phase, while having the advantage of generality creates many computational complications. With the current state of computational facilities, modeling single and multiphase flow in porous media is a concern [26]. Discretizing the void space into pores and throats simplifies the problem, as the transporting elements can be analyzed separately. Furthermore, such a representation has some theoretical advantages [27].

Numerous studies have repeatedly explained the influence of local heterogeneity on flow patterns and stressed the changes in macroscopic properties, such as permeability according to microscale details [28,29,30,31,32,33]. The literature on the quantitative evaluation of the phenomena at pore and throat levels is almost rare.

This study is the continuation of previous research. One of the authors stated the limitations of pore network simulation to predict the permeability of porous media [34]. It was shown that the most significant restriction of current pore network simulation is the lack of a rigorous model for transport phenomena in pore bodies. The pore body’s contribution to total permeation was studied by changing the size of pores relative to throats in model porous media. It was concluded that pore network simulation is reliable as long as this ratio is large. However, suppose the size of pores and throats are at the same order of magnitude. In that case, pore network simulation overestimates permeability since the correlations used for pore bodies cannot model the actual flow patterns in pores. Veyskarami et al. (2016) utilized the pore network simulation to analyze the non-Darcian flow in porous media [35]. Baychev et al. (2019) stated pore network simulation’s reliability mainly depends on the accuracy of the techniques constructed the pore network [36].

Image-based direct numerical simulation methods show that streamlines are highly correlated at high velocities due to inertial core flow and eddies [19,28,34,37]. The jet flow leaving a throat changes the flow patterns dictated by the pore walls. In a purely viscous flow, the fluid follows the walls of the pores, and one expects the maximum shear stresses and energy dissipation near the walls. The creation of jet flow induces the zones of high dissipated energy alongside the jet flow. The incident of such a jet with a wall would also create a high-pressure spot. As a result, the flow entering the downstream throats is influenced by the upstream pore flow.

The correlations describing the flow in each element are limited to the boundary conditions of that element. While mass conservation is guaranteed in pore network simulation using Kirchhoff’s rule, the flow patterns of each flow element are assumed uncorrelated from other components, and the short-range correlation is ignored.

This research’s motivation is to examine pores and throats’ importance numerically. The purpose of this study is to determine which element (pore or throat) is more prone to initiate non-Darcian flow. Forchheimer equation is used as a reference model to select the criteria for termination of Darcy’s law. Each element’s role in the onset of the Forchheimer regime in some model porous media has been investigated. Although the model porous media investigated here is simplistic, the methodology introduced is remarkable. This study highlighted the significance of fluid flow in the pore body in two-dimensional model porous media. The dissipated energy in the pores and throats of the networks was analyzed. Using a vortex identification method, the size of eddies and the related energy dissipation were also studied. The main contribution of the current study is determining the exponent for the pore body when the viscous forces are dominant and its variation by changing the inertia force.

2. Macroscale Analysis

In this section, some macro-scale fluid flow parameters in the networks are investigated. Traditionally, fluid flow in porous media is simulated using the equation of motion, i.e., Darcy law [38]. It assumes a linear dependency of the superficial velocity (U) to the pressure gradient $(\mathrm{U}=\mathrm{m}⁰/(\mathsf{\rho }{\mathrm{A}}_{\mathrm{Network}}$)).

$$-\frac{\Delta \mathrm{p}}{\mathrm{L}}=\frac{\mathsf{\mu }\mathrm{U}}{{\mathrm{K}}_{\mathrm{D}}}$$

(1)

Here, $\mathrm{m}⁰$ is the mass flow rate (kg/s), $\mathsf{\rho }$ is the density (kg/m³), ${\mathrm{A}}_{\mathrm{Network}}$ is the cross section area of the network (m²), $\Delta \mathrm{p}$ is the pressure difference across the network (Pa), L is the length of the network in the main flow direction, $\mathsf{\mu }$ is the fluid’s dynamic viscosity (Pa·s), and ${\mathrm{K}}_{\mathrm{D}}$ is the Darcy permeability (m²). We postulated isotropic porous medium and neglected the gravitational effects.

This linear relationship is not valid at high fluid velocity, and it is necessary to employ the higher-order approximations, such as the Forchheimer equation, i.e., Equation (2) [39].

$$-\frac{\Delta \mathrm{p}}{\mathrm{L}\mathsf{\mu }\mathrm{U}}=\frac{1}{{\mathrm{K}}_{\mathrm{app}}}=\frac{1}{{\mathrm{K}}_{\mathrm{F}}}+\mathsf{\beta }\frac{\mathsf{\rho }\mathrm{U}}{\mathsf{\mu }}$$

(2)

where ${\mathrm{K}}_{\mathrm{app}}$ is the apparent permeability, and ${\mathrm{K}}_{\mathrm{F}}$ and $\mathsf{\beta }$ are parameters of the Forchheimer equation. Knowing $\Delta \mathrm{p}$ and $\mathrm{U}$, all these parameters can be calculated. Considering theoretical deficiencies in deriving Equation (3) [3,24,40,41], we assume Forchheimer equation is authorized. For the sake of generality, the reversibility criteria should be met in isotropic porous media: by reversing the flow direction, the same permeability should be obtained. On this basis, Firdaouss and Guermond showed that permeability must obey odd powers for U [41]. This would be the most significant shortcomings in deriving the Forchheimer equation.

Reynolds ($\mathrm{Re}$) and Forchheimer (Fo) numbers are two parameters frequently used to describe the flow regime in porous media. Different characteristic lengths used to calculate Reynolds number are investigated in [19]. Some authors used the grain size for unconsolidated porous media as the characteristic length (${\mathrm{l}}_{\mathrm{c}}$) to define the Reynolds number [9,19,42]:

$${\mathrm{Re}}_{\mathrm{l}}=\frac{{\mathsf{\rho }\mathrm{Ul}}_{\mathrm{c}}}{\mathsf{\mu }}$$

(3)

For consolidated materials, ${\mathrm{l}}_{\mathrm{c}}$ can be extracted as $\mathsf{\pi }/{\mathrm{S}}_{\mathrm{v}}$, where ${\mathrm{S}}_{\mathrm{v}}$ shows the specific surface of porous media [19].

Some others benefitted the Brinkman screening length ($\sqrt{{\mathrm{K}}_{\mathrm{D}}}$) as the characteristic length [19,43]:

$${\mathrm{Re}}_{\mathrm{d}}=\frac{\mathsf{\rho }\mathrm{U}\sqrt{{\mathrm{K}}_{\mathrm{D}}}}{\mathsf{\mu }}$$

(4)

Some studies have defined the Reynolds number based on the aperture of fracture [5,44].

Subsequently, the Forchheimer number is usually used to investigate the Darcy flow-cessation criterion [23,43,45]:

$${\mathrm{F}}_{\mathrm{O}}=\frac{{\mathsf{\rho }\mathrm{UK}}_{\mathrm{D}}\mathsf{\beta }}{\mathsf{\mu }}$$

(5)

The Forchheimer number shows the ratio of liquid–solid interaction to viscous resistance. Zeng and Grigg proposed a modified version of the Forchheimer number [46]. As Equation (6) shows, the non-Darcy effect parameter (E) has been introduced as the ratio of pressure gradient consumed in overcoming fluid–solid interactions to the overall pressure gradient:

$$\mathrm{E}=\frac{{\mathrm{K}}_{\mathrm{F}}\mathsf{\rho }\mathrm{U}\mathsf{\beta }}{\mathsf{\mu }+{\mathrm{K}}_{\mathrm{F}}\mathsf{\rho }\mathrm{U}\mathsf{\beta }}=\frac{{\mathsf{\rho }\mathrm{U}}^2\mathsf{\beta }}{-\Delta \mathrm{P}/\mathrm{L}}$$

(6)

Then the modified Forchheimer number is:

$${\mathrm{F}}_{\mathrm{O}}^{*}=\frac{\mathrm{E}}{1-\mathrm{E}}$$

(7)

This Forchheimer modified number is more consistent in investigating the non-Darcy onset [46].

By inserting Equations (6) into (7) and dividing ${\mathrm{K}}_{\mathrm{app}}$ by ${\mathrm{K}}_{\mathrm{D}}$, it is easy to prove Equation (8):

$$\frac{{\mathrm{K}}_{\mathrm{app}}}{{\mathrm{K}}_{\mathrm{D}}}={\left(\frac{{\mathsf{\rho }\mathrm{U}\mathsf{\beta }\mathrm{K}}_{\mathrm{D}}}{\mathsf{\mu }}\right)}^{-1}\frac{{\mathrm{F}}_{\mathrm{O}}^{*}}{1+{\mathrm{F}}_{\mathrm{O}}^{*}}={\left({\mathrm{Re}}_{\mathrm{D}}\sqrt{{\mathrm{K}}_{\mathrm{D}}}\mathsf{\beta }\right)}^{-1}\frac{{\mathrm{F}}_{\mathrm{O}}^{*}}{1+{\mathrm{F}}_{\mathrm{O}}^{*}}=\frac{{\mathrm{F}}_{\mathrm{O}}^{*}/{\mathrm{F}}_{\mathrm{O}}}{1+{\mathrm{F}}_{\mathrm{O}}^{*}}=\frac{\mathrm{E}}{{\mathrm{F}}_{\mathrm{O}}}$$

(8)

Equation (7) asserts that plotting $\frac{{\mathrm{K}}_{\mathrm{app}}}{{\mathrm{K}}_{\mathrm{D}}}$ against $\frac{{\mathrm{F}}_{\mathrm{O}}^{*}/{\mathrm{F}}_{\mathrm{O}}}{1+{ \mathrm{F}}_{\mathrm{O}}^{*}}$, should result in a line of slope 1. This equation can serve as a consistency check for data in the Forchheimer regime.

The point where $\frac{{\mathrm{K}}_{\mathrm{app}}}{{\mathrm{K}}_{\mathrm{D}}}=0.99$ is considered as a criterion for the commencement of the non-Darcian flow in this study. The point at which the nonlinear pressure drop contribution is as much as 10% of the total pressure drop is considered by some other researchers as the beginning of the non-Darcian flow (i.e., ${\mathrm{E}}_{\mathrm{cD}}\approx 0.1$) [5,10,46].

3. Modeling

3.1. Geometry of Solid and Void Domains

The geometry used in this article is quite similar to that was considered in our previous research work [34]. A regular two-dimensional pore network with a 30 × 30 size was generated. Consider a periodic array of regular octagons, as shown in Figure 1a. The space between these octagons is the void space for fluid flow. The characteristic dimensions (${\mathrm{d}}_{\mathrm{R}}$, ${\mathrm{l}}_{\mathrm{t}}$ and ${\mathrm{d}}_{\mathrm{t}}$) of pore and throat are shown in Figure 1b. As this figure shows, the expansions and contractions at the throat entry equal half of the pore sizes. A flow path (FP) was formed by a throat and two adjacent pore bodies. The solid phase is shown in white, and the void phase is colored.

The pores’ walls and throats are smooth, and no off-axial links are applicable. In this work, ${\mathrm{l}}_{\mathrm{t}}$ is considered constant (${\mathrm{l}}_{\mathrm{t}}=6 \mathrm{mm}$). Some throats are blocked randomly (27%) to make the network more tortuous. The average coordination number is 2.65.

The dimensional analysis extracts two non-dimensional morphological parameters, i.e., ${\mathrm{L}}_{\mathrm{t}}={\mathrm{l}}_{\mathrm{t}}/{\mathrm{d}}_{\mathrm{t}}$ and ${\mathrm{D}}_{\mathrm{R}}={\mathrm{d}}_{\mathrm{R}}/{\mathrm{d}}_{\mathrm{t}}$.

Table 1. Geometrical parameters of the networks.

Network	${\mathbf{L}}_{\mathbf{t}}{=\mathbf{l}}_{\mathbf{t}}{/\mathbf{d}}_{\mathbf{t}}$	${\mathbf{D}}_{\mathbf{R}}{=\mathbf{d}}_{\mathbf{R}}/{\mathbf{d}}_{\mathbf{t}}$	${\mathbf{D}}_{\mathbf{C}}{=\mathbf{d}}_{\mathbf{C}}/{\mathbf{d}}_{\mathbf{t}}$	Pore Area/Total Void Area	Throat Area/Total Void Area	Pore Area/Throat Area	$$\begin{gathered} {\mathbf{l}}_{\mathbf{C}}/{\mathbf{d}}_{\mathbf{t}} \\ \left(\mathbf{mm}\right) \end{gathered}$$
H1	$6/1$	9.48	7.71	0.856	0.144	5.94	18
H2	$6/2$	5.24	4.58	0.806	0.194	4.14	20
H3	$6/3$	3.82	3.55	0.782	0.218	3.58	23
H4	$6/4$	3.12	3	0.771	0.229	3.36	25

${\mathrm{l}}_{\mathrm{C}}=\pi /{\mathrm{S}}_{\mathrm{v}}$.

summarizes the morphologies used to investigate the effect of different variables. It should be noted that all networks have the same topology. This means that the same throats are blocked in all networks so that the local heterogeneity is the same. Remarkably, ${\mathrm{l}}_{\mathrm{c}}$ are greater than the flow paths.

Each pore is considered as a circular space surrounded by the pore walls. Equation (9) estimates the pore diameter (${\mathrm{d}}_{\mathrm{C}}$):

$${\mathrm{d}}_{\mathrm{C}}=4\frac{\mathrm{A}}{\mathsf{\Pi }}$$

(9)

where A and Π show the pore area and its perimeter, respectively.

3.2. Governing Equations

The single-phase flow of an incompressible fluid inside the generated geometry was modeled through solving mass (Equation (10)) and momentum (Equation (11)) conservation equations in the steady-state situation:

Mass balance equation:

$$\nabla .\mathrm{V}=0$$

(10)

$$V.\nabla V=-\frac{1}{\mathsf{\rho }}\nabla \mathrm{p}+\frac{\mathsf{\mu }}{\mathsf{\rho }}{\nabla }^{2}V$$

(11)

Here, p is the local static pressure, and $V=\mathrm{u}i+\mathrm{w}j$ is a vector that shows the fluid velocity. The viscosity ($\mathsf{\mu }$) and density ($\mathsf{\rho }$) of water as the working fluid are assumed $0.001 \mathrm{Pa}·\mathrm{s}$ and ${998 \mathrm{kg}/\mathrm{m}}^3$, respectively. No-slip boundary condition was applied on the fluid–solid interfaces, i.e., $V=0$. Both upper and lower faces of the generated networks shown in Figure 1a,b are exposed to symmetry boundary conditions. The mass flow rate is assumed constant; a uniform profile is applied for velocity upstream the contraction half pore to model fluid flow in an isolated FP.

The semi-implicit method for pressure-linked equations was employed to solve the pressure and velocity inside the pore space. Steady-state convergence is achieved when the residual for each velocity component becomes less than 1 × 10⁻⁶. The current study uses the quadratic mesh to divide the domain into computing elements.

3.3. Dissipated Energy

Let us consider a control volume in the fluid phase. The dissipation rate in the fluid due to internal viscous forces leads to a pressure drop in the system along the main flow direction [24,47,48].

$$\Delta \mathrm{p}=\frac{\mathsf{\rho }{{\displaystyle \int }}^{ }\mathsf{\phi }\mathrm{dv}}{\mathrm{Q}}$$

(12)

where viscous dissipation is calculated by:

$$\mathsf{\phi }=\mathsf{\mu }\left[{\left({\mathsf{\partial }}_{\mathrm{x}}\mathrm{u}\right)}^2+{\left({\mathsf{\partial }}_{\mathrm{y}}\mathrm{w}\right)}^2+{\left({\mathsf{\partial }}_{\mathrm{x}}\mathrm{w}+{\mathsf{\partial }}_{\mathrm{y}}\mathrm{u}\right)}^2\right]$$

(13)

Q is the mass flow rate throughout the network. Some other methods for averaging are introduced in the literature [20,48]. With enough mesh numbers, $\Delta \mathrm{p}$ should be independent of method of averaging [34]. This study uses $\mathsf{\Phi }=\mathsf{\rho }{{\displaystyle \int }}^{ }\mathsf{\phi }\mathrm{dv}$ as the total dissipation.

3.4. Eddy Identification

It was described in previous sections that the induced energy dissipation by eddies created in the pore volume is one of the phenomena responsible for the non-Darcian flow. To better understand the numbers of the eddy and how they influence the micro and macro flow, it is necessary to compare the core and boundary vortices directly. Indeed, the vortex identification method has been utilized in this regard [49,50,51]. Equations (14) and (15) introduce two scalars ${\mathsf{\Gamma }}_1$ (core vortex) and ${\mathsf{\Gamma }}_2$ (boundary vortex), respectively [49,50,51].

$${\mathsf{\Gamma }}_1\left(x_{\mathrm{P}}\right)=\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{s}}\frac{\left[r_{\mathrm{PM}}\times v_{\mathrm{M}}\right].n}{\left|r_{\mathrm{PM}}\right|.\left|v_{\mathrm{M}}\right|}$$

(14)

$${\mathsf{\Gamma }}_2\left(x_{\mathrm{P}}\right)=\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{s}}\frac{\left[r_{\mathrm{PM}}\times \left(v_{\mathrm{M}}-v_{\mathrm{P}}\right)\right].n}{\left|r_{\mathrm{PM}}\right|.\left|v_{\mathrm{M}}-v_{\mathrm{P}}\right|}$$

(15)

An area (s) with 10 × 10 grid sizes is surrounded by any desired point, such as $x_{\mathrm{P}}$. Moreover, a vector n is normal to the described area. $r_{\mathrm{PM}}$ shows the vector connecting $x_{\mathrm{P}}$ and any other point $x_{\mathrm{M}}$ in an area of s.$v_{\mathrm{P}}$ and $v_{\mathrm{M}}$ are associated with the fluid velocities in the related cells. N is number of cells in domain s, Γ₁ and Γ₂ are parameters defining core and borders of eddies. The iso-counters $\left|{\mathsf{\Gamma }}_2\right|=0.5$ are utilized to detect the vortexes’ boundary. Vortices’ core is identified as the cells of $\left|{\mathsf{\Gamma }}_1\right|=0.8$.

3.5. Dissipation in Throats

Total energy dissipation has two component: pores and throats. In the lights of direct numerical simulation, the magnitude of both parts can be decomposed and analyzed separately. Such a study can be done for more complicated geometries provided that a firm definition for boundaries of pores and throats is used to distinguish the boundaries of pores and throats [36]. At this step, the focus is devoted on determining the relative importance of pores and throats’ resistance against permeation.

Let us assume the Hagen–Poiseuille equation explains fluid flow inside the idealized throats so that [52]

$$\mathsf{\delta }\mathrm{p}=\frac{3\mathsf{\mu }{\mathrm{l}}_{\mathrm{t}}}{{2\mathsf{\rho }\mathrm{b}\mathrm{d}}_{\mathrm{t}}^3}\mathsf{\delta }\mathrm{Q}$$

(16)

where $\mathsf{\delta }\mathrm{Q}$ is the mass flow rate of the fluid that passes through a single throat of length ${\mathrm{l}}_{\mathrm{t}}$, width ${\mathrm{d}}_{\mathrm{t}}$, and depth b. $\mathsf{\delta }\mathrm{p}$ is the corresponding pressure drop. $\mathsf{\delta }\mathrm{Q}$ can be scaled as ${\mathsf{\rho }\mathrm{bd}}_{\mathrm{t}}{\mathrm{v}}_{\mathrm{t}},$ where ${\mathrm{v}}_{\mathrm{t}}$ is velocity in the throat. It should be noted that the Hagen–Poiseuille flow is fully developed but the throats are too short for a fully-developed flow to set in. Thus, this is a simplifying assumption that introduces errors to the analysis.

The Hagen–Poiseuille equation is basically defined to predict the pressure drop in channels under fully developed regime condition. Adloo and Abbasi addressed two deficiencies in using this relation in pore network simulation [34]. Firstly, Hagen–Poiseuille is not authorized to model the pressure drop for short throats and high-velocity fluid. Secondly, it is shown that the local heterogeneity in the pores affects the streamlines entering a throat. Therefore, the fluid may not enter the throat parallel to its main axis. Despite the forgoing shortcomings, pore network simulators benefit from this simplified model [53,54,55].

Rearranging Equation (16) yields:

$$\mathsf{\delta }\mathrm{p}=\frac{3\mathsf{\mu }{\mathrm{l}}_{\mathrm{t}}}{{2\mathrm{d}}_{\mathrm{t}}^2}{\mathrm{v}}_{\mathrm{t}}$$

(17)

which can be written in dimensionless form as:

$$\left(\frac{\mathsf{\delta }\mathrm{p}}{{\mathsf{\rho }\mathrm{v}}_{\mathrm{t}}^2}\right)=12{\mathrm{L}}_{\mathrm{t}}{\mathrm{Re}}^{-1}$$

(18)

Multiplying both sides of Equation (17) by $\mathsf{\delta }\mathrm{Q}$ results in Equation (19).

$$\mathsf{\delta }\mathrm{Q}\mathsf{\delta }\mathrm{p}=\frac{3\mathrm{b}\mathsf{\rho }\mathsf{\mu }{\mathrm{l}}_{\mathrm{t}}}{{2\mathrm{d}}_{\mathrm{t}}}{\mathrm{v}}_{\mathrm{t}}^2$$

(19)

which is the energy dissipation in each throat.

In the same way as Equation (19), since the entire flow passes through the throats, it is plausible to say:

$${\mathsf{\Phi }}_{\mathrm{thrt}}={\left({{\displaystyle \int }}^{ }\mathsf{\phi }\mathrm{dv}\right)}_{\mathrm{thrt}}={\mathsf{\alpha }\mathrm{U}}^{{\mathrm{n}}_{\mathrm{t}}}$$

(20)

where ${\mathsf{\Phi }}_{\mathrm{thrt}}$ is the dissipated energy due to passing fluid through all accessible throats. $\mathsf{\alpha }$ is a proportionality constant that depends on the number of accessible throats, fluid properties, and throat geometry. Use of this general equation would also mimic the possible errors due to use of the Hagen–Poiseuille equation.

3.6. Dissipation in Pores

The pore body is a junction of some throats. In viscous flow, the shear stresses control the flow patterns so that the streamlines are parallel to the pore wall. At high velocities, any disturbance in streamlines coming from throats would disturb the flow in the pore. Pore geometry and fluid separation also severely affect the flow.

To the best of our knowledge, there is no closed model that considers the short-range correlations in the pore body. Still, let us suppose Equation (21) pretends the dissipation in the pore body:

$${\mathsf{\Phi }}_{\mathrm{pore}}={\left({{\displaystyle \int }}^{ }\mathsf{\phi }\mathrm{dv}\right)}_{\mathrm{pore}}={\mathsf{\lambda }\mathrm{U}}^{{\mathrm{n}}_{\mathrm{p}}}$$

(21)

where n_p is the exponent, and $\mathsf{\lambda }$ is a parameter that depends on the number of pores, average coordination number, and flow patterns as well as fluid property and pore characteristic length. Unlike throats for which defining an average velocity is plausible, flow in pores is much more complicated, and determining an average velocity would be too simplistic. Therefore, we postulate that superficial velocity U scales pore velocity.

4. Results and Discussion

4.1. Mesh Study

The computation cells was created by dividing each solid wall into several vertices. The size of the faces was then set. By continuously refining the grids, it was found that increasing the number of vertices above 25/wall and reducing the face size below 3 × 10⁻⁴ m had no significant impact on the permeability and the total dissipation for geometry H4.

4.2. Overall Dissipation

Figure 2a is sketched according to Equation (8) and shows that the extracted permeability data reasonably obey the Forchheimer equation and are consistent.

Table 2. Darcy and Forchheimer parameters for networks.

Network	KD × 109 (m2)	KF × 109 (m2)	β × 10−2 (m−1)
H1	5.4	5.3	167
H2	38.7	38.1	35.4
H3	118	119	15.1
H4	256	257	8.51

reports the parameters for Darcy and Forchheimer’s equations.

Figure 2b–e show the variation of $\frac{{\mathrm{K}}_{\mathrm{app}}}{{\mathrm{K}}_{\mathrm{D}}}$ against various parameters. First of all, Figure 2b,c plots dimensionless permeability versus ${\mathrm{Re}}_{\mathrm{l}}$ and ${\mathrm{Re}}_{\mathrm{D}}$, respectively.

Figure 2b asserts that the larger the permeability of the network is, the smaller the $\mathrm{R}{\mathrm{e}}_{\mathrm{l}}$ for the transition to non-Darcian flow.

In other words, the transition to non-Darcian flow occurs more easily for larger-permeability networks. This observation contradicts some other research observations [9,19]. In contrast, Figure 2c states that the higher the permeability is, the larger the critical $\mathrm{R}{\mathrm{e}}_{\mathrm{D}}$ at which Darcy’s law ceases.

Regardless of the differences in the definition of Reynolds number, both fail to predict a single value for regime change. The exact value for shifting the non-Darcian flow is reported in

Table 3. The exponents of Equations (20) and (21) and the criteria for the onset of non-Darcy regime.

	Throats			Pores				Total
Network	${\mathrm{n}}_{\mathrm{t},\mathrm{D}}$ 1	${\mathrm{n}}_{\mathrm{t},\mathrm{n}\mathrm{D}}$ 2	${\mathrm{E}}_{\mathrm{c}\mathrm{H}\mathrm{P}}$	${\mathrm{n}}_{\mathrm{p},\mathrm{D}}$	${\mathrm{n}}_{\mathrm{p},\mathrm{n}\mathrm{D}}$	$\mathrm{R}{\mathrm{e}}_{\mathrm{p}}$	${\mathrm{E}}_{\mathrm{c}\mathrm{P}}$	${\mathrm{n}}_{\mathrm{n}\mathrm{D}}$	$\mathrm{R}{\mathrm{e}}_{\mathrm{d}}$	$\mathrm{R}{\mathrm{e}}_{\mathrm{l}}$	${\mathrm{F}}_{\mathrm{O}}^{*}$	${\mathrm{E}}_{\mathsf{Í}}$
H1	2	2.07	0.124	2	2.48	0.98	0.0093	2.18	0.025	6.07	0.030	0.030
H2	2	2.11	0.093	2	2.35	1.09	0.0133	2.17	0.047	4.33	0.031	0.031
H3	2	2.14	0.075	2	2.32	1.2	0.023	2.21	0.071	4.27	0.037	0.036
H4	2	2.21	0.07	2	2.3	1.6	0.023	2.3	0.091	4.19	0.039	0.037
$\mathsf{\sigma }/m$ 3			0.23			0.19	0.35		0.43	0.17	0.093	0.096

¹ Exponent in Darcy regime; ² exponent in non-Darcy regime; ³ standard deviation over mean value. Exponents are extracted with more than 99.9% accuracy.

. The last row in Table 3 shows the normalized standard deviation, which is the standard deviation divided by the average value in each column.

In contrast to $\mathrm{R}{\mathrm{e}}_{\mathrm{L}}$ and $\mathrm{R}{\mathrm{e}}_{\mathrm{D}}$, Figure 2d,e shows that ${\mathrm{F}}_{\mathrm{O}}^{*}$ and E are more reliable in foreseeing Darcy’s cessation in this study. Unlike Figure 2b,c, variations in data are more confined, showing that ${\mathrm{F}}_{\mathrm{O}}^{*}$ and E are less affected by morphological parameters, at least far from the percolation threshold [56]. Comparing normalized standard deviations, it was found that parameter $\mathrm{E}$ is more consistent in predicting both the Darcy-cessation and permeability data. Therefore, hereafter, we used $\mathrm{E}$ in our analysis.

Figure 2b,e show the cessation of Darcy flow ($\frac{{\mathrm{K}}_{\mathrm{app}}}{{\mathrm{K}}_{\mathrm{D}}}=0.99$). Circles in Figure 2a indicate the data series that are compared in the next section.

4.3. Throat Dissipation

As reported in Table 3, ${\mathrm{n}}_{\mathrm{t},\mathrm{D}}=2$ in the Darcy regime following the Hagen–Poiseuille model (Equation (19)). As mentioned before, the Hagen–Poiseuille is valid for fully developed systems, but the throats are too short for a fully-developed flow to establish. However, this simplifying assumption may result in errors in the obtained results. The exponent in Equation (20) increases with decreasing ${\mathrm{L}}_{\mathrm{t}}$ in non-Darcy regime, reminding that the accuracy of Hagen–Poiseuille model decreases with decreasing ${\mathrm{L}}_{\mathrm{t}}$ [57].

Figure 3a exhibits the ratio of dissipation in the throat to the total dissipation (${\mathsf{\Phi }}_{\mathrm{t}}/{\mathsf{\Phi }}_{\mathrm{N}}$) as a function of $\mathrm{E}$. It shows that the importance of throat resistance decreases as throat width increases. On the other hand, the significance of throat resistance declines in the non-Darcy regime as $\mathrm{E}$ increases.

In this study, the number and the relative location of open throats in all geometries are similar. However, the fractions of the fluid that passes through them are not equal according to network resistance and tortuosity differences.

Combining parallel tube theory with the Hagen–Poiseuille model suggests that throat dissipation should obey ${\mathsf{\Omega }}_{\mathrm{t}}{\mathrm{U}}^2$ in Darcy regime. Here, ${\mathsf{\Omega }}_{\mathrm{t}}$ is a fitting parameter that is a function of fluid property and the geometry of the porous media. ${\mathsf{\Omega }}_{\mathrm{t}}$ is independent of U in Darcy regime. If ${\mathsf{\Omega }}_{\mathrm{t}}{}_0$ is the value of ${\mathsf{\Omega }}_{\mathrm{t}}$ at infinitesimal velocity, the ratio ${\mathsf{\Omega }}_{\mathrm{t}}/{\mathsf{\Omega }}_{\mathrm{t}}{}_0$ reflects the throats deviation from Hagen–Poiseuille model.

Figure 3b depicts ${\mathsf{\Omega }}_{\mathrm{t}}/{\mathsf{\Omega }}_{\mathrm{t}}{}_0$ as a function of $\mathrm{E}$. As it is seen, for sufficiently small velocities, $\frac{{\mathsf{\Omega }}_{\mathrm{t}}}{{\mathsf{\Omega }}_{\mathrm{t}}{}_0}=1$. However, it increases as E increases and becomes a function of U (${\mathsf{\Omega }}_{\mathrm{t}}={\mathsf{\alpha }\mathrm{U}}^{{\mathrm{n}}_{\mathrm{t}}-2}$). The cessation of Hagen–Poiseuille approximation (${\mathrm{E}}_{\mathrm{cHP}}$) as the point where ${\mathsf{\Omega }}_{\mathrm{t}}/{\mathsf{\Omega }}_{\mathrm{t}}{}_0=0.99$ is indicated. Notably, ${\mathrm{E}}_{\mathrm{cHP}}$ is considerably greater than ${\mathrm{E}}_{\mathrm{cD}}$, which asserts that the cessation of Hagen–Poiseuille equation in throats does not provoke the non-Darcian flow. Moreover, the larger the ${\mathrm{L}}_{\mathrm{t}}$, the more rapidly the Hagen–Poiseuille equation fails (with respect to $\mathrm{E}$).

4.4. Pore Dissipation

Figure 4a,b displays the relative significance of pore body dissipation to total dissipation (${\mathsf{\Phi }}_{\mathrm{p}}/{\mathsf{\Phi }}_{\mathrm{N}}$) and throat dissipation (${\mathsf{\Phi }}_{\mathrm{p}}/{\mathsf{\Phi }}_{\mathrm{t}}$), respectively.

According to Table 1, ${\mathrm{D}}_{\mathrm{R}}$ decreases as the throats become wider, showing an equal weight of pores and throats. Figure 4a,b support this thought. Furthermore, the fraction of pore dissipation increases in the Forchheimer regime with E, in agreement with Figure 3a.

By analogy to throats, it is assumed that ${\mathsf{\Omega }}_{\mathrm{p}}{\mathrm{U}}^2$ simulates pore dissipation in Darcy regime, where ${\mathsf{\Omega }}_{\mathrm{p}}$ is a function of fluid velocity and pore geometry. Defining ${\mathsf{\Omega }}_{\mathrm{p}}{}_0$ as the value of ${\mathsf{\Omega }}_{\mathrm{p}}$ at infinitesimal velocity, the ${\mathsf{\Omega }}_{\mathrm{p}}/{\mathsf{\Omega }}_{\mathrm{p}0}=0.99$ as the termination of applicability of quadratic equation in pores (${\mathrm{E}}_{\mathrm{c}-\mathrm{pore}}$). Above this criterion, the viscous interactions are controlling. Below that, inertial forces become significant. Figure 4c represents the variation of ${\mathsf{\Omega }}_{\mathrm{p}}/{\mathsf{\Omega }}_{\mathrm{p}0}$ as a function of E. At low E, ${\mathsf{\Omega }}_{\mathrm{p}}/{\mathsf{\Omega }}_{\mathrm{p}0}$ equals 1. However, increasing the velocity increases ${\mathsf{\Omega }}_{\mathrm{p}}/{\mathsf{\Omega }}_{\mathrm{p}0}$ severely. Compared with Figure 3b, departure of ${\mathsf{\Omega }}_{\mathrm{p}}/{\mathsf{\Omega }}_{\mathrm{p}0}$ is more observable at high velocities. Moreover, it is evident that for all networks, ${\mathrm{E}}_{\mathrm{c}-\mathrm{pore}}<{\mathrm{E}}_{\mathrm{cD}}$ in contrast to ${\mathrm{E}}_{\mathrm{cHP}}$.

Figure 4c shows that the larger ${\mathrm{D}}_{\mathrm{C}}$ is, the sooner the quadratic equation in the pores fails, in agreement with Figure 2e. At this criterion, a jet flow would leave a throat without influencing by the pore’s wall. The quadratic equation ${\mathsf{\lambda }\mathrm{U}}^2$ stems from viscous forces, which are influenced by the wall’s friction. Larger pores do not meet these criteria at elevated velocities.

Table 3 lists the values of the exponent ${\mathrm{n}}_{\mathrm{p}}$ for Equation (21). Regarding boundary layer theory without flow separation, ${\mathrm{n}}_{\mathrm{p}}=5/2$ [24]. Furthermore, the dissipation in the diverging channel in the vicinity of infinite Reynolds obeys ${\mathrm{C}}_1{\mathrm{U}}^{7/2}+{\mathrm{C}}_2{\mathrm{U}}^{9/4},$ which is higher than $5/2$ [24]. As shown in Table 3, ${\mathrm{n}}_{\mathrm{p}}$ increases with ${\mathrm{D}}_{\mathrm{R}}$, in agreement with Figure 4a.

In the absence of any explicit definition for the pore velocity, the Reynolds number in the pores is defined as follows:

$${\mathrm{Re}}_{\mathrm{p}}=\frac{{\mathsf{\rho }\mathrm{d}}_{\mathrm{C}}\mathrm{U}}{\mathsf{\mu }}$$

(22)

Table 3 shows that in all networks, the departure from the viscous flow starts at $\mathrm{R}{\mathrm{e}}_{\mathrm{p}}\approx 1$. This would support the idea that inertia forces and viscous forces are about the same order of magnitude at this point.

To the best of our knowledge, no study has yet examined the magnitude of energy loss at sub-pore segments (pore and throat separately). However, numerical simulations in somewhat similar geometries have revealed that the inertial effects concentrate in the vicinity of expansion/contraction sections [27].

4.5. Effect of Eddies

Figure 5a shows the ratio of eddy’s cells to the total computation cells (${\mathrm{Area}}_{\mathrm{e}}/{\mathrm{Area}}_{\mathrm{v}}$), and Figure 5b indicates the proportion of eddy’s dissipation to the total energy loss (${\mathsf{\Phi }}_{\mathrm{e}}/{\mathsf{\Phi }}_{\mathrm{N}}$). It can be perceived that eddies occupy less than 2% of the total computation area, and the dissipation due to eddies is negligible to total dissipation, following previous theoretical predictions [24]. Interestingly, vortex dissipation rises drastically even before Darcy law ceases. Besides, the larger the ${\mathrm{D}}_{\mathrm{C} }$, the larger are the eddies and their related dissipation. Therefore, eddies’ effects should be searched by increasing the tortuosity and redirecting the streamlines.

Comparing Figure 5b and Figure 4c, regime changes in pores (increasing the exponent ${\mathrm{n}}_{\mathrm{p}}$) coincide with sudden growth of eddies.

Figure 5b shows that increasing ${\mathrm{D}}_{\mathrm{C}}$ increases ${\mathsf{\Phi }}_{\mathrm{e}}/{\mathsf{\Phi }}_{\mathrm{N}}$ in agreement with exponents of pores (Table 3). However, as shown in Table 3, the exponent for total dissipation varies in the opposite direction and increases as ${\mathrm{D}}_{\mathrm{C}}$ decreases. This observation should be explained according to the throat’s dissipation, which is controlled by the shear layers attached to the wall. By increasing ${\mathrm{D}}_{\mathrm{C}}$, the dissipation shifts towards throats.

4.6. Pore-Scale Analysis

Figure 6 compares the streamlines and dissipation in different networks with approximately the same $\frac{\mathrm{E}}{{\mathrm{F}}_{\mathrm{O}}}$ (series I–IV in Figure 2a). Note that the superficial velocity increases with increasing the series number.

At smaller velocities (smaller $\frac{\mathrm{E}}{{\mathrm{F}}_{\mathrm{O}}}$), the viscous forces control the fluid flow, and as a result, the dissipation’s intensity concentrates in the throats and near the throat entrances. At such a condition, streamlines tend to spread over the entire pore space since inertial forces are not that high to regulate them. They are approximately parallel to each other and follow the solid walls. It can be argued that the flow patterns in the flow paths are independent of upstream flow.

With increasing the velocity, inertia forces become more assertive and are dominant. Therefore, the streamlines associate to form a jet-like flow, which tends to touch the walls less. The jet-flow incident with the walls of the pore body (even in such regular networks) is notable, as shown by purple circles. Moreover, the intensity of dissipation increases in the pore bodies, especially near the jet flow borders. For more elaborated details of the physics, one is referred to [23,27], which insists that the distortion in the pressure and velocity fields causes non-Darcian flow.

Eddy forms when a jet finds the least resistance path, bypassing other portions of the void space. Circulation happens because of shear stresses on the boundaries of a dead volume. Eddies play two crucial roles: (i) decreasing the effective void volume for transportation and (ii) changing the flow direction in flow paths with weak inertia effects. The black circles show how eddies appearing at the entrance of the throats influence the flow direction.

Evidently, as inertia forces become more dominant, the flow patterns in the pore body become more asymmetric. Streamlines do not follow the walls but are affected by the streams entering and leaving the pores. The flow patterns of different networks are not similar to each other when compared at similar $\frac{\mathrm{E}}{{\mathrm{F}}_{\mathrm{O}}}$.

4.7. Discussion

The significance of upstream flow at a throat entrance on the deviation from the Hagen–Poiseuille model has been shown in another study [34]. The departure is meaningful at high velocities when inertial core flow controls the streamlines in the pore bodies. As a result, a large angle between the throat axis and the upstream streamlines would contradict the Hagen–Poiseuille postulation. At lower velocities, streamlines spread over the throat entrance more equally. The information in Table 3 supports the idea that local heterogeneity affects the throat’s conductance at high velocities. Strictly speaking, wider throats are more prone to disobey the Hagen–Poiseuille model.

Because of complicated geometry of real porous media, it would be not easy to put a definite boundary between pores and throats and separate their contribution to fluid flow, subsequently. Nevertheless, some universally accepted schemes, such as medial axis, maximal ball, and watershed, are widely used in the literature [26,58,59,60,61,62].

The literature in pore network simulation includes the correlations to calculate the pressure drop in the half-pore body. The half-pore body means that the diverging/converging segments connected with a throat. Although different geometries are yet considered, and some sophisticated correlations are suggested [26,34,35,58,63,64,65,66], they are all backed up with a critical assumption that the half-pore entry is subjected to a uniform pressure or velocity. Supposing such an assumption simplifies resolving the flow appreciably. Consequently, the local solution and the resulting streamlines are notably independent of other transporting elements in a short distance.

As is proved numerically [20,27], pressure drag has the most crucial role in the pressure gradient. Flow separation downstream from a diverging section creates low pressure and a low-velocity layer attached to the surface [24]. The jet-like flow in the upper layer bends when impinging on an obstacle, which causes a high-pressure, low-velocity zone on the collision area. These complications are ignored in the isolated correlations mentioned above.

5. Conclusions

This study highlighted the significance of fluid flow in the pore body in two-dimensional model porous media. The dissipated energy in the pores and throats of the networks was analyzed. Using a vortex identification method, the size of eddies and the related energy dissipation were also studied.

The Forchheimer equation was utilized to find the onset of non-Darcian flow quantitatively.

1.

The parameter E has less uncertainty in predicting the cessation of Darcy flow. Generally, in the networks studied here, the larger ${\mathrm{D}}_{\mathrm{C}}$ is, the sooner Darcy-flow cessation occurs.

2.

Investigations show that the exponent of the power-law equation for pore body, ${\mathrm{n}}_{\mathrm{p}}$, equals 2 when viscous forces are dominant. The exponent increases as the inertia forces increase. The geometrical parameter ${\mathrm{D}}_{\mathrm{C}}$ has an inverse effect on ${\mathrm{n}}_{\mathrm{p}}$. Moreover, the onset of the increase in the exponent happens earlier than that of the throats. Compared with the beginning of non-Darcian flow in the entire network, it is found that pore flow induces the termination of Darcy flow.

3.

The dissipation due to pore bodies is more apparent when the size of pore and throats are of the same order of magnitude. The relative losses of pore body increase as the velocity increases in contrast to throats.

4.

The area and dissipation due to eddies are almost negligible compared to the total area and dissipation.

5.

The sudden increase in the dissipation due to eddies coincides with the regime change in pores.

6.

As a jet flow appears in a pore, the magnitude of dissipated energy increases in the jet’s borders. The dissipation could be as high as the dissipation in a throat.

7.

To investigate the fluid flow in porous media, the pores and throats should be discretely analyzed using the following formulation:

$${\left({{\displaystyle \int }}^{ }\mathsf{\Phi }\mathrm{dv}\right)}_{\mathrm{total}}={\left({{\displaystyle \int }}^{ }\mathsf{\Phi }\mathrm{dv}\right)}_{\mathrm{pore}}+{\left({{\displaystyle \int }}^{ }\mathsf{\Phi }\mathrm{dv}\right)}_{\mathrm{throat}}={\mathsf{\alpha }\mathrm{U}}^{{\mathrm{n}}_{\mathrm{t}}}+{\mathsf{\lambda }\mathrm{U}}^{{\mathrm{n}}_{\mathrm{p}}}$$

(23)

where ${\mathrm{n}}_{\mathrm{t}}=\mathrm{f}\left({\mathrm{L}}_{\mathrm{t}}\right)\ge 2,$ and ${\mathrm{n}}_{\mathrm{p}}=\mathrm{f}\left({\mathrm{D}}_{\mathrm{C}}\right)\ge 2$. $\mathsf{\alpha }$ and $\mathsf{\lambda }$ are functions of fluid properties and geometries of throats and pores, respectively.

It was acknowledged that utilizing the Hagen–Poiseuille model adds some errors to the achieved results. The future numerical studies could involve Navier–Stokes and continuity equations on either simplified network geometries to reduce the level of observed error.