-
Notifications
You must be signed in to change notification settings - Fork 30
Simulating Flow in Digital Rock Images
In this example, permeability and relative permeability are performed based on a 3D Bennetheimer sandstone sample that is available from the Digital Rocks Portal
https://www.digitalrocksportal.org/projects/11/origin_data/17/
The image can be downloaded using the linux command line based on the following
wget https://www.digitalrocksportal.org/projects/11/images/65559/download/
The image is 1000x1000x1000 with a porosity of 0.22. The physical length of the system is 3.0035 mm in each direction. In the original image, the porespace is labeled as 0 and the solid is labeled 255 (assuming unsigned 8-bit integer type).
Domain decomposition is performed to distribute the digital image to the the MPI processes assigned to perform the simulation in parallel. LBPM uses a regular process grid with equally-sized sub-domains assigned to each MPI process. In this case, the original 1000 x 1000 x 1000 image will be distributed using 4 x 4 x 4 process grid to assign 250 x 250 x 250 sub-domains to each MPI process.
Domain {
Filename = "bentheimer.raw"
nproc = 4, 4, 4 // Number of processors (Npx,Npy,Npz)
n = 250, 250, 250 // Size of local domain (Nx,Ny,Nz)
N = 1000, 1000, 1000 // size of the input image
L = 3.0035e-3, 3.0035e-3, 3.0035e-3 // Length of domain (x,y,z)
BC = 0 // Boundary condition type
Sw = 0.20
ReadType = "8bit"
ReadValues = 0, 255
WriteValues = 2, 0
}
Units are not assumed for the length of the domain. In this case, the sample length L = 3.0035e-3, 3.0035e-3, 3.0035e-3 is provided in meters.
All that is done by lbpm_serial_decomp is to
- read the binary input file (
"8bit"and "16bit" are supported forReadType); - extract the sub-domains associated with each MPI process;
- re-label the images on the mapping between
ReadValuesandWriteValues; and - write 8-bit binary files specifying the local sub-domains associated with each MPI process
For the 64 MPI process this case the files will ID.00000 -- ID.00063). Irrespective of the process grid to be used for simulation, the domain decomposition tool should be launched with one MPI process,
mpirun -np 1 $LBPM_DIR/lbpm_serial_decomp input.db
Output resulting from successful application of the domain decomposition tool is as follows:
Input media: bentheimer.raw
Relabeling 2 values
oldvalue=0, newvalue =2
oldvalue=-1, newvalue =0
Dimensions of segmented image: 1000 x 1000 x 1000
Reading 8-bit input data
Read segmented data from bentheimer.raw
Distributing subdomains across 64 processors
Process grid: 4 x 4 x 4
Subdomain size: 250 x 250 x 250
Size of transition region: 0
Label=0, Count=222000731
Label=-1, Count=802191773
Note that the convention within LBPM is to support signed integer types. For this reason the input value of 255 is mapped to -1. For two-fluid simulations, it is assumed that the non-wetting fluid is labeled 1 and the wetting fluid is labeled 2. Any fluid labels in the input image should be mapped accordingly.
The single fluid simulater lbpm_permeability_simulator uses a multi-relaxation time (MRT) collision operator to recover the Navier-Stokes equations. The parameters for the model are specified by adding a section to th input database, labeled as MRT
MRT {
timestepMax = 100000
tau = 0.7
F = 0, 0, 1.0e-5
Restart = false
din = 1.0
dout = 1.0
flux = 0.0
}
The value of timestepMax controls how many timesteps will be executed. The permeability is obtained based on the steady-state velocity field, which can be generated by applying an external body force F = 0, 0 1.0e-5 with fully periodic boundary conditions. To launch the simulation,
mpirun -np 64 $LBPM_DIR/lbpm_permeability_simulator input.db
********************************************************
Running Single Phase Permeability Calculation
********************************************************
Read input media...
Read input media...
Initialize from segmented data: solid=0, NWP=1, WP=2
Media porosity = 0.238141
Reading in signed distance function...
Domain set.
Create ScaLBL_Communicator
Set up memory efficient layout
Allocating distributions
Setting up device map and neighbor list
Initializing distributions
Beginning AA timesteps...
********************************************************
Additionally, a time log of the average velocity and measured permeability will be written to the file Permeability.csv. The entry in the last column gives the permeability, with the units determined based on the units provided for length L in the Domain section of the input database. Given a sufficient number of timesteps, the value of the permeability should stabilize as a steady-state flow field is achieved at the pore-scale. To view a time history of the measured permeability, issue the following command from the linux command line: awk '{print $1, $10}' Permeability.csv. This time history can be used to verify whether or not steady-state has been achieved based on a particular simulation.
LBPM includes two pre-processing tools to establish initial conditions based for two-fluid flow simulations based on a morphological analysis of the porespace. As input, the pre-processors require as input 8-bit binary input files ID.xxxxx (generated by the domain decomposition) and the input file database input.db. The tool lbpm_morphopen_pp carries out a morphological opening operation based on a distance map generated from the binary input image. The tool lbpm_morphdrain_pp carries out a morphological drainage operation, which represents a morphological opening operation where disconnected portions of the non-wetting fluid are removed. For the morphological drainage process, fluid is injected from the z boundary. For both morphological drainage and morphological opening, a target saturation value (e.g. Sw = 0.20) is provided in the Domain section of the input file database. The morphological pre-processing tools should be launched in parallel, and will rely on the same domain decomposition to be used for simulation.
mpirun -np 64 $LBPM_DIR/lbpm_morphopen_pp input.db
The input files ID.00000 - ID.00063 will be updated with the final fluid distributions instantiated by the morphological approach. Example output from successful application is as follows
Target saturation 0.200000
Initialized solid phase -- Converting to Signed Distance function
Media Porosity: 0.221999
Maximum pore size: 33.619191
Performing morphological opening with target saturation 0.200000
0.973235 31.938231
0.966348 30.341320
0.963179 28.824254
0.957378 27.383041
0.950695 26.013889
0.945010 24.713195
0.938086 23.477535
0.930908 22.303658
0.923435 21.188475
0.915979 20.129051
0.905179 19.122599
0.893375 18.166469
0.882140 17.258146
0.868941 16.395238
0.853263 15.575476
0.834218 14.796703
0.813206 14.056867
0.790634 13.354024
0.767067 12.686323
0.743739 12.052007
0.714445 11.449406
0.684729 10.876936
0.655375 10.333089
0.627946 9.816435
0.592800 9.325613
0.563609 8.859332
0.532081 8.416366
0.504145 7.995547
0.477285 7.595770
0.447399 7.215982
0.424752 6.855183
0.390426 6.512423
0.375153 6.186802
0.352279 5.877462
0.329382 5.583589
0.297794 5.304410
0.297562 5.039189
0.279914 4.787230
0.249926 4.547868
0.234465 4.320475
0.228681 4.104451
0.206594 3.899228
0.195826 3.704267
Final saturation=0.195826
Final critical radius=3.704267
Writing ID file
In addition to establishing an initial condition for two-fluid flow simulation, the morphological analysis tools provide information on the resolution of the porespace. The final critical radius is reported in voxels, which corresponds to the sphere radius yielding the target saturation value. Also reported by the tool are the maximum pore size and a list of radii and saturation values obtained during the execution of the code. This information can be used to identify situations where the pore structure is insufficiently resolved to support accurate numerical simulations.
The two-fluid color lattice Boltzmann simulator lbpm_color_simulator can be used to simulate two-fluid flow based on a variety of different boundary conditions. The most straightforward way to measure the relative permeability is by simulating steady-state flow with a fixed volume fraction for each fluid, using periodic boundary conditions and an external driving force to limit boundary effects. The parameters for the color lattice Boltzmann model are specified in the Color section of the input database
Color {
tauA = 0.7;
tauB = 0.7;
rhoA = 1.0;
rhoB = 1.0;
alpha = 1e-3;
beta = 0.95;
F = 2.2681e-5, 0, 0.0e-5
Restart = false
pBC = 0
din = 1.0
dout = 1.0
timestepMax = 2500000
interval = 1000
tol = 1e-5;
das = 0.1
dbs = 0.9
flux = 0.0
capillary_number = 1.0e-4
ComponentLabels = 0
ComponentAffinity = -1.0
target_saturation = 0.05, 0.10, 0.15, 0.25, 0.35, 0.45, 0.55, 0.60, 0.65
}
There are a variety of options that determine the particular simulation protocol to be used within LBPM. In this case, we consider the specific options used to obtain steady-state relative permeability measurements. relying on an internal morphological algorithm to vary the volume fraction of the fluids. In this case, simulation parameters will be automatically modified to match a target capillary number (based on the total flow rate),
capillary_number = 1.0e-4
The fluid volume fractions will updated internally based on a morphological algorithm designed to match a sequence of target saturation values
target_saturation = 0.05, 0.10, 0.15, 0.25, 0.35, 0.45, 0.55, 0.60, 0.65
The morphological approach is applied to alter the fluid geometries after steady-state is achieved. The criterion for steady-state are determined based on the Analysis section of the input database
Analysis {
tolerance = 0.01
morph_interval = 10000
morph_delta = -0.1
blobid_interval = 1000 // Frequency to perform blob identification
analysis_interval = 1000 // Frequency to perform analysis
restart_interval = 20000 // Frequency to write restart data
visualization_interval = 1000000 // Frequency to write visualization data
restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
The error tolerance for applied to determine steady state is provided by tolerance, which is measured by comparing the measured capillary number to a previous value based on the lag specified as morph_interval. This corresponds to the number of timesteps to simulate before assessing steady-state. For this case, steady-state portion of the simulation will continue until the variation in the capillary number is less than 0.01 over a period of 10000 timesteps. Once this condition is met, the morphological approach will kick so that the subsequent value in the target_saturation list can be considered.