Geophysics: compute the permeability of a 3D Porous Medium¶

Main author: Wim Degruyter

This tutorial illustrates, step-by-step, how to compute the permeability of a given porous media. This process consists of three major steps:

Read the geometry, defined by a stack of binary (black and white) images. For this tutorial, we use an artificial geometry, consisting of two hemispheres which partially overlap. This represents a simple media with just a single tube trhough which the fluid flows from the inlet to the outlet. It is however esaily replaced by any complex media, by substituting the bitmap images with user-supplied data. A typical example would be a porous media, with a structure obtained from experimental data.
Simulate a stationary (time-independent), pressure-driven flow through this media, by imposing a constant pressure at the inlet, and a constant, lower pressure at the outlet. All physical quantities (velocity, pressure, and viscosity) are, for convenience, expressed in a system of lattice units. The final quantity of interest, the permability, has dimensions of length squared. Therefore, the actual permeability is the lattice permeability times the spatial resolution squared. If you would like to convert other physical quantities into different units, you should read the tutorial on unit conversion.
Compute the permeability which, according to Darcy’s law, is proportional to the ratio between the flow rate through the media and the applied pressure gradient. To be able to apply Darcy’s law one has to make sure the flow is laminar. Therefore it is recommended to run each simulation several times at differnt pressure differnces. The permeability should stay a constant.

The code developed in this tutorial can be found in the directory palabos/examples/tutorial/permeability. This Palabos code has been developed as part of a permeability study of vesicular volcanic rocks [1].

[1]	Degruyter W., Bachmann O., Burgisser A., Malaspinas O. A synchrotron perspective on gas flow through pumices. submitted to Geosphere

Pre-processing¶

Requirements¶

We assume you already have downloaded and installed Palabos. Installation support is provided in the user’s guide. We also recommend to go through the other Palabos tutorials to get acquainted with compiling and running of Palabos programs.

The input required for permeability are a series of black and white images in .bmp format, with filename namexxxx.bmp with xxxx numbers starting from 0001. In the directory twoSpheres one can find a series of bw images defining two connected hemispheres. To inspect the stack of images one can use e.g. ImageJ. Here are two example slices, inlet (left) and halfway (right):

Convert image stack to .dat file¶

The first step is to create an input file from the images, so the code is able to read in the geometry. This conversion is done by a Matlab script called createDAT.m. Open Matlab and go to the directory of createDAT.m. At the prompt type createDAT(number-of-files, path/to/inputprefix, path/to/output.dat) e.g. try for the twoSpheres files and type createDAT(48, 'twoSpheres/', 'twoSpheres', 'twoSpheres.dat'). Every voxel is given a value: 0 for a fluid voxel (blue), 1 for a material voxel that touches (26-connected) a pore voxel (green), and 2 for an interior material voxel (red) illustrated again by the inlet slice (left) and the halfway slice (right):

The Matlab script visualizes the conversion fo illustration purposes. If you are converting a large number of files this process can take some time and it is better to switch of the display by commenting out lines 113-115, 258-260, and 372-374.

Once the .dat file is created we can start a simulation.

Simulation¶

Brief outline of the code¶

The permeability.cpp code is listed below:

#include "palabos3D.h"
#include "palabos3D.hh"
#include <vector>
#include <cmath>

using namespace plb;
using namespace std;


typedef double T;
#define DESCRIPTOR descriptors::D3Q19Descriptor


// This function object returns a zero velocity, and a pressure which decreases
//   linearly in x-direction. It is used to initialize the particle populations.
class PressureGradient {
public:
    PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_)
    { }
    void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const
    {
        velocity.resetToZero();
        density = 1. - deltaP*DESCRIPTOR<T>::invCs2 / (T)(nx-1) * (T)iX;

    }
private:
    T deltaP;
    plint nx;
};

void porousMediaSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
    OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition,
    MultiScalarField3D<int>& geometry, T deltaP)
{
    const plint nx = lattice.getNx();
    const plint ny = lattice.getNy();
    const plint nz = lattice.getNz();

    pcout << "Definition of inlet/outlet." << endl;
    Box3D inlet (0,0,      1,ny-2, 1,nz-2);
    boundaryCondition->addPressureBoundary0N(inlet, lattice);
    setBoundaryDensity(lattice, inlet, 1.);

    Box3D outlet(nx-1,nx-1, 1,ny-2, 1,nz-2);
    boundaryCondition->addPressureBoundary0P(outlet, lattice);
    setBoundaryDensity(lattice, outlet, 1. - deltaP*DESCRIPTOR<T>::invCs2);

    pcout << "Definition of the geometry." << endl;
    // Where "geometry" evaluates to 1, use bounce-back.
    defineDynamics(lattice, geometry, new BounceBack<T,DESCRIPTOR>(), 1);
    // Where "geometry" evaluates to 2, use no-dynamics (which does nothing).
    defineDynamics(lattice, geometry, new NoDynamics<T,DESCRIPTOR>(), 2);

    pcout << "Initilization of rho and u." << endl;
    initializeAtEquilibrium( lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx) );

    lattice.initialize();
    delete boundaryCondition;
}

void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
    const plint nx = lattice.getNx();
    const plint ny = lattice.getNy();
    const plint nz = lattice.getNz();

    const plint imSize = 600;
    ImageWriter<T> imageWriter("leeloo");

    // Write velocity-norm at x=0.
    imageWriter.writeScaledGif( createFileName("ux_inlet", iter, 6),
        *computeVelocityNorm(lattice, Box3D(0,0, 0,ny-1, 0,nz-1)),
        imSize, imSize );

    // Write velocity-norm at x=nx/2.
    imageWriter.writeScaledGif( createFileName("ux_half", iter, 6),
        *computeVelocityNorm(lattice, Box3D(nx/2,nx/2, 0,ny-1, 0,nz-1)),
        imSize, imSize );

}


void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
        VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), 1.);
        vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", 1.);
        vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", 1.);
}

T computePermeability (
    MultiBlockLattice3D<T,DESCRIPTOR>& lattice, T nu, T deltaP, Box3D domain )
{
    pcout << "Computing the permeability." << endl;

    // Compute only the x-direction of the velocity (direction of the flow).
    plint xComponent = 0;
    plint nx = lattice.getNx();

    T meanU = computeAverage (
        *computeVelocityComponent (lattice, domain, xComponent ) );


    pcout << "Average velocity     = " << meanU             << endl;
    pcout << "Lattice viscosity nu = " << nu                << endl;
    pcout << "Grad P               = " << deltaP/(T)(nx-1)            << endl;
    pcout << "Permeability         = " << nu*meanU / (deltaP/(T)(nx-1)) << endl;

    return meanU;
}

