Origin of the inertial deviation from Darcy's law: An investigation from a microscopic flow analysis on two-dimensional model structures

Inertial flow in porous media occurs in many situations of practical relevance among which one can cite flows in column reactors, in filters, in aquifers, or near wells for hydrocarbon recovery. It is characterized by a deviation from Darcy's law that leads to a nonlinear relationship between the pressure drop and the filtration velocity. In this work, this deviation, also known as the nonlinear, inertial, correction to Darcy's law, which is subject to controversy upon its origin and dependence on the filtration velocity, is studied through numerical simulations. First, the microscopic flow problem was solved computationally for a wide range of Reynolds numbers up to the limit of steady flow within ordered and disordered porous structures. In a second step, the macroscopic characteristics of the porous medium and flow (permeability and inertial correction tensors) that appear in the macroscale model were computed. From these results, different flow regimes were identified: (1) the weak inertia regime where the inertial correction has a cubic dependence on the filtration velocity and (2) the strong inertia (Forchheimer) regime where the inertial correction depends on the square of the filtration velocity. However, the existence and origin of those regimes, which depend also on the microstructure and flow orientation, are still not well understood in terms of their physical interpretations, as many causes have been conjectured in the literature. In the present study, we provide an in-depth analysis of the flow structure to identify the origin of the deviation from Darcy's law. For accuracy and clarity purposes, this is carried out on two-dimensional structures. Unlike the previous studies reported in the literature, where the origin of inertial effects is often identified on a heuristic basis, a theoretical justification is presented in this work. Indeed, a decomposition of the convective inertial term into two components is carried out formally allowing the identification of a correlation between the flow structure and the different inertial regimes. These components correspond to the curvature of the flow streamlines weighted by the local fluid kinetic energy on the one hand and the distribution of the kinetic energy along these lines on the other hand. In addition, the role of the recirculation zones in the occurrence and in the form of the deviation from Darcy's law was thoroughly analyzed. For the porous structures under consideration, it is shown that (1) the kinetic energy lost in the vortices is insignificant even at high filtration velocities and (2) the shape of the flow streamlines induced by the recirculation zones plays an important role in the variation of the flow structure, which is correlated itself to the different flow regimes.

Figure 1

2D model structures of porous media corresponding to ordered (a, b) and disordered (c, d) arrays of parallel cylinders of square (a, c) and circular (b, d) cross sections. (a) $\mathrm{OS}$, (b) $\mathrm{OC}$, (c) $\mathrm{DS},$ and (d) $\mathrm{DC}$. Representative elementary volume (REV) made of $1\times 1$ (a, b) and $30\times 30$ (c, d) geometrical unit cells. Porosity $\epsilon =75\%$.

Figure 2

Large (macroscopic) scale structure ($\mathrm{OS}$) and the periodic unit cell that corresponds to the REV for steady flow.

Figure 3

Normalized intrinsic average of the $x$ component of the velocity, ${\left({〈v_x^{*}〉}^{\beta }\right)}_n$, versus the number of grid blocks. Model porous structures of Fig. 1. Normalization with respect to ${〈v_x^{*}〉}^{\beta }$ obtained with the finest mesh. $\mathrm{OS}, {\mathrm{Re}}^{*}=17500$ (${\mathrm{Re}}_{\mathrm{k}}=23.41$); $\mathrm{OC}, {\mathrm{Re}}^{*}=17500$ (${\mathrm{Re}}_{\mathrm{k}}=21.09$); $\mathrm{DS}, {\mathrm{Re}}^{*}=30000$ (${\mathrm{Re}}_{\mathrm{k}}=18.58$); $\mathrm{DC}, {\mathrm{Re}}^{*}=30000$ (${\mathrm{Re}}_{\mathrm{k}}=19.87$); $\epsilon =75\%$; $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x$.

Figure 4

$x$ component of the inertial correction vector, $f_{cx}$, versus the Reynolds number ${\mathrm{Re}}_{\mathrm{k}}$ for the model porous structures of Fig. 1 ($\mathrm{OS}, \mathrm{OC}, \mathrm{DS},$ and $\mathrm{DC}$) and two different orientations of $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }$: $\theta =0^{\circ }$ ($\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x$) and $\theta ={45}^{\circ }$ [$\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }=\left({\widehat{\mathbf{e}}}_x+{\widehat{\mathbf{e}}}_y\right)/\sqrt{2}$]. Porosity $\epsilon =75\%$.

Figure 5

Fluid kinetic energy color map, ${\mathbf{v}}^{*2}$, and steady flow streamlines in the $\mathrm{OS}$ of Fig. 1. The thick red streamlines delimit recirculation zones. The dotted line represents an arbitrary cross section $S\perp$ to $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }$. (a) ${\mathrm{Re}}_{\mathrm{k}}=0$, (b) ${\mathrm{Re}}_{\mathrm{k}}=23.41$, (c) ${\mathrm{Re}}_{\mathrm{k}}=0$, (d) ${\mathrm{Re}}_{\mathrm{k}}=6.76$. (a and b) $\theta =0^{\circ }$ ($\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x$). (c and d), $\theta ={45}^{\circ } [\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }=\left({\widehat{\mathbf{e}}}_x+{\widehat{\mathbf{e}}}_y\right)/\sqrt{2}]$. Porosity $\epsilon =75\%$.

Figure 6

First derivative of the $x$ component of the inertial correction normalized by its maximum value ${\left(\partial f_{cx}/\partial {\mathrm{Re}}_{\mathrm{k}}\right)}_n$ versus ${\mathrm{Re}}_{\mathrm{k}}$ for the model porous structures of Fig. 1 ($\mathrm{OS}, \mathrm{OC}, \mathrm{DS},$ and $\mathrm{DC}$) and two different orientations of $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }$: $\theta =0^{\circ } (\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x)$ and $\theta ={45}^{\circ } [\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }=\left({\widehat{\mathbf{e}}}_x+{\widehat{\mathbf{e}}}_y\right)/\sqrt{2}]$. Porosity $\epsilon =75\%$.

Figure 7

Ratio between the fluid kinetic energy in the vortices and the total fluid kinetic energy, $E{c}_{\mathrm{Vort}}^{*}/E{c}^{*}$, versus the Reynolds number ${\mathrm{Re}}_{\mathrm{k}}$. Ordered model porous structures $\mathrm{OS}$ and $\mathrm{OC}$ of Fig. 1. Two different orientations of the macroscopic pressure gradient, $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$. Identification of flow regimes from data in Table 1. Porosity $\epsilon =75\%$.

Figure 8

Ratio between the volume occupied by vortices and the total volume of fluid $\beta , V_{\beta \mathrm{Vort}}/V_{\beta }$, versus the Reynolds number ${\mathrm{Re}}_{\mathrm{k}}$. Ordered model porous structures $\mathrm{OS}$ and $\mathrm{OC}$ of Fig. 1. Two different orientations of the macroscopic pressure gradient, $\theta =0^{\circ }$ and $\theta ={45}^{\circ }$. Identification of flow regimes from data in Table 1. Porosity $\epsilon =75\%$.

Figure 9

Streamlines tortuosity, $T$, versus the Reynolds number ${\mathrm{Re}}_{\mathrm{k}}$ for the model porous structures of Fig. 1 ($\mathrm{OS}, \mathrm{OC}, \mathrm{DS},$ and $\mathrm{DC}$) and two different orientations of $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }$: $\theta =0^{\circ }$ ($\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x$), and $\theta ={45}^{\circ } [\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }=\left({\widehat{\mathbf{e}}}_x+{\widehat{\mathbf{e}}}_y\right)/\sqrt{2}]$. Identification of flow regimes from data in Table 1. Porosity $\epsilon =75\%$.

Figure 10

Variance of the fluid kinetic energy, $\mathtt{Var}\left({\mathbf{v}}^{*2}\right)$, versus the Reynolds number ${\mathrm{Re}}_{\mathrm{k}}$ for the model porous structures of Fig. 1 ($\mathrm{OS}, \mathrm{OC}, \mathrm{DS},$ and $\mathrm{DC}$) and two different orientations of $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }$: $\theta =0^{\circ }$ ($\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x$) and $\theta ={45}^{\circ } [\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }=\left({\widehat{\mathbf{e}}}_x+{\widehat{\mathbf{e}}}_y\right)/\sqrt{2}]$. Identification of flow regimes from data in Table 1. Porosity $\epsilon =75\%$.

Figure 11

Intrinsic average of the absolute value of the streamlines curvature weighted by the local kinetic energy, ${〈\left|{\mathbf{v}}^{*2}\kappa \widehat{\mathbf{n}}\right|〉}^{\beta }$, versus the Reynolds number ${\mathrm{Re}}_{\mathrm{k}}$ for the model porous structures of Fig. 1 ($\mathrm{OS}, \mathrm{OC}, \mathrm{DS},$ and $\mathrm{DC}$) and two different orientations of $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }$: $\theta =0^{\circ }$ ($\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x$) and $\theta ={45}^{\circ } [\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }=\left({\widehat{\mathbf{e}}}_x+{\widehat{\mathbf{e}}}_y\right)/\sqrt{2}]$. Identification of flow regimes from data in Table 1. Porosity $\epsilon =75\%$.

Figure 12

Variation of the kinetic energy along the streamlines in the direction of flow, $\frac{1}{2}{〈\left|\frac{d{\mathbf{v}}^{*2}}{ds}\widehat{\mathbf{t}}\right|〉}^{\beta }$, versus the Reynolds number ${\mathrm{Re}}_{\mathrm{k}}$ for the model porous structures of Fig. 1 ($\mathrm{OS}, \mathrm{OC}, \mathrm{DS},$ and $\mathrm{DC}$) and two different orientations of $\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }$: $\theta =0^{\circ }$ ($\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }={\widehat{\mathbf{e}}}_x$) and $\theta ={45}^{\circ } [\mathbf{\nabla }{\langle p^{*}\rangle }^{\beta }=\left({\widehat{\mathbf{e}}}_x+{\widehat{\mathbf{e}}}_y\right)/\sqrt{2}]$. Identification of flow regimes from data in Table 1.