PoreFoam Single-Phase Flow Simulation

Single-phase electrical and flow simulation on 3D images of porous media.

This document presents the details of the porefoam codes and scripts and the instructions for using them to run single-phase flow simulations on micro-CT images of porous media. The porefoam package consists of pre-processing, processing and and post-processing applications, which includes standard C++ codes as well as applications written using OpenFOAM. Several bash scripts are developed to simplify the workflow. This document presents a description of the keywords, file-name conventions and the scripts used to do direct single-phase flow simulations.

Installation

Prerequisites

Except standard Linux compilers and libraries (g++, cmake and an mpi library), all other prerequisites are provided in a folder named thirdparty. The thirdparty includes zlib, libtiff and a minified openfoam.

The minified openfoam library, after being compiled, can reside side-by-side with other openfoam installations withought any conflict. This means you can install and use other openfoam versions alongside the porefoam codes, but if you do so, you better delete the /works/thirdparty/foam-miniext-3.4/applications contents so that you don’t have mutiple copies of the same application.

Downloading and extracting files

Create a directory in your home (~) folder named works and download and extract the codes inside it. You can choose any other folder and the scripts will work, but this ducument assumes you extract the source codes in the ~/works floder for simplicity. The source codes, will reside in subfolders ~/works/thirdparthy and ~/works/apps (called psSrc in the scripts). When compiling the codes, the executables, libraries and the header files will be generated in ~/works/bin ~/works/lib ~/works/include folders.

To check that the directories are created correctly, type in a terminal

#  Important: replace ~/works/apps with the
#  directory in which you have extracted the codes
ls   ~/works/apps/scripts   ~/works/apps/porefoam* ~/works/apps/voxelImage

and it should not show the error message No such file or directory.

Compiling the codes

The codes requires a recent C++ compiler, which support -std=c++11 flag. The compiler is set through the variable psCXX in file ~/works/apps/Makefile.in. The default (g++) most likely will work so you don’t need to do any change.

To compile the codes, open a terminal and type the following commands:

(cd ~/works/ && make -j)

After everything compiled and working, you can run the following commands to clean the temporary files

(cd ~/works/ && make clean)

If instead of the clean argument you use distclean, the ~/works/lib, ~/works/bin, ~/works/include and ~/works/shared folders will be deleted, so be extra carefull when using this option.

Setting the Environmental variables:

Edit your ~/.bashrc file (type gedit ~/.bashrc to open it) and add the following line at its end:

source ~/works/apps/bashrc

This makes the porefoam scripts available in any new terminal you open.

Usage

Input file format

The required input files are:

An input header file.
A micro-CT image, which can be in ascii (suffix should be .dat), or binary (better to have .raw suffix, the data should be in 8bit unsigned char), gzip compressed raw, for which the file suffix should be .raw.gz, or a 3D .tif file. You can use ImageJ, or the imageFileConvert from the voxelImage library to convert between different image types.

The input file for the flow simulation codes is a .mhd header file, which is a standard format compatible with Paraview and Fiji (ImageJ with plugins), with additional optional image manipulation commands specific to porefoam package.

The main keywords of the .mhd header are as follows. (First to third keywords should not be changed):

ObjectType = Image
NDims = 3
ElementType = MET_UCHAR
keyword: DimSize — used to assign the dimensions of the image: Nx, Ny and Nz
keyword: ElementSpacing — used for assigning voxel size: , and should be equal
keyword: ElementDataFile — specifies the name of binary 8bit data file (.raw), ascii (.dat) files are supported too, but only by the porefoam codes. For compressed images, the file suffix should be .raw.gz. The suffix for tiff files should be ‘.tif’. The script can be run on the .tif file directly if the ‘XResolution’ tiff tag is set correctly.

    ObjectType =  Image
    NDims =       3
    ElementType = MET_UCHAR

    DimSize =        400    400    400
    ElementSpacing = 5.345  5.345  5.345
    Offset =         0       0       0

    ElementDataFile = Berea.raw

    ####optional keywords, uncomment to activate:
    #BinaryData = True
    #threshold  0 0   # void space range in the image
    #direction  z       # turn image in y or z direction
    #cropD       0 0 0 200 200 200   # crop the image
    #resample   2   # coarsen (>1) or refine (<1)

