Consider the Cartesian coordinates in 3D space: ℝ^3 \Big |_{ \{x, y, z \} } and its infinitesimal volumetric element: \delta \Omega = \{ (x, x+\delta), (y, y+\delta y), (z, z+\delta z) \} \in ℝ^3 with volume \delta V = \delta x \, \delta y \, \delta z bounded by six faces: \{ (\delta \Sigma_x, \, \delta \Sigma_{x+\delta x}), \, (\delta \Sigma_y, \, \delta \Sigma_{y+\delta y}), \, (\delta \Sigma_z, \, \delta \Sigma_{z+\delta z}) \} which have the same area along corresponding axis:
(1) | \delta A(\delta \Sigma_x) = \delta A(\Sigma_{x+\delta x }) = \delta A_{yz} = \delta y \cdot \delta z |
(2) | \delta A(\delta \Sigma_y) = \delta A(\Sigma_{y+\delta y }) = \delta A_{xz} = \delta x \cdot \delta z |
(3) | \delta A(\delta \Sigma_z) = \delta A(\Sigma_{z+\delta z }) = \delta A_{xy} = \delta x \cdot \delta y |
Consider the volumetric element \delta \Omega is filled with porous media with porosity \phi(x,y,z) saturated by fluid with density \rho(x,y,z).
The pore volume is going to be \delta V_{\phi} = \phi \cdot \delta V and the fluid mass contained in this volume is \delta m = \rho \cdot \delta V_{\phi} = \rho \cdot \phi \cdot \delta V.
The mass flowrate through any face \delta \Sigma with area \delta A is defined as:
(4) | \frac{dm}{dt} \Big|_{\delta \Sigma} = {\bf j} \, {\bf \delta A} |
where
{\bf \delta A} = \delta A \cdot {\bf n} | vector area |
{\bf n} | normal vector to elementary area \delta A |
{\bf j} = \rho \cdot {\bf u} | mass flux vector |
{\bf u} | fluid flow velocity |
The total mass balance of the volumetric element \delta \Omega honours the mass conservation:
(5) | \frac{dm}{dt} \Big|_{\delta \Omega} = \sum_{\alpha} j_{\alpha}A_{\alpha} + \delta \dot m_q |
(6) | \frac{dm}{dt} \Big|_{\delta \Omega} = j_x|_{x}\cdot \delta A_{yz} - j_x|_{x+\delta x}\cdot \delta A_{yz} +j_y|_{y}\cdot \delta A_{xz} - j_y|_{y+\delta y}\cdot \delta A_{xz} + j_z|_{z}\cdot \delta A_{xy} - j_z|_{z+\delta z}\cdot \delta A_{xy} + \delta \dot m_q |
where
\delta \dot m_q | the rate of the mass variation which happens inside the volumetric element \delta \Omega |
Dividing the (5) by the volume \delta V:
(7) | \frac{dm}{dt \, \delta V} \Big|_{\delta \Omega} = \frac {\partial (\rho \, \phi)}{\partial t} = \frac{j_x|_x - j_x|_{x+\delta x}}{\delta x} + \frac{j_y|_y - j_y|_{y+\delta y}}{\delta y} + \frac{j_z|_z - j_z|_{z+\delta z}}{\delta z} + \frac{\delta \dot m_q}{\delta V} |
or in differential form:
(8) | \frac{\partial (\rho \phi)}{\partial t} = - \nabla \, {\bf j} + \frac{\delta \dot m_q}{\delta V} |
(9) | \frac{\partial (\rho \phi)}{\partial t} + \nabla \, {\bf j} = \frac{\delta \dot m_q}{\delta V} |
The mass rate generated/consumed by a finite number of well-reservoir contacts can be expressed as:
(10) | \frac{\delta \dot m_q}{\delta V} = \sum_k \rho_k \cdot q_k(t) \cdot \delta({\bf r}-{\bf r}_k) |
where
q_k(t) | volumetric flowrate of the source/stock at the reservoir point {\bf r}_k |
\rho_k (t) = \rho(p(t, {\rm r}_k)) | fluid density at the reservoir point {\bf r}_k |
The next step is to re-right
(10) in equivalent form:
(11) | \frac{\delta \dot m_q}{\delta V} = \sum_k \rho_k \cdot q_k(t) \cdot \delta({\bf r}-{\bf r}_k) = \rho(t,{\rm r}) \cdot \sum_k q_k(t) \cdot \delta({\bf r}-{\bf r}_k) |
which turns (9) into:
(12) | \frac{\partial (\rho \phi)}{\partial t} + \nabla \, {\bf j} = \rho \cdot \sum_k q_k(t) \cdot \delta({\bf r}-{\bf r}_k) |
Substituting the mass flux {\bf j} = \rho \cdot {\bf u} into (12):
(13) | \frac{\partial (\rho \phi)}{\partial t} + \nabla \, ( \rho \, {\bf u}) = \rho \cdot \sum_k q_k(t) \cdot \delta({\bf r}-{\bf r}_k) |
Take into account that:
(14) | \frac{\partial (\rho \phi)}{\partial t} = \rho \cdot \frac{\partial \phi}{\partial t} + \frac{\partial \rho}{\partial t} \cdot \phi = \rho \cdot \phi \left( \frac{1}{\rho} \cdot \frac{\partial \phi}{\partial t} + \frac{1}{\phi} \cdot \frac{\partial \rho}{\partial t} \right) = \rho \cdot \phi \cdot ( c + c_r) = \rho \cdot \phi \cdot c_t |
and
(15) | \nabla \, ( \rho \, {\bf u}) = \rho \cdot \nabla \,{\bf u} + \nabla \, \rho \cdot {\bf u} = \rho \cdot \nabla \,{\bf u} + \rho \cdot c \cdot {\bf u} \, \nabla \, p |
where
\displaystyle c_r = \frac{1}{\phi} \, \frac{\partial \phi}{\partial p} | reservoir pore compressibility |
\displaystyle c = \frac{1}{\rho} \, \frac{\partial \rho}{\partial p} | fluid compressibility |
c_t = c_r + c | total compressibility |
Substituting the (14) and (15) in (13) and cancelling the fluid density \rho one arrives to:
(16) | \phi \, c_t \, \frac{\partial p}{\partial t} + \nabla \, {\bf u} + c \cdot {\bf u} \, \nabla \, p = \sum_k q_k(t) \cdot \delta({\bf r}-{\bf r}_k) |
See also
Physics / Mechanics / Continuum mechanics / Fluid Mechanics / Fluid Dynamics / Pressure Diffusion / Pressure Diffusion @model / Single-phase pressure diffusion @model