int main(int argc, char **argv)
{
    plbInit(&argc, &argv);

    if (argc!=7)
    {
        pcout << "Error missing some input parameter\n";
        pcout << "The structure is :\n";
        pcout << "1. Input file name.\n";
        pcout << "2. Output directory name.\n";
        pcout << "3. number of cells in X direction.\n";
        pcout << "4. number of cells in Y direction.\n";
        pcout << "5. number of cells in Z direction.\n";
        pcout << "6. Delta P .\n";
        exit (EXIT_FAILURE);
    }
    std::string fNameIn  = argv[1];
    std::string fNameOut = argv[2];

    const plint nx = atoi(argv[3]);
    const plint ny = atoi(argv[4]);
    const plint nz = atoi(argv[5]);
    const T   deltaP = atof(argv[6]);

    global::directories().setOutputDir(fNameOut+"/");

    const T omega = 1.0;
    const T nu    = ((T)1/omega-0.5)/DESCRIPTOR<T>::invCs2;

    pcout << "Creation of the lattice." << endl;
    MultiBlockLattice3D<T,DESCRIPTOR> lattice(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega));
    // Switch off periodicity.
    lattice.periodicity().toggleAll(false);

    pcout << "Reading the geometry file." << endl;
    MultiScalarField3D<int> geometry(nx,ny,nz);
    plb_ifstream geometryFile(fNameIn.c_str());
    if (!geometryFile.is_open()) {
            pcout << "Error: could not open geometry file " << fNameIn << endl;
            return -1;
    }
    geometryFile >> geometry;

    pcout << "nu = " << nu << endl;
    pcout << "deltaP = " << deltaP << endl;
    pcout << "omega = " << omega << endl;
    pcout << "nx = " << lattice.getNx() << endl;
    pcout << "ny = " << lattice.getNy() << endl;
    pcout << "nz = " << lattice.getNz() << endl;

    porousMediaSetup(lattice, createLocalBoundaryCondition3D<T,DESCRIPTOR>(), geometry, deltaP);

    // The value-tracer is used to stop the simulation once it has converged.
    // 1st parameter:velocity
    // 2nd parameter:size
    // 3rd parameters:threshold
    // 1st and second parameters ae used for the length of the time average (size/velocity)
    util::ValueTracer<T> converge(1.0,1000.0,1.0e-4);

    pcout << "Simulation begins" << endl;
    plint iT=0;

    const plint maxT = 30000;
    for (;iT<maxT; ++iT)
    {
        if (iT % 20 == 0) {
            pcout << "Iteration " << iT << endl;
        }
        if (iT % 500 == 0 && iT>0) {
            writeGifs(lattice,iT);
        }

        lattice.collideAndStream();
        converge.takeValue(getStoredAverageEnergy(lattice),true);

        if (converge.hasConverged())
        {
            break;
        }
    }

    pcout << "End of simulation at iteration " << iT << endl;

    pcout << "Permeability:" << endl << endl;
    computePermeability(lattice, nu, deltaP, lattice.getBoundingBox());
    pcout << endl;

    pcout << "Writing VTK file ..." << endl << endl;
    writeVTK(lattice,iT);
    pcout << "Finished!" << endl << endl;
}

Line 1-12: General definitions and includes are discussed in the other tutorials of Palabos. Here we use the D3Q19 lattice, but other lattices (e.g. D3Q15 or D3Q27) are also possible.
Line 16-29: This functional defines the initial conditions. It is used in function porousMediaSetup: the fluid has initially zero velocity, with a linear pressure gradient in the x-direction.
Line 31-59: Here, the boundary conditions are set, i.e. all voxels read in as 1 are defined to have bounce back boundary conditions, and all voxels read in with value two are set to carry no dynamics. The inlet and outlet are defined to have a fixed pressure difference.
Line 61-109: Various outputs are created. During the simulation .gif files showing the velocity distribution within a slice are made (line 61-80). A .vti file with the velocity distribution is written, which can be read by Paraview (line 83-88). The permeability is computed by applying Darcy’s law to the simulated velocity data (line 90-109).
Line 111-201: Main part of the code. The code requires 6 input values from the user: the input.dat file, the output directory, the size of the file in each orthogonal drection, and finally the pressure difference between the in- and outlet. From these input values, the whole run is instantiated and the simulation starts (line 174). The ValueTracer (line 168) is used to stop the simulation when steady state is reached. A maximum of number of iterations is defined (line 173). Every 500 steps a .gif file is written to the output directory (line 179). Once converged the permeability of the medium is written to the screen and a .vti file is created.

Running a simulation¶

First the code needs to be compiled. Check if the Makefile is in order. Open a terminal and type make, once in the permeability directory. To run a simulation type ./permeability path/to/input.dat path/to/outputdir/ nx ny nz deltaP in a terminal in the permeability directory. Let’s test it on our example by typing ./permeability twoSpheres.dat tmp/ 48 64 64 0.00005. The progress of the simulation is written to the screen. One can find the .gif and .vti files in the tmp subdirectory.

Post-processing¶

Use our favorite image program to visualize the .gif files. Here we see the velocity distribution at the inlet (left) and halfway (right) after 4500 iterations: