PETRONODE - Iterative Density Porosity

Iterative Density Porosity

(Python 3.x)

Algorithm for computing total porosity (Phi_t) in formations with variable matrix density.

The "classic" density porosity algorithm presumes that matrix density, \(\rho _{\mathit{ma}}\), remains constant within the lithofacies:

\[\phi _t=\frac{\rho _{\mathit{ma}}-\rho _b}{\rho _{\mathit{ma}}-\rho _f},\rho _{\mathit{ma}}=\mathit{Const}\]

In many cases, however, such matrix density can be a function of total porosity:

\[\phi _t=\frac{\rho _{\mathit{ma}}(\phi _t)-\rho _b}{\rho _{\mathit{ma}}(\phi _t)-\rho _f}\]

This algorithm resolves the above equation numerically. It must be stressed that the dependency between the matrix density and total porosity is empirical and relies on laboratory measurements performed on full-bore or side-wall cores.

Download source code: Iterative_Density_Porosity.py

import numpy as np


def Iterative_Density_Porosity( Rho_b, Rho_fl, Rho_ma=lambda p: 2.65,
                                accuracy=0.001, max_iter=100,
                                rho_ma_initial=3.0, verbose=False):
    """
    Solves the density iteratively into total porosity
    Rho_b - array-like, bulk density, in g/cc
    Rho_fl - array-like, density of fluid in flushed zone, in g/cc
    Rho_ma - equation for matrix density, dependency from total porosity Phi_t
    accuracy = 0.001 - desired solution accuracy, in v/v
    max_iter = 100 - maximum number of iterations
    rho_ma_initial = 3.0 - initial matrix density, g/cc
    verbose - print intermediate results
    returns:
        total porosity Phi_t as an array, in v/v
        estimated matrix density as an array, in g/cc
    """
    density = np.array(Rho_b)
    l = len(density)
    fluid = np.array(Rho_fl)
    if len(fluid) != l:
        fluid = np.ones(l) * fluid[0]
    Phi_t = np.clip( (rho_ma_initial-density) / (rho_ma_initial-fluid), 0, 1)
    Rho_ma_est = np.ones(l) * rho_ma_initial
    for i, d in enumerate(density):
        f = fluid[i]
        p = Phi_t[i]
        for j in range(max_iter):
            Rho_ma_est[i] = Rho_ma(p)
            Phi_t[i] = np.clip((Rho_ma_est[i] - d) / (Rho_ma_est[i] - f), 0, 1)
            if abs( Phi_t[i] - p) >= accuracy:
                p = Phi_t[i]
                continue    
            if verbose:
                print("Rho_b={:.3f}, Rho_f={:.3f} : Phi_t={:.3f}, Rho_ma={:.3f}"
                    " (after {:d} iterations)".format( d, f, Phi_t[i], Rho_ma_est[i], j))
            break
    return Phi_t, Rho_ma_est

if __name__ == "__main__":
    print("Polymictic feldspar sandstone with matrix density (from core)")
    print("\t2.56 + Phi_t*0.55:")
    print("Flushed zone fluid density is assumed 1.03 g/cc;")
    print("Bulk density is environmentally corrected;")    
    print("Rho_b\t\tPhi_t(classic)\tPhi_t(iterative)")
    Rho_b = np.linspace(2.0, 2.6, 12)
    Phi_t_classic = (2.65-Rho_b) / (2.65-1)
    Phi_t, Rho_ma_apparent = Iterative_Density_Porosity( Rho_b, [1.03],
        Rho_ma=lambda p: np.clip(2.56+p*0.55, 2.56, 2.67), verbose=False)
    for rhob, phit_c, phit in zip(Rho_b, Phi_t_classic, Phi_t):
        print("{:.3f}\t\t{:.3f}\t\t{:.3f}".format(rhob, phit_c, phit))


Home Software Manuals Training Codelets Iterative Density Porosity Archie S_wa Triangle Solver Height by Triangulation Object Size by Triangulation Nelder-Mead Optimization Letters and numbers Mission Papers Links Contacts	Iterative Density Porosity (Python 3.x) Algorithm for computing total porosity (Phi_t) in formations with variable matrix density. The "classic" density porosity algorithm presumes that matrix density, \(\rho _{\mathit{ma}}\), remains constant within the lithofacies: \[\phi _t=\frac{\rho _{\mathit{ma}}-\rho _b}{\rho _{\mathit{ma}}-\rho _f},\rho _{\mathit{ma}}=\mathit{Const}\] In many cases, however, such matrix density can be a function of total porosity: \[\phi _t=\frac{\rho _{\mathit{ma}}(\phi _t)-\rho _b}{\rho _{\mathit{ma}}(\phi _t)-\rho _f}\] This algorithm resolves the above equation numerically. It must be stressed that the dependency between the matrix density and total porosity is empirical and relies on laboratory measurements performed on full-bore or side-wall cores. Download source code: Iterative_Density_Porosity.py import numpy as np def Iterative_Density_Porosity( Rho_b, Rho_fl, Rho_ma=lambda p: 2.65, accuracy=0.001, max_iter=100, rho_ma_initial=3.0, verbose=False): """ Solves the density iteratively into total porosity Rho_b - array-like, bulk density, in g/cc Rho_fl - array-like, density of fluid in flushed zone, in g/cc Rho_ma - equation for matrix density, dependency from total porosity Phi_t accuracy = 0.001 - desired solution accuracy, in v/v max_iter = 100 - maximum number of iterations rho_ma_initial = 3.0 - initial matrix density, g/cc verbose - print intermediate results returns: total porosity Phi_t as an array, in v/v estimated matrix density as an array, in g/cc """ density = np.array(Rho_b) l = len(density) fluid = np.array(Rho_fl) if len(fluid) != l: fluid = np.ones(l) * fluid[0] Phi_t = np.clip( (rho_ma_initial-density) / (rho_ma_initial-fluid), 0, 1) Rho_ma_est = np.ones(l) * rho_ma_initial for i, d in enumerate(density): f = fluid[i] p = Phi_t[i] for j in range(max_iter): Rho_ma_est[i] = Rho_ma(p) Phi_t[i] = np.clip((Rho_ma_est[i] - d) / (Rho_ma_est[i] - f), 0, 1) if abs( Phi_t[i] - p) >= accuracy: p = Phi_t[i] continue if verbose: print("Rho_b={:.3f}, Rho_f={:.3f} : Phi_t={:.3f}, Rho_ma={:.3f}" " (after {:d} iterations)".format( d, f, Phi_t[i], Rho_ma_est[i], j)) break return Phi_t, Rho_ma_est if __name__ == "__main__": print("Polymictic feldspar sandstone with matrix density (from core)") print("\t2.56 + Phi_t0.55:") print("Flushed zone fluid density is assumed 1.03 g/cc;") print("Bulk density is environmentally corrected;") print("Rho_b\t\tPhi_t(classic)\tPhi_t(iterative)") Rho_b = np.linspace(2.0, 2.6, 12) Phi_t_classic = (2.65-Rho_b) / (2.65-1) Phi_t, Rho_ma_apparent = Iterative_Density_Porosity( Rho_b, [1.03], Rho_ma=lambda p: np.clip(2.56+p0.55, 2.56, 2.67), verbose=False) for rhob, phit_c, phit in zip(Rho_b, Phi_t_classic, Phi_t): print("{:.3f}\t\t{:.3f}\t\t{:.3f}".format(rhob, phit_c, phit))
(C) M.Yakimov, 2006	Contact us