The Darcy law gives

$u_\alpha= -\frac1{\mu_\alpha} K_\alpha\nabla p_\alpha$

where $\mu_\alpha$ and $K_\alpha$ represent the viscosity and the permeability tensor for each phase $\alpha$ . In the two phase case, we have either $\alpha=w$ or $\alpha=o$ . In this tutorial, we do not take into account capillary pressure so that $p=p_w=p_o$ and gravity effects. We denote by $K$ the absolute permeability tensor and each phase permeability is defined through its relative permeability by the expression

$K_\alpha=k_{r\alpha}K.$

The phase mobility are defined as

$\lambda_\alpha=\frac{k_{r\alpha}}{\mu_\alpha}$

so that the Darcy law may be rewritten as

$u_\alpha= -\lambda_\alpha K\nabla p.$

The conservation of mass for each phase writes:

$\frac{\partial}{\partial t}(\phi\rho_\alpha s_\alpha)+\nabla\cdot(\rho_\alpha u_\alpha)=q_\alpha$

where $s_\alpha$ denotes the saturation of the phase $\alpha$ and $q_\alpha$ is a source term. Let us consider a two phase flow with oil and water. We assume that the rock and both fluid phases are incompressible. Since $s_w+s_o=1$ , we may add the conservation equations to get

$\nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o}.$

where we define

$u=u_w+u_o.$

Let the total mobility be equal to

$\lambda=\lambda_w+\lambda_o$

Then, we have

$u=-\lambda K\nabla p.$

The set of equations

$\nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o},\quad u=-\lambda K\nabla p.$

is referred to as the pressure equation. We introduce the fractional flow $f_w$ as

$f_w=\frac{\lambda_w}{\lambda_w+\lambda_o}$

and obtain

$\phi\frac{\partial s_w}{\partial t}+\nabla\cdot(f_w u)=\frac{q_w}{\rho_w}$

which is referred to as the transport equation. The pressure and transport equation are coupled. In this tutorial, we implement a splitting scheme, where, at each time step, we decouple the two equations. We solve first the pressure equation and then update the water saturation by solving the transport equation assuming that $u$ is constant in time in the time step interval we are considering.

Program walk-through.

Main function

int main ()
try
{