View Source

We start with the reservoir flow continuity equation:

\frac{\partial (\rho \phi)}{\partial t} + \nabla \, ( \rho \, {\bf u}) = \sum_k \dot m_k(t) \cdot \delta({\bf r}-{\bf r}_k)

percolation model:

{\bf u} = - M \cdot ( \nabla p - \rho \, {\bf g})

and the reservoir boundary flow condition:

{\rm F}_{\Gamma}(p, {\bf u}) = 0

where

	well-reservoir contact of the -th well
	normal vector of differential area on the well-reservoir contact, pointing inside wellbore
	mass flowrate at -th well
	sandface flowrate at -th well
	fluid density as function of reservoir fluid pressure

Then use the following equality:

d(\rho \, \phi) = \rho \, d \phi + \phi \, d\rho = \rho \, \phi \, \left( \frac{d \phi }{\phi} +  \frac{d \rho }{\rho}  \right) 
= \rho \, \phi \, \left( \frac{1}{\phi} \frac{d \phi}{dp} \, dp +  \frac{1}{\rho} \frac{d \rho}{dp} \, dp  \right) 
= \rho \, \phi \, (c_{\phi} \, dp + c \, dp) = \rho \, \phi \, c_t \, dp

where

total compressibility

to arrive at:

\rho \, \phi \, c_t  \cdot \frac{\partial p}{\partial t} + \nabla \, ( \rho \, {\bf u}) = \sum_k \rho \, q_k(t) \cdot \delta({\bf r}-{\bf r}_k)

{\rm F}_{\Gamma}(p, {\bf u}) = 0

The left-hand side of equation can be transformed in the following way:

\nabla \, ( \rho \, {\bf u}) =  \rho \, \nabla \, {\bf u} + (\nabla  \rho, \, {\bf u}) =  \rho \, \nabla \, {\bf u} + \frac{d\rho}{dp} \cdot (\nabla  p, \, {\bf u}) = \rho \, \nabla \, {\bf u} + \rho \, c \cdot (\nabla  p, \, {\bf u})

where is fluid compressibility.

By using the Dirac delta function property: the right-hand side of equation can be transformed in the following way:

\sum_k  \rho(p(t, {\bf r}))  \cdot q_k(t) \cdot \delta({\bf r}-{\bf r}_k) 
= \sum_k  \rho(p(t, {\bf r}_k)) \cdot  q_k(t) \cdot \delta({\bf r}-{\bf r}_k) 
 = \sum_k  \rho(p(t, {\bf r})) \cdot  q_k(t)  \cdot \delta({\bf r}-{\bf r}_k)  
= \rho(p) \cdot \sum_k q_k(t) \cdot \delta({\bf r}-{\bf r}_k)

Substituting and into and reducing the density one arrives to:

 \phi \, c_t  \cdot \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)

{\rm F}_{\Gamma}(p, {\bf u}) = 0