Finite-difference formulation

The equations that govern the fluid flow through a porous medium can be solved numerically using a finite-difference method. For three dimensions, the finite difference approximations of equation (12) through (14) can be written as

${V_{ijk}\over \Delta t} \Delta_t \: (\phi \rho_g S_g + \phi \rho_w (1-S_g-S_h) ... ...{ijk} = [ (\Delta T_g \rho_g \Delta P)_i - (\Delta T_g {\rho_g}^2 g \Delta z)_i$

$+ (\Delta T_w \rho_w X_{g,w} \Delta P)_i - (\Delta T_w {\rho_w}^2 X_{g,w} g \Delta z)_i] + [ (\Delta T_g \rho_g \Delta P)_j- (\Delta T_g {\rho_g}^2 g \Delta z)_j$

$+ (\Delta T_w \rho_w X_{g,w} \Delta P)_j - (\Delta T_w {\rho_w}^2 X_{g,w} g \De... ...z)_j] + [ (\Delta T_g \rho_g \Delta P)_k - (\Delta T_g {\rho_g}^2 g \Delta z)_k$
 

$$+ (\Delta T_w \rho_w X_{g,w} \Delta P)_k - (\Delta T_w {\rho_w}^2 X_{g,w} g \Delta z)_k] + \dot{m_g}_{ijk},$$ (15)

${V_{ijk}\over \Delta t} \Delta_t \: (\phi \rho_w (1-S_g-S_h) (1-X_{g,w}) + \phi \rho_h S_h (1-X_{g,h}))_{ijk} = [ (\Delta T_w \rho_w (1-X_{g,w}) \Delta P)_i$

$- (\Delta T_w {\rho_w}^2 (1-X_{g,w}) g \Delta z)_i] + [ (\Delta T_w \rho_w (1-X_{g,w}) \Delta P)_j - (\Delta T_w {\rho_w}^2 (1-X_{g,w}) g \Delta z)_j]$
 

$$+ [(\Delta T_w \rho_w (1-X_{g,w}) \Delta P)_k - (\Delta T_w {\rho_w}^2 (1-X_{g,w}) g \Delta z)_k] + \dot{m_w}_{ijk},$$ (16)

 

$${V_{ijk}\over \Delta t} \Delta_t \: (\phi \rho_h S_h)_{ijk} = - \dot{m_h}_{ijk}$$ (17)

where V_ijk is the grid block volume and

$$T_g = {A \over \Delta d} \kappa_{abs} {\kappa_{rg} \over \mu_g},$$ (18)

$$T_w = {A \over \Delta d} \kappa_{abs} {\kappa_{rw} \over \mu_w}.$$ (19)

are the gas and water transmissibilities. The cross-sectional area normal to the direction of the flow at the block-boundary is denoted by A, and the distance between two connecting grid points by $\Delta d$.

The accumulation terms in equations (15) to (17) are given as

$\Delta_t \: (\phi \rho_g S_g + \phi \rho_w (1-S_g-S_h) X_{g,w}+ \phi \rho_h S_h X_{g,h})_{ijk}$

$= (\phi \rho_g S_g + \phi \rho_w (1-S_g-S_h) X_{g,w} + \phi \rho_h S_h X_{g,h})_{ijk}^{t+1}$ 

$$- (\phi \rho_g S_g + \phi \rho_w (1-S_g-S_h) X_{g,w}+ \phi \rho_h S_h X_{g,h})_{ijk}^{t},$$ (20)

$\Delta_t \:(\phi \rho_w (1-S_g-S_h) (1-X_{g,w}) + \phi \rho_h S_h (1-X_{g,h}))_{ijk}$

$= (\phi \rho_w (1-S_g-S_h) (1-X_{g,w}) + \phi \rho_h S_h (1-X_{g,h}))_{ijk}^{t+1}$ 

$$- (\phi \rho_w (1-S_g-S_h) (1-X_{g,w}) + \phi \rho_h S_h (1-X_{g,h}))_{ijk}^t,$$ (21)

 

$$\Delta_t \: (\phi \rho_h S_h)_{ijk} = (\phi \rho_h S_h)_{ijk}^{t+1} - (\phi \rho_h S_h)_{ijk}^t.$$ (22)

The finite-difference operators in equations (15) to (17) are given by

 

$$(\Delta T_g \rho_g \Delta P)_l = T_{g_{l+{1 \over 2}}} \rho_... ...T_{g_{l-{1 \over 2}}} \rho_{g_{l-{1 \over 2}}} (P_l - P_{l-1}),$$ (23)

 

$$(\Delta T_g {\rho_g}^2 g \Delta z)_l = T_{g_{l+{1 \over 2}}}... ...l-{1 \over 2}}} {\rho^2_{g_{l-{1 \over 2}}}} g (z_l - z_{l-1}),$$ (24)

 

$$(\Delta T_w \lambda_w \Delta P)_l = T_{w_{l+{1 \over 2}}} \l... ...w_{l-{1 \over 2}}} \lambda_{w_{l-{1 \over 2}}} (P_l - P_{l-1}),$$ (25)

 

$$(\Delta T_w \Gamma_w g \Delta z)_l = T_{w_{l+{1 \over 2}}} \... ...w_{l-{1 \over 2}}} \Gamma_{w_{l-{1 \over 2}}} g (z_l - z_{l-1})$$ (26)

where subscript l = i,j,k and

$\lambda_w = \rho_w X_{g,w}$ for equation (15),
$\lambda_w = \rho_w (1-X_{g,w})$ for equation (16),
$\Gamma_w = {\rho_w}^2 X_{g,w}$ for equation (15) and
$\Gamma_w = {\rho_w}^2 (1-X_{g,w})$ for equation (16).

Furthermore, it is:  

$$T_{g_{l+{1 \over 2}}} = {A_{l+{1 \over 2}} \over \Delta d_l} (\kappa_{abs} {\kappa_{rg} \over \mu_g} )_{l+{1 \over 2}},$$ (27)

 

$$T_{w_{l+{1 \over 2}}} = {A_{l+{1 \over 2}} \over \Delta d_l} (\kappa_{abs} {\kappa_{rw} \over \mu_w} )_{l+{1 \over 2}}$$ (28)

where $A_{l+{1 \over 2}}$ is the cross-sectional area normal to the direction of the flow at the block-boundary $l+{1 \over 2}$. The distance between grid points l and l+1 is given by $\Delta d_j$. The absolute permeability is taken in the direction of the flow.

Equations (20) to (28) use coefficients that must be evaluated at intercell boundaries between two grid points. Simple arithmetic mean is used to approximate the densities $\rho_g$, $\rho_w$, and $\rho_h$, and the mass fractions X_g,w and X_g,h. The absolute permeability $\kappa_{abs}$ as a function of S_h is determined using harmonic weighting. Upstream weighting is used for the relative relative permeabilities and viscosities. Defining the phase potential as

$$\Phi_p = P_p - \rho_p g h$$ (29)

where the subscript p represents either gas phase g or water phase w, upstream weighting means:

$\Phi_{p_l} \gt \Phi_{p_{l-1}}$ $\longrightarrow$ $({\kappa_{rp} \over \mu_p})_{l-{1 \over 2}} = ({\kappa_{rp} \over \mu_p})_l$

$\Phi_{p_l} < \Phi_{p_{l-1}}$ $\longrightarrow$ $({\kappa_{rp} \over \mu_p})_{l-{1 \over 2}} = ({\kappa_{rp} \over \mu_p})_{l-1}$

Equations (15) to (28) represent the finite-difference set-up for a three-dimensional fluid-flow simulator of a gas-water-hydrate system, which can be solved using the Newton-Raphson method.