2 Asymptotic expansion of the governing equations for low-Mach-number flow

Pore-scale flow of a compressible fluid is governed by the continuity and compressible Navier–Stokes equations (Anderson Reference Anderson1995)

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{v})=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D70C}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}=-\unicode[STIX]{x1D735}p+\left(\unicode[STIX]{x1D707}_{b}+\frac{1}{3}\unicode[STIX]{x1D707}\right)\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v})+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{v},~p,~\unicode[STIX]{x1D70C}$ are the velocity, pressure and density, and $\unicode[STIX]{x1D707},~\unicode[STIX]{x1D707}_{b}$ are the shear and bulk viscosity of the fluid, respectively. Temperature change is assumed to be negligible. We follow Klainerman & Majda (Reference Klainerman and Majda1982) in making the governing equations dimensionless. Let $v_{ref}$ be the magnitude of a reference fluid particle velocity, $L_{ref}$ a reference length, and a reference time is defined as $t_{ref}=L_{ref}/v_{ref}$ . The reference pressure and density are $p_{ref}$ and $\unicode[STIX]{x1D70C}_{ref}$ , respectively. A reference speed of sound can be defined as $c_{ref}=\sqrt{p_{ref}/\unicode[STIX]{x1D70C}_{ref}}$ (a constant factor is dropped for convenience; see Klainerman & Majda (Reference Klainerman and Majda1982)). We define dimensionless velocity, pressure, density, coordinates and time as $\bar{\boldsymbol{v}}=\boldsymbol{v}/v_{ref},~\bar{p}=p/p_{ref},~\bar{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{ref},~\bar{\boldsymbol{x}}=\boldsymbol{x}/L_{ref}$ and $\bar{t}=t/t_{ref}$ . A global Mach number is defined as $M=v_{ref}/c_{ref}$ , and it will be assumed small. A Reynolds number is defined as $Re=\unicode[STIX]{x1D70C}_{ref}v_{ref}L_{ref}/\unicode[STIX]{x1D707}$ . The dimensionless equations are

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}\bar{t}}+\bar{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }(\bar{\unicode[STIX]{x1D70C}}\bar{\boldsymbol{v}})=0, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\bar{\boldsymbol{v}}}{\unicode[STIX]{x2202}\bar{t}}+\bar{\unicode[STIX]{x1D70C}}\bar{\boldsymbol{v}}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D735}}\bar{\boldsymbol{v}}=-\frac{1}{M^{2}}\bar{\unicode[STIX]{x1D735}}\bar{p}+\frac{\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+1/3}{Re}\bar{\unicode[STIX]{x1D735}}(\bar{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\bar{\boldsymbol{v}})+\frac{1}{Re}\bar{\unicode[STIX]{x1D6FB}}^{2}\bar{\boldsymbol{v}}. & \displaystyle\end{eqnarray}$$

Using these equations and for inviscid isentropic flow, Klainerman & Majda (Reference Klainerman and Majda1982) have rigorously proved that, for given initial data, the linearized acoustics is a uniformly valid principal correction in the deviation of the compressible flow solution from the incompressible solution as the global Mach number $M\rightarrow 0$ , i.e.

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle |\bar{\boldsymbol{v}}-(\bar{\boldsymbol{v}}^{I}+M\,\tilde{\boldsymbol{v}}^{A})|\leqslant \unicode[STIX]{x1D6E5}_{5}M^{2}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle |\bar{p}-(\bar{p}_{0}+M^{2}(\tilde{p}^{I}+\tilde{p}^{A}))|\leqslant \unicode[STIX]{x1D6E5}_{6}M^{3}, & \displaystyle\end{eqnarray}$$

where $|~|$  stands for the magnitude, and superscripts ‘ $I$ ’ and ‘ $A$ ’ stand for incompressible and acoustic quantities, respectively; $\bar{\boldsymbol{v}}^{I}=\lim _{M\rightarrow 0}\bar{\boldsymbol{v}}$ with $\bar{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\bar{\boldsymbol{v}}^{I}=0$ is the incompressible solution at vanishing Mach number; $\unicode[STIX]{x1D6E5}_{5}>0,~\unicode[STIX]{x1D6E5}_{6}>0$ are fixed constants, and $\bar{p}_{0}$ is a background thermodynamic pressure. This asymptotic result has been extended by Klein (Reference Klein1995) to non-isentropic flows and by Munz et al. (Reference Munz, Roller, Klein and Geratz2003) and Munz, Dumbser & Rolle (Reference Munz, Dumbser and Rolle2007) to the general heat-conducting viscous compressible Navier–Stokes equations. These asymptotic results have been used to construct efficient numerical schemes such as the multiple-pressure-variable method to overcome slow convergence, even failure, of fully compressible solvers for small-Mach-number flows (Klein Reference Klein1995; Klein & Munz Reference Klein and Munz1995; Bijl & Wesseling Reference Bijl and Wesseling1998; Bailly, Bogey & Juve Reference Bailly, Bogey and Juve1999; Wesseling Reference Wesseling2001; Munz et al. Reference Munz, Roller, Klein and Geratz2003, Reference Munz, Dumbser and Rolle2007).

The above asymptotic theory allows us to decompose a low-Mach-number flow into the sum of two fields, the limiting incompressible field $\bar{\boldsymbol{v}}^{I}$ and a linearized acoustic field $\bar{\boldsymbol{v}}^{A}$ ,

(2.7) $$\begin{eqnarray}\bar{\boldsymbol{v}}=\bar{\boldsymbol{v}}^{I}+\bar{\boldsymbol{v}}^{A}.\end{eqnarray}$$

For primary recovery, the outer boundary of the reservoir is non-penetrable. Thus, except for isolated dead pores, all pores in the reservoir are half-sealed, with the sealed ends located on the outer boundary or in the interior of the reservoir, and the open ends located at the wellbore where the fluid is being produced. Therefore, an incompressible velocity field does not contribute to the mass flow rate at the wellbore, since integration of the incompressibility condition over the connected pore space gives the corresponding volumetric flow rate at the wellbore as zero:

(2.8) $$\begin{eqnarray}Q\;^{I}=\int _{wellbore}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{v}^{I}\,\text{d}a=-\int _{outer\;boundary}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{v}^{I}\,\text{d}a=0,\end{eqnarray}$$

where $\boldsymbol{n}$ is the unit outward normal vector on the corresponding boundaries. Since the flow in half-sealed pores relies entirely on the volumetric expansion or compressibility of the fluid, the limiting incompressible velocity $\bar{\boldsymbol{v}}^{I}=\lim _{M\rightarrow 0}\bar{\boldsymbol{v}}$ must be identically zero in primary recovery, $\bar{\boldsymbol{v}}^{I}=0$ . For $\bar{\boldsymbol{v}}^{I}$ to be non-zero, the pores must have two open ends. Physically, an incompressible fluid cannot expand in volume; and in half-sealed pores, there is no fresh fluid to take the place of the volume in the pore space that would have been vacated by the produced fluid. In other words, to leading order, flow in primary recovery is completely characterized by the acoustic field $\bar{\boldsymbol{v}}^{A}$ (since $\bar{\boldsymbol{v}}^{I}=0$ ), which is coupled to the density change of the fluid. Primary recovery is determined by how a fluid expands in volume in semi-sealed pores, not by how it is pushed out of the pores such as in a tube with open ends.

The acoustic velocity field $\bar{\boldsymbol{v}}^{A}$ can be further decomposed into the sum of an irrotational part $\bar{\boldsymbol{v}}_{IR}$ and a rotational part $\bar{\boldsymbol{v}}_{RT}$ using the Helmholtz decomposition theorem (Aris Reference Aris1989; Leal Reference Leal2010; Panton Reference Panton2013),

(2.9) $$\begin{eqnarray}\bar{\boldsymbol{v}}^{A}=\bar{\boldsymbol{v}}_{IR}+\bar{\boldsymbol{v}}_{RT},\end{eqnarray}$$

where the irrotational part is a potential flow and the rotational part is solenoidal (divergence-free):

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\times \bar{\boldsymbol{v}}_{IR}=0,\quad \bar{\boldsymbol{v}}_{IR}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}, & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\boldsymbol{v}}_{RT}=0. & \displaystyle\end{eqnarray}$$

In the above, $\unicode[STIX]{x1D6F7}$ is the scalar velocity potential for the irrotational part of the velocity. Thus, for primary recovery, $\bar{\boldsymbol{v}}^{I}=0$ , and the asymptotic results (2.5) and (2.6) suggest the following small-Mach-number expansions for the velocity and pressure:

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\boldsymbol{v}}=\bar{\boldsymbol{v}}^{A}+O(M^{2})=M\,\tilde{\boldsymbol{v}}_{RT}+M\,\tilde{\boldsymbol{v}}_{IR}+O(M^{2}), & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{p}-\bar{p}_{0}=M^{2}\tilde{p}_{1}+M^{2}\tilde{p}^{A}+O(M^{3}), & \displaystyle\end{eqnarray}$$

where $\bar{\boldsymbol{v}}_{RT}=M\,\tilde{\boldsymbol{v}}_{RT}$ and $\bar{\boldsymbol{v}}_{IR}=M\,\tilde{\boldsymbol{v}}_{IR}$ . The pressure $\bar{p}_{1}=M^{2}\,\tilde{p}_{1}$ is the Lagrangian multiplier associated with the solenoidal field $\bar{\boldsymbol{v}}_{RT}$ to ensure that the incompressibility condition for $\bar{\boldsymbol{v}}_{RT}$ is satisfied. The pressure $\bar{p}_{1}$ is thus a hydrodynamic pressure induced by the solenoidal field $\bar{\boldsymbol{v}}_{RT}$ and it is not described by the equation of state. The necessity to include $\tilde{p}_{1}$ in the pressure expansion has been discussed by Klein (Reference Klein1995) and Munz et al. (Reference Munz, Roller, Klein and Geratz2003, Reference Munz, Dumbser and Rolle2007). The estimates (2.5) and (2.6) guarantee the convergence of the series (2.12) and (2.13) in the limit $M\rightarrow 0$ . Density can be similarly expanded as

(2.14) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}=1+M^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}+O(M^{4}).\end{eqnarray}$$