Notes:

The order of the first 7 keywords should not be changed. The remaining keywords are optional.

Characters # and // are used for marking the rest of a line as comments

All keyword and its data should be given in a single line

For more advanced features you can edit the scripts in the ~/works/apps/scripts/singlePhase folder, other input parameters which are not handled through the scripts can be edited manually by editing the files in ~/works/apps/scripts/singlePhase/base folder that follows, more or less, the OpenFOAM input file format conventions.

Running simulations

A sample input file is provided in docs/Image.mhd, which is a ascii header file that should be used in combination with a micro-CT image, such as those available from Imperial College Pore-scale modelling website.

To run a single-phase flow simulation on a micro-CT image, first the header file should be prepared as discussed in the previous section, then they should be copied in a directory where there is enough disk space; the 3D simulation results take disk-spaces on the order of Gigabytes, depending on the size of the image. Then in this directory run one of the following bash commands:

# 1. Set the Environmental variables if you have not done so far:
#   First, you have to add the scripts for single-phase
#   flow calculation to the system path. To do so, type:
bash # just in case
export PATH=$PATH:$HOME/works/apps/scripts/singlePhase
export nProcX=2  # optional
export nProcY=2  # optional
export nProcZ=3  # optional
# 2. Running the single-phase flow solver,

# to get help for how to run the scripts, type:
AllRunImagePar -h

# to run on image XXX.mhd in serial mode, type:
AllRunImage XXX.mhd
# or to run on all .raw files in the current directory in serial mode:
AllRunImage
#  or to run on XXX.raw file in parallel mode on a single machine:
AllRunImagePar XXX.mhd
#  or to run on all .mhd files in parallel on a single machine
AllRunImagePar
#  or, for parallel mode on a distributed cluster,
#  edit ~/works/apps/scripts/singlePhase/machines.txt first and type:
AllRunImageParDistributed XXX.mhd

The same approach can be used to run the simulation on ascii files, provided that you have .dat as the file suffix (e.g. XXX.dat); of course you should replace “XXX.raw” to the name of ascii file, XXX.dat, or set the keyword BinaryData to False in the .mhd file.

Simulation results and visualization

After running the simulation scrips the results will be saved in a subdirectory inside a directory with the same name as the micro-CT image (without the suffix), which includes the log files of individual applications, openfoam simulation results and the pressure and velocity fields converted into image files, in either tif, .dat (ascii) or .raw format: To change the default behavior, you have to open the script files AllRunImage/AllRunImagePar (located in folder ~/works/apps/scripts/singlePhase/) and change : ${outPutFormat:=tif} to : ${outPutFormat:=raw} or : ${outPutFormat:=raw.gz} or : ${outPutFormat:=dat}.

In binary format, the velocity and pressure fields are written in small-endian 32bit floating-point numbers. When oldAscii/oldBinary format is chosen the velocity files will have the same size as the original image, but otherwise the (default) Ufx, Ufy and Ufz files will have one additional layers in x, y and z directions respectively (to match the number of faces between all pairs of neighbouring voxels). They can be opened in image-processing software like Avizo and ImageJ.

The openfoam simulation parameters and results are saved inside a folder called openfoam case, which includes subdirectories named constant and system and in many other time directories which have floating-point numbers as their names such as 0 or 0.1, etc. To visualize the openfoam files, you can use Paraview and open the system/controlDict file. The permeability and effective (connected) porosity are reported in a file with prefix summary_ in the case folder. Velocity distributions and formation factor are also computed and written into the same file.

Upscaling and flow statistics

The direct single-phase flow simulation results are post-processed by calc_distributions code and its results are written in a file named summary_CASENAME.txt, where case name is the name of the folder in which the simulations are performed, which is in the format of IMAGENAME-DelP-FLOWDIR, e.g. Berea-1-X. The following are short descriptions of various data written in this file:

effPorosity: effective (largest connected path) porosity, computed based on the connected path bounding box, not the total image. Bounding box used are also provided in case a different definition is desired.
V_pore: pore volume, of connected pores
L_x, L_y, L_z: lengths of connected pores bounding box, in x, y and z directions respectively.
K_x, K_y, K_z: permeability of connected pores, in x, y and z directions respectively.
FF_x, FF_x, FF_x: formation factor of connected pores, in x, y and z directions respectively.
DarcyVelocity: Darcy velocity, average of inlet and (=) outlet fluxes divided by
U_D: Darcy velocity, average of inlet and (=) outlet fluxes divided by
DelP: average pressure drop
Umax: maximum of magnitude of velocity vector field
Re: Reynold number =

Tabular data with header x=mag(U)/U_D PDF dV/Vdx PDF(log10(x)) dV/Vd(log10(x)):

x=mag(U)/U_D: x-axis in velocity distribution, magnitude of velocity at bin centres, normalized by Darcy velocity ()
PDF: probably density function of
dV/Vdx: same as PDF, probably density function of , volume weighted but gives same values as PDF above in typical (no sub-resolution porosity) cases.
PDF(log10(x)): probably density function of ,
dV/Vd(log10(x)): same as PDF(log10(x)) in typical (no sub-resolution porosity) cases.

Tabular data with header x=U_x/U_D PDF dV/Vdx:

x=U_x/U_D: x axis in velocity distribution, x-component of velocity at bin centres, normalized by Darcy velocity
PDF: probably density function of x=U_x/U_D
dV/Vdx: probably density function of x=U_x/U_D, volume weighted.

Tabular data with header x=U_y/U_D PDF dV/Vdx: Same as U_x distribution above, but for y-component of velocity.

Tabular data with header x=U_z/U_D PDF dV/Vdx: Same as U_x distribution above, but for z-component of velocity.

| Tabular data with header x=U_m/U_D dummy PDF dV/Vdx distConstDelCbrtU: | Velocity distributions obtained by binning the velocity in uniform intervals in the domain of cubic-root of velocity. These can then be plotted more accurately as most velocities are centred near zero. These data has not been used in any publication, so double-checking for consistency with results above are needed.

x=U_m/U_D: Velocity magnitude at bin centres, normalized by Darcy velocity
dummy: To be ignored, this is just a reminder for lack of testing of these data
PDF: probably density function of x=U_x/U_D
dV/Vdx: probably density function of x=U_x/U_D, volume weighted
distConstDelCbrtU: To be ignored

| Tabular data with header x=U_x/U_D dummy PDF dV/Vdx: | Same as table above for x=U_m/U_D, but for x-component of velocity field, U_x instead of U_m

| Tabular data with header x=U_y/U_D dummy PDF dV/Vdx: | Same as table above for x=U_m/U_D, but for x-component of velocity field, U_x instead of U_m

| Tabular data with header x=U_z/U_D dummy PDF dV/Vdx: | Same as table above for x=U_m/U_D, but for x-component of velocity field, U_x instead of U_m

References

The single-phase flow solver is derived from the two-phase direct simulation code discussed in the paper:

Raeini, A.Q., Blunt, M. J. and Bijeljic B. (2012). Modelling Two-Phase Flow in Porous Media at the Pore Scale Using the Volume-of-Fluid Method, Journal of Computational Physics, 231, 5653-5668 http://dx.doi.org/10.1016/j.jcp.2012.04.011

Essentially the main difference with the openfoam-provided codes is the treatment of the no-slip boundary condition for improving the accuracy of permeability prediction. The details of the pre/post-processing tools are not documented, but you can have a look at the source code or send an email (see below) if you want to know more.

Applications of the single-phase direct simulation code can be found in following papers:

Bijeljic B, Raeini, A. Mostaghimi P., and Blunt, M.J. “Predictions of non-Fickian solute transport in different classes of porous media using direct simulation on pore-scale images”, Physical Review E. 87, 013011 (2013), http://dx.doi.org/10.1103/PhysRevE.87.013011
Bijeljic, B., P. Mostaghimi, and M. J. Blunt “Insights into non-Fickian solute transport in carbonates”, Water Resources Research, 49, 2714-2728, (2013) http://dx.doi.org/10.1002/wrcr.20238
Pereira Nunes, J.P., Bijeljic, B. and Blunt, M.J. “Time-of-Flight Distributions and Breakthrough Curves in Heterogeneous Porous Media Using a Pore-Scale Streamline Tracing Algorithm”, Transport in Porous Media (2015) 109: 317. https://doi.org/10.1007/s11242-015-0520-y
B.P. Muljadi, , M.J. Blunt, A.Q. Raeini, B. Bijeljic, “The impact of porous media heterogeneity on non-Darcy flow behaviour from pore-scale simulation”, Advances in Water Resources, (2016) http://dx.doi.org/10.1016/j.advwatres.2015.05.019
\…

For more information and recent publications, please refer to our website:

http://www.imperial.ac.uk/earth-science/research/research-groups/perm/research/pore-scale-modelling/publications/

If you have any questions, please contact us by email:

Ali Qaseminejad Raeini: a.qaseminejad-raeini09@imperial.ac.uk

Please also report any problems, feedback or suggestions using the email adress above.

.. porefoam1f:

Overview the applications used in the scripts

OpenFOAM utilities:

renumberMesh:

: renumbers mesh to improve performance (optional)

redistributePar:

: used to redistribute mesh between processors (deactivated doesn’t work always)

Solvers and utilities not available in OpenFOAM:

voxelToFoam(Par):

: converts .raw/.tif/.am files file to OpenFOAM format, used for single-phase flow simulation. It also does simple image processing tasks like cropping and thresholding. Moreover, it lables the disconnected parts of the image and only keeps the largest connected pathway (to avoid linear solver crash). A parallel version of this code is developed but not included in all versions of poreFoam.

iInterFoam101SP:

: single-phase flow solver, trimmed version of the two-phase flow solver, the main differences with the OpenFOAM single-phase flow solvers is the implementation of the boundary conditions and the pressure-velocity coupling.

iPotentialFoam:

: used for initialization of pressure and velocity for improving the convergence of simulations (optional, effective when the Reynold number is low and the image is a typical micro-CT image of a porous rock).

calc_distributions:

: post-processing, calculates effective porosity and permeability, formation factor and velocity distributions.

FOAM2Voxel:

: post-processing, converts Openfoam simulation results into .dat, .raw and .tif files.

Additional useful utilities:

Ufraw2Uc:

: converts face-centred Uf*.raw files into cell-centred velocities, useful for importing velocities in Avizo \…

voxelImageConvert:

: converts between image formats, applies cropping, resampling etc, supports basic data types other than unsigned char (default); see .mhd header file specifications: https://itk.org/Wiki/ITK/MetaIO/Documentation for more details.

Using zlib

The voxelImage library can read from .gz compressed raw binary files directly, and write .gz compressed files directly if zlib is installed in the thirdparty folder. i.e. ~/works/thirdparty/zlib-1.2.11. The original zlib package can be downloaded from http://zlib.net/. The zlib directory is set in the ~/works/apps/Makefile.in file. Note that the voxelImage library needs the contrib/iostream3/zfstream.cc file to be compiled with zlib library, so it is easier if you use the one provided zlib in the thirdparty folder.

Using libtiff

The voxelImage library can read from .tif 3D images, and write .tif files if libtiff is installed in the thirdparty folder. i.e. ~/works/thirdparty/libtiff. The libtiff directory is set in the ~/works/apps/Makefile.in file. The default tif compression is set to LZW, in the /works/apps/voxelImage/voxelTiff.h file. The libtiff package can be downloaded from https://gitlab.com/libtiff/libtiff and is provided in the thirdparty folder just for convinience.