Density change is only related to the acoustic pressure $\bar{p}^{A}=M^{2}\tilde{p}^{A}$ via the equation of state

(2.15) $$\begin{eqnarray}\tilde{p}^{A}=\tilde{\unicode[STIX]{x1D70C}}_{2}.\end{eqnarray}$$

When the expansions (2.12) and (2.13) are substituted into the governing equations, to leading order we obtain the following equations for the irrotational field:

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70C}}_{2}}{\unicode[STIX]{x2202}\tilde{t}}+\bar{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\tilde{\boldsymbol{v}}_{IR}=0, & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{\boldsymbol{v}}_{IR}}{\unicode[STIX]{x2202}\tilde{t}}=-\bar{\unicode[STIX]{x1D735}}\tilde{p}^{A}+\frac{\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3}{Re_{AC}}\bar{\unicode[STIX]{x1D735}}(\bar{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\tilde{\boldsymbol{v}}_{IR}), & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{p}^{A}=\tilde{\unicode[STIX]{x1D70C}}_{2}, & \displaystyle\end{eqnarray}$$

where the fast-time variable $\tilde{t}=\bar{t}/M=t/t_{ref,AC}$ , with $t_{ref,AC}=L_{ref}/c_{ref}$ being the acoustic time scale; and the acoustic Reynolds number $Re_{AC}=\unicode[STIX]{x1D70C}_{ref}c_{ref}L_{ref}/\unicode[STIX]{x1D707}$ , which is related to the Reynolds number by $Re=MRe_{AC}$ . The leading-order solenoidal field equations are

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\bar{\boldsymbol{v}}_{RT}=0, & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{\boldsymbol{v}}_{RT}}{\unicode[STIX]{x2202}\bar{t}}=-\bar{\unicode[STIX]{x1D735}}\tilde{p}_{1}+\frac{1}{Re}\bar{\unicode[STIX]{x1D6FB}}^{2}\bar{\boldsymbol{v}}_{RT}. & \displaystyle\end{eqnarray}$$

In the acoustic literature, the irrotational field and the solenoidal field are called the longitudinal mode and the transverse mode, respectively (Morse & Ingard Reference Morse and Ingard1968; Pierce Reference Pierce1981).

It is noted that the solenoidal field is absent in the leading-order irrotational velocity equations (2.16)–(2.18). The irrotational velocity, however, should not be expected to satisfy the no-slip condition on a solid boundary. The solenoidal field $\bar{\boldsymbol{v}}_{RT}$ is induced by the irrotational velocity $\bar{\boldsymbol{v}}_{IR}$ at the boundary so that the no-slip condition for the overall velocity $\bar{\boldsymbol{v}}^{A}=\bar{\boldsymbol{v}}_{RT}+\bar{\boldsymbol{v}}_{IR}$ on the pore wall is satisfied:

(2.21) $$\begin{eqnarray}\bar{\boldsymbol{v}}_{RT}=-\bar{\boldsymbol{v}}_{IR}.\end{eqnarray}$$

For primary recovery, the solenoidal field $\bar{\boldsymbol{v}}_{RT}$ does not generate a net flow rate (see (2.8)).

It is noted that mode decomposition of the linear acoustic field has been used for a long time and it was examined theoretically in detail by Lagerstrom, Cole & Trilling (Reference Lagerstrom, Cole and Trilling1949) and Wu (Reference Wu1956). However, there is one major difference between these earlier works and the one presented here: in the earlier studies, the pressure $\tilde{p}_{1}$ in the solenoidal field momentum equation (2.20) has been set to zero identically, $\tilde{p}_{1}=0$ (Lagerstrom et al. Reference Lagerstrom, Cole and Trilling1949; Wu Reference Wu1956; Morse & Ingard Reference Morse and Ingard1968; Pierce Reference Pierce1981; Temkin Reference Temkin1981; Friend & Yeo Reference Friend and Yeo2011). However, without this induced pressure, which serves as a Lagrangian multiplier to enforce the incompressibility condition for the solenoidal field, the solenoidal equations (2.19) and (2.20) are, in general, overdetermined, as there would be four equations for three unknowns (the three velocity components). Thus, in the absence of $\tilde{p}_{1}$ , except for a few special cases, the set of equations (2.19) and (2.20) is likely to give rise to solutions that are internally incompatible due to overconstraint. In fact, a key result of the asymptotic theory of Klainerman & Majda (Reference Klainerman and Majda1982) is the splitting of the pressure into the sum of a thermodynamic pressure, a hydrodynamic pressure and an acoustic pressure for small-Mach-number flows. The distinct multiple roles that pressure plays physically in the small-Mach-number limit of a compressible flow and the necessity to include $\tilde{p}_{1}$ in the pressure expansion have been discussed in detail in Klein (Reference Klein1995), Klein & Munz (Reference Klein and Munz1995) and Munz et al. (Reference Munz, Roller, Klein and Geratz2003, Reference Munz, Dumbser and Rolle2007). These results form the basis for the multiple-pressure-variable numerical method for small-Mach-number flows (Klein Reference Klein1995; Klein & Munz Reference Klein and Munz1995; Bijl & Wesseling Reference Bijl and Wesseling1998; Bailly et al. Reference Bailly, Bogey and Juve1999; Klein et al. Reference Klein, Botta, Schneider, Munz, Roller, Meister, Hoffmann and Sonar2001; Wesseling Reference Wesseling2001; Munz et al. Reference Munz, Roller, Klein and Geratz2003, Reference Munz, Dumbser and Rolle2007).

3 The damped wave equation for density and its asymptotic form at large times

For primary fluid recovery, integration of the continuity equation (2.1) in the entire interconnected pore space gives the mass flow rate at the wellbore as

(3.1) $$\begin{eqnarray}{\dot{m}}=-\int _{pore\;space}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}\,\text{d}V,\end{eqnarray}$$

where the property that the outer boundary is non-penetrable has been used. Thus, density change in the pore space determines the mass production rate at the wellbore. Density change appears only in the equations (2.16)–(2.18) for the irrotational field, which leads to a decoupled damped wave equation for the density:

(3.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}}{\unicode[STIX]{x2202}\tilde{t}^{2}}=\bar{\unicode[STIX]{x1D6FB}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}+\frac{\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3}{Re_{AC}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\tilde{t}}\bar{\unicode[STIX]{x1D6FB}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}.\end{eqnarray}$$

This equation is based on the fast acoustic time scale $\tilde{t}$ . The dimensional form of the damped wave equation is

(3.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t^{2}}=c_{ref}^{2}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70C}+D_{\unicode[STIX]{x1D70C},0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70C},\end{eqnarray}$$

where

(3.4) $$\begin{eqnarray}D_{\unicode[STIX]{x1D70C},0}=\frac{\unicode[STIX]{x1D707}_{b}+4\unicode[STIX]{x1D707}/3}{\unicode[STIX]{x1D70C}_{ref}}=(\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3)\unicode[STIX]{x1D708}_{ref},\end{eqnarray}$$

with $\unicode[STIX]{x1D708}_{ref}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}_{ref}$ being the kinematic viscosity. Here $D_{\unicode[STIX]{x1D70C},0}$ is a damping coefficient, as the last term on the right-hand side of (3.3) represents viscous damping of the acoustic (density) waves; $D_{\unicode[STIX]{x1D70C},0}$ is also the self-diffusion coefficient of the fluid, as the damping stems from the self-diffusion of the fluid mass due to a local density gradient. When the bulk viscosity is zero as in the Stokes hypothesis, $D_{\unicode[STIX]{x1D70C},0}=4\unicode[STIX]{x1D708}_{ref}/3$ . This result is consistent with the kinetic theory of gases, which shows that the self-diffusivity of a gas is proportional to the momentum diffusion coefficient $\unicode[STIX]{x1D708}$ (Hirschfelder, Curtiss & Bird Reference Hirschfelder, Curtiss and Bird1954). Damping presented by the last term on the right-hand side of (3.3) is also called sound absorption, as the acoustic waves are absorbed by the bulk of the fluid. In the absence of a source for continuous excitation, an acoustic wave is expected to be absorbed by the fluid as well as the boundaries and damped out over large times (here ‘large’ is relative to the acoustic time scale). The linear damped wave equation (3.3), called the Stokes equation by Rayleigh (Reference Rayleigh1945), is well known in the acoustic literature and it has appeared in many textbooks (e.g. Morse & Ingard Reference Morse and Ingard1968; Temkin Reference Temkin1981).

When the damped wave equation (3.2) is expressed in terms of the slow time variable $\bar{t}=M\,\tilde{t}$ , which is suitable for large times, we have

(3.5) $$\begin{eqnarray}M^{2}\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}}{\unicode[STIX]{x2202}\bar{t}^{2}}=\bar{\unicode[STIX]{x1D6FB}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}+M^{2}\frac{\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3}{Re}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\bar{t}}\bar{\unicode[STIX]{x1D735}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2},\end{eqnarray}$$

where the relation $Re=MRe_{AC}$ has been used. Thus, when viewed on the large time scale, to leading order in Mach number,

(3.6) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D6FB}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}=0.\end{eqnarray}$$

This is the equation for $\bar{t}\rightarrow \infty$ that corresponds to the final equilibrium state. At the next order, i.e. for finite $\bar{t}$ , we have a diffusion equation,

(3.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}}{\unicode[STIX]{x2202}\bar{t}^{2}}=\frac{\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3}{Re}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\bar{t}}\bar{\unicode[STIX]{x1D735}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}.\end{eqnarray}$$

Consider an initial value problem starting from the state of rest. Then, the initial condition $\tilde{t}=0:~\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70C}}_{2}/\unicode[STIX]{x2202}\tilde{t}=0$ can be approximated as $\bar{t}=0:~\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70C}}_{2}/\unicode[STIX]{x2202}\bar{t}=0$ for the diffusion equation (3.7) if the acoustic waves are damped rapidly when measured on the large time scale $\bar{t}$ . Integration of the diffusion equation (3.7) yields the standard diffusion equation

(3.8) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70C}}_{2}}{\unicode[STIX]{x2202}\bar{t}}=\frac{\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3}{Re}\bar{\unicode[STIX]{x1D6FB}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2}.\end{eqnarray}$$

Thus, at large times, the density wave becomes purely diffusive, controlled by diffusion with a diffusive velocity scale $v_{ref}=D_{\unicode[STIX]{x1D70C},0}/L_{ref}$ . The Reynolds number is then $Re=\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3$ , and the dimensionless diffusion equation (3.8) becomes

(3.9) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70C}}_{2}}{\unicode[STIX]{x2202}\bar{t}}=\bar{\unicode[STIX]{x1D735}}^{2}\tilde{\unicode[STIX]{x1D70C}}_{2},\end{eqnarray}$$

where $\bar{t}=t/t_{ref}=D_{\unicode[STIX]{x1D70C},0}t/L_{ref}^{2}$ . Once the density is solved from the diffusion equation (3.9), the irrotational velocity can be obtained from the continuity equation (2.16) and the solenoidal field from solving (2.19) and (2.20).

It should be noted that, with an error of the order of $O(M^{4})$ , equation (3.9) can be written as

(3.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}\bar{t}}=\bar{\unicode[STIX]{x1D6FB}}^{2}\ln \bar{\unicode[STIX]{x1D70C}}.\end{eqnarray}$$

In dimensional form,

(3.11) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(D_{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}),\end{eqnarray}$$

where

(3.12) $$\begin{eqnarray}D_{\unicode[STIX]{x1D70C}}=(\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3)\unicode[STIX]{x1D708},\end{eqnarray}$$

with $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}$ being the variable kinematic viscosity, is a density-dependent diffusion coefficient. It can be easily shown from (3.11) and the continuity equation that the mass flux of the irrotational field $\dot{\boldsymbol{m}}_{IR}=\unicode[STIX]{x1D70C}\boldsymbol{v}_{IR}$ for primary recovery is given by a Fick’s law type of relation,

(3.13) $$\begin{eqnarray}\dot{\boldsymbol{m}}_{IR}=-D_{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}.\end{eqnarray}$$

Since the solenoidal field does not contribute to the mass flow rate in primary recovery, the irrotational mass flux $\dot{\boldsymbol{m}}_{IR}$ is also the pore-scale mass flux $\dot{\boldsymbol{m}}=\unicode[STIX]{x1D70C}\boldsymbol{v}$ for primary recovery, i.e. $\dot{\boldsymbol{m}}=\dot{\boldsymbol{m}}_{IR}$ .

The above analysis shows that, at large times, density in the pore space diffuses according to the diffusion equation (3.11) with the diffusion coefficient given by (3.12). We emphasize that the diffusive time scale $t_{ref}=L_{ref}^{2}/D_{\unicode[STIX]{x1D70C},0}$ is always much larger than the acoustic time scale $t_{ref,AC}=L_{ref}/c_{ref}$ ; and ‘large’ or ‘long’ time is measured relative to the acoustic time scale. The density diffusion equation (3.11) is derived from an asymptotic analysis of the compressible Navier–Stokes equations at low Mach numbers, and it is the fundamental pore-scale transport equation from which upscaling can be performed.

For a semi-sealed cylindrical pore with a radius $R$ and length $L$ , Chen & Shen (Reference Chen and Shen2018a ,Reference Chen and Shen b ) obtained the complete solution for drainage flow using the linearized damped wave equation (3.3). A prominent feature of the pore-scale flow is that density relaxes quickly in the transverse direction after the start-up and density becomes uniform in the tube cross-section. The velocity field satisfies the no-slip condition and the fluid production rate at large times is given by

(3.14) $$\begin{eqnarray}{\dot{m}}_{e}(t)=2M_{L}\frac{D_{\unicode[STIX]{x1D70C},0}}{L^{2}}\mathop{\sum }_{n=0}^{\infty }\exp \left[-\frac{(2n+1)^{2}\unicode[STIX]{x03C0}^{2}D_{\unicode[STIX]{x1D70C},0}t}{4L^{2}}\right],\end{eqnarray}$$

where $M_{L}=\unicode[STIX]{x03C0}R^{2}L(\unicode[STIX]{x1D70C}_{i}-\unicode[STIX]{x1D70C}_{e})$ is the amount of producible fluid for the given density drop; and $\unicode[STIX]{x1D70C}_{i}$ and $\unicode[STIX]{x1D70C}_{e}$ are the initial density and exit density, respectively. The same result can be obtained by solving the linearized density diffusion equation (3.9). Clearly, the fluid production rate exhibits a slip-like behaviour as it is proportional to $R^{2}$ . The fluid production rate is proportional to the fluid’s kinematic viscosity instead of the reciprocal of the kinematic viscosity as provided by the classical Poiseuille’s law. These pore-scale results are indicative of the invalidity of the Poiseuille–Darcy framework for flow in primary recovery at the pore scale, and consequently at the macroscale.

4 Macroscopic equations for primary recovery

At the pore scale, density change is controlled by self-diffusion and it obeys the diffusion equation (3.11). Thus, to compute the fluid production rate for primary recovery using the macroscopic equation (A 4) given in appendix A, we only need to upscale the density diffusion equation (3.11) subject to the zero-mass-flux condition on the pore wall,

(4.1) $$\begin{eqnarray}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}=0.\end{eqnarray}$$

The spatial deviation of the density, $\unicode[STIX]{x1D70C}-\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ , as defined in Whitaker (Reference Whitaker1986) with $\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ being the intrinsic average of the density (appendix A), also satisfies the boundary condition (4.1).

The problem of upscaling a diffusion equation with a heterogeneous reaction in porous media has been studied in detail by Whitaker (Reference Whitaker1999) and more recently by Lugo-Mendez et al. (Reference Lugo-Mendez, Valdes-Parada, Porter, Wood and Ochoa-Tapia2015). Other works on upscaling of a diffusion equation using the same methodology include, but are not limited to, Wood & Whitaker (Reference Wood and Whitaker1997, Reference Wood and Whitaker1999), Valdes-Parada et al. (Reference Valdes-Parada, Porter, Narayanaswamy, Ford and Wood2009), Valdes-Parada & Aguilar-Madera (Reference Valdes-Parada and Aguilar-Madera2011), Valdes-Parada, Aguilar-Madera & Alvarez-Ramirez (Reference Valdes-Parada, Aguilar-Madera and Alvarez-Ramirez2011), Benitez-Olivares, Valdes-Parada & Saucedo-Castaneda (Reference Benitez-Olivares, Valdes-Parada and Saucedo-Castaneda2016), Santos-Sanchez, Valdes-Parada & Chirino (Reference Santos-Sanchez, Valdes-Parada and Chirino2016) and Valdes-Parada, Lasseux & Whitaker (Reference Valdes-Parada, Lasseux and Whitaker2017). In the absence of reaction, the microscopic-scale diffusion problem studied by Whitaker (Reference Whitaker1999) and Lugo-Mendez et al. (Reference Lugo-Mendez, Valdes-Parada, Porter, Wood and Ochoa-Tapia2015) is exactly the same as the diffusion equation (3.11) with the no-penetration condition on the fluid–solid surface (4.1). Thus, the upscaled density diffusion equation for the present problem is a special case of the upscaled diffusion equation of Lugo-Mendez et al. (Reference Lugo-Mendez, Valdes-Parada, Porter, Wood and Ochoa-Tapia2015, equation (41)), which takes the form

(4.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D70C}\rangle ^{f}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(D_{{<}\unicode[STIX]{x1D70C}{>}^{f}}\unicode[STIX]{x1D735}\langle \unicode[STIX]{x1D70C}\rangle ^{f}),\end{eqnarray}$$

where the diffusion coefficient is valued at the averaged density,

(4.3) $$\begin{eqnarray}D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}=\frac{\unicode[STIX]{x1D707}_{b}+4\unicode[STIX]{x1D707}/3}{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}.\end{eqnarray}$$

The upscaled equation (4.2) is derived under the condition that the closure problem for the spatial deviation of the density is local, periodic, quasi-steady and linear. Under this condition and in the absence of reaction, the local closure problem for the spatial deviation of the density (equations (18a–e) of Lugo-Mendez et al. (Reference Lugo-Mendez, Valdes-Parada, Porter, Wood and Ochoa-Tapia2015)), is homogeneous when subject to the no-penetration boundary condition and zero initial deviation. Thus, the spatial deviation of the density is identically zero. The specific time- and length-scale constraints corresponding to the above condition for simplifying the closure problem have been examined in detail by Lugo-Mendez et al. (Reference Lugo-Mendez, Valdes-Parada, Porter, Wood and Ochoa-Tapia2015), which become (note that the so-called Thiele modulus for the closure problem is zero in the absence of reaction)

(4.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\ell _{\unicode[STIX]{x1D6FE}}\ll r_{0}\ll L,\\ {\displaystyle \frac{\ell _{\unicode[STIX]{x1D6FE}}}{L}}\ll 1,\\ 1\ll {\displaystyle \frac{D_{\unicode[STIX]{x1D70C},0}t^{\ast }}{\ell _{\unicode[STIX]{x1D6FE}}^{2}}}.\end{array}\right\}\end{eqnarray}$$

In (4.4), $\ell _{\unicode[STIX]{x1D6FE}}$ is the characteristic length for the fluid phase, which can be taken as the pore throat diameter; $r_{0}$ is the characteristic size of the averaging volume; $L$ is the characteristic length of the porous system, which can be taken as the size of the reservoir; and $t^{\ast }$ is a characteristic time for the fluid production. The characteristic time for primary recovery is the diffusive time scale $L^{2}/D_{\unicode[STIX]{x1D70C},0}$ ; thus, the third condition in (4.4) coincides with the second condition. The diffusion coefficient (4.3) is valued at the mean density $\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ of the averaging cell. This can be obtained by a Taylor-series expansion and neglecting the nonlinear terms under the local theory, as shown by Lugo-Mendez et al. (Reference Lugo-Mendez, Valdes-Parada, Porter, Wood and Ochoa-Tapia2015). In cases where the length-scale constraints prevent such a linearization, Lugo-Mendez et al. (Reference Lugo-Mendez, Valdes-Parada, Porter, Wood and Ochoa-Tapia2015) have also developed an iterative procedure for the nonlinear closure problem.

The new diffusion equation (4.2) for $\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ can be compared to the classical diffusion equation (A 6) for $\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ based on Darcy’s law listed in appendix A. The diffusion coefficient $D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}$ for (4.2) differs significantly, both physically and magnitude-wise, from the classical diffusion coefficient $D_{Darcy}$ in (A 6) given by (A 7), $D_{Darcy}=\unicode[STIX]{x1D705}c_{T}^{2}\langle \unicode[STIX]{x1D70C}\rangle ^{f}/(\unicode[STIX]{x1D707}\unicode[STIX]{x1D719})$ , with $\unicode[STIX]{x1D705},~\unicode[STIX]{x1D719},~c_{T}$ being the permeability and porosity of the porous medium and the isothermal speed of sound of the fluid, respectively. The diffusion coefficient $D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}$ is directly related to the self-diffusion coefficient of the fluid,

(4.5) $$\begin{eqnarray}D_{11}=(\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}+4/3)\unicode[STIX]{x1D708}.\end{eqnarray}$$

When the bulk viscosity is zero, $D_{11}=1.33\unicode[STIX]{x1D708}$ for a dilute gas (Hirschfelder et al. Reference Hirschfelder, Curtiss and Bird1954). The appearance of the self-diffusion coefficient in the pore-scale fluid density diffusion equation (3.11) as well as the upscaled fluid density diffusion equation (4.2) for the intrinsic average density $\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ is natural and expected, as the density diffusion equation describes the mass transport process in the diffusion-dominated flow regime. On the other hand, $D_{Darcy}$ is inversely proportional to the self-diffusion coefficient $D_{11}$ , since $D_{Darcy}$ is inversely proportional to the kinematic viscosity $\unicode[STIX]{x1D708}$ . Darcy’s law is only appropriate for incompressible flows and it is irrelevant to the self-diffusion process. Therefore, the classical diffusion coefficient $D_{Darcy}$ cannot be the correct diffusion coefficient for the intrinsic average density $\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ .

In addition, the magnitude of $D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}$ can be several orders of magnitude higher than $D_{Darcy}$ . If we evaluate the diffusivities $D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}$ and $D_{Darcy}$ for methane at typical shale conditions at temperature $80\,^{\circ }\text{C}$ and pressure of 25 MPa (appendix B), we have $D_{Darcy}=4.53\times 10^{-8}~\text{m}^{2}~\text{s}^{-1}$ when the porosity is $\unicode[STIX]{x1D719}=0.1$ and the permeability is one nanodarcy ( $\unicode[STIX]{x1D705}\approx 10^{-21}~\text{m}^{2}$ ); while $D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}=4.82\times 10^{-5}~\text{m}^{2}~\text{s}^{-1}$ . Thus, for methane in shale, $D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}$ is three orders of magnitude higher than $D_{Darcy}$ . At a lower pressure, the difference between these two diffusion coefficients becomes even larger, since, as the pressure is decreased, $D_{\langle \unicode[STIX]{x1D70C}\rangle ^{f}}$ increases while $D_{Darcy}$ decreases.

To find the mass flux $\boldsymbol{J}$ , which is the mass flow rate per unit porous medium area, the density diffusion equation (4.2) needs to be written in the form

(4.6) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}m}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}=0,\end{eqnarray}$$

where $m$ is fluid mass per unit porous medium volume, i.e. the superficial average density $\langle \unicode[STIX]{x1D70C}\rangle$ , which is related to intrinsic averaged density $\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ by $\langle \unicode[STIX]{x1D70C}\rangle =\unicode[STIX]{x1D719}\langle \unicode[STIX]{x1D70C}\rangle ^{f}$ (appendix A). Thus,

(4.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D70C}\rangle }{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\unicode[STIX]{x1D719}\frac{\unicode[STIX]{x1D707}_{b}+4\unicode[STIX]{x1D707}/3}{\langle \unicode[STIX]{x1D70C}\rangle }\unicode[STIX]{x1D735}\langle \unicode[STIX]{x1D70C}\rangle \right)=0.\end{eqnarray}$$

Equation (4.7) is the effective (i.e. homogenized) medium equation, which shows that, in the effective medium, mass (density) diffuses with a mass diffusion coefficient

(4.8) $$\begin{eqnarray}D_{m}=\unicode[STIX]{x1D719}\frac{\unicode[STIX]{x1D707}_{b}+4\unicode[STIX]{x1D707}/3}{\langle \unicode[STIX]{x1D70C}\rangle }.\end{eqnarray}$$

Clearly, $D_{m}$ depends on the porosity, but not on the permeability. The mass flux $\boldsymbol{J}$ is given by

(4.9) $$\begin{eqnarray}\boldsymbol{J}=-D_{m}\unicode[STIX]{x1D735}\langle \unicode[STIX]{x1D70C}\rangle ,\end{eqnarray}$$

which is also independent of permeability. Permeability is a property introduced by Darcy’s law for incompressible flows and it plays no role in determining the mass flow rate in primary recovery. Thus, Darcy’s law should be replaced by the mass flux equation (4.9) for primary fluid recovery.

Figure 3. Top: Microscale simulation for density in drainage flow (primary recovery) with 100 horizontally repeated identical unit cells. The solid circles (totalling 49 for $\unicode[STIX]{x1D719}=0.5840$ ) are randomly distributed in the unit cell. Bottom: Macroscale effective medium simulation for density distribution. Three porosities $\unicode[STIX]{x1D719}=0.5840,~0.4630,~0.3202$ are considered.

5 Comparison of macroscale equation with pore-scale simulation

In this section, we solve the upscaled and the microscale density diffusion equations for drainage of a compressible fluid from a semi-sealed two-dimensional porous medium channel and compare the average density distributions (figure 3). The porous medium is made of 100 identical unit cells with a compressible fluid saturating the pore space with an initial density $\unicode[STIX]{x1D70C}_{i}$ . The solids in a unit cell with a size of 0.02 m $\times$ 0.02 m are randomly distributed circles with various diameters between 0.0004 m and 0.0025 m obeying the normal size distribution (several such random generated patterns have been tested and they produce the same cell averaged densities). The left end of the channel is sealed while the right end is opened at $t=0$ , with an exit density $\unicode[STIX]{x1D70C}_{e}<\unicode[STIX]{x1D70C}_{i}$ . We solve the microscale density diffusion equation (3.11) in the pore space, together with the no-penetration condition (4.1) imposed on the solid surfaces, including the channel walls. We then compute the average density in each unit cell and compare the result to the corresponding macroscale density. The macroscale density is obtained by solving the effective medium density diffusion equation (4.7) subject to the no-penetration condition on the channel walls. Both calculations are carried out using COMSOL (COMSOL 2016). In these computations, $\unicode[STIX]{x1D70C}_{i}=1~\text{kg}~\text{m}^{-3}$ ,  $\unicode[STIX]{x1D70C}_{e}=0.9~\text{kg}~\text{m}^{-3}$ ,  $\unicode[STIX]{x1D707}=10^{-3}~\text{Pa}~\text{s}$ ,  $r_{0}=0.02~\text{m}$ , $L=2~\text{m}$ and three porosities $\unicode[STIX]{x1D719}=0.5840,0.4630,0.3202$ have been used. A large viscosity value is chosen in order to shorten the computational time required for draining out the fluid from the channel. Convergence tests have been conducted by mesh refinements.

Figure 4 shows the local density distributions for each system for the case $\unicode[STIX]{x1D719}=0.5840$ in three segments of the channel at time $t=1000$  s, for cells 6–15, cells 46–55 and cells 85–94. It is observed that the microscale density distribution in each unit cell is fairly uniform despite the presence of the solid fractions, and there is little difference between the microscale and macroscale densities in these three segments. Figure 5 shows a specific comparison of the cell-averaged dimensionless density $\unicode[STIX]{x1D70C}_{avg}/\unicode[STIX]{x1D70C}_{i}$ for these two systems for cells 10, 50 and 90 (marked as A, B and C, respectively). For the microscale model, the average density $\unicode[STIX]{x1D70C}_{avg,mic}$ is the intrinsically averaged density in the unit cell; while for the macroscale model, $\unicode[STIX]{x1D70C}_{avg,mac}$ is averaged over the entire unit cell area since the medium is homogenized. The relative error of the cell-averaged density between these two models, defined as $|\unicode[STIX]{x1D70C}_{avg,mac}-\unicode[STIX]{x1D70C}_{avg,mic}|/\unicode[STIX]{x1D70C}_{avg,mic}$ , is plotted against time for each of these three cells in figure 6(a). The maximum relative errors for these three cells (A, B and C) are 0.6 %, 0.4 % and 0.1 %, respectively. Similarly, for $\unicode[STIX]{x1D719}=0.4630$ , the maximum relative errors for these three cells are 0.65 %, 0.45 % and 0.15 %, respectively, as shown in figure 6(b). For $\unicode[STIX]{x1D719}=0.3202$ , the maximum relative errors for these three cells are 1 %, 0.7 % and 0.3 %, respectively (figure 6 c). This range of relative error is acceptable for many practical applications. Together, these comparisons provide a numerical validation of the upscaled density equation.

Figure 4. Computed local density distributions in the various cells indicated for the two systems, microscale (top) and macroscale (effective medium) (bottom) at time $t=1000~\text{s}$ , for $\unicode[STIX]{x1D719}=0.5840$ and $\unicode[STIX]{x1D707}=10^{-3}~\text{Pa}~\text{s}$ . Initial density $\unicode[STIX]{x1D70C}_{i}=1~\text{kg}~\text{m}^{-3}$ and exit density $\unicode[STIX]{x1D70C}_{e}=0.9~\text{kg}~\text{m}^{-3}$ .

Figure 5. Comparison of the dimensionless average density for the three cells (A, B and C) computed from the macroscale equation and the microscale equation, for $\unicode[STIX]{x1D719}=0.5840$ .

Figure 6. Relative error $|\unicode[STIX]{x1D70C}_{avg,mac}-\unicode[STIX]{x1D70C}_{avg,mic}|/\unicode[STIX]{x1D70C}_{avg,mic}$ versus time for the three representative cells: (a) $\unicode[STIX]{x1D719}=0.5840$ , (b) $\unicode[STIX]{x1D719}=\text{0.4630}$ and (c) $\unicode[STIX]{x1D719}=\text{0.3202}$ .