Auto calculate and plot the representative elementary volume (REV) of a porous medium using porosity
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
#from joblib import Parallel, delayed
sample_array=[]
n_slice=30
image_shape=(35,35)
for i in range(n_slice):
i= np.random.randint(2, size=image_shape)
sample_array.append(i)
sample_array=np.asarray(sample_array)
print ('generating a 3D volume of size',sample_array.shape)
generating a synthetic 3D sample volume of size (30, 35, 35)
Generate a random coordinate for subsampling
import random
def gen_pos(shape,size):
x=random.randint(0, min(shape)-size)
y=random.randint(0, min(shape)-size)
z=random.randint(0, min(shape)-size)
return x,y,z #slice, row, colomn
gen_pos(sample_array.shape,10)
(12, 3, 4)
Generate a 3D subsampling window with given coordinate and window size
def window3d(x,y,z,size,sample_array):
subsample=[]
for i in range(size):
for j in range(size):
subsample.append(sample_array[i+x][j+y][z:z+size])
return np.asarray(np.split(np.asarray(subsample),size))
Calculate volume of voids and total volume in the subsampled window and return porosity
def porosity(a):
void=np.count_nonzero(a==0)
volume=a.shape[0]*a.shape[1]*a.shape[2]
phi=void/volume
return phi
Iteratively generate random subsampling windows and calculate porosity for 10 times, with fixed window size
def rev_sampling(sample_array,size,n_iter=10):
revlist=[]
for i in range(n_iter): #number of iterations for each window size
x,y,z=gen_pos(sample_array.shape,size)
sample=window3d(x,y,z,size,sample_array)
poro=porosity(sample)
revlist.append(poro)
return revlist
Increase window size step-wisely and plot the porosity measurements
def plot(size,data,tp,fsize=20):
plt.figure(figsize=(15,9))
for ind,i in enumerate(data):
x=np.full(shape=10,fill_value=((ind+1)*size)**3,dtype=np.int)
plt.plot(x,i,'.')
plt.xlabel('Subsample volume',fontsize=fsize)
plt.ylabel('Porosity',fontsize=fsize)
plt.xticks(fontsize=fsize)
plt.yticks(fontsize=fsize)
plt.axhline(tp,0,1,label='sample total porosity')
plt.legend(fontsize=fsize)
return
def rev(sample_array,step=10): #step must be smaller than the sample lenghts
size=round(min(sample_array.shape)/step)-1
rev_all_list=[]
for i in tqdm(range(step)):
rev_all_list.append(rev_sampling(sample_array,(i+1)*size))
total_porosity=porosity(sample_array)
plot(size,rev_all_list,tp=total_porosity)
return rev_all_list
rev0=rev(sample_array,step=15)
100%|█████████████████████████████████████████| 15/15 [00:00<00:00, 468.75it/s]�[A
Because synthetic data are randomly generated with 0s and 1s, therefore the porosity should be close to 0.5. From the plot we can see porosity converging to 0.5 as the subsample volume gets larger. The representative elementary for the synthetic data is 750 because from that point on, the porosity measurements converged and stabilised around 0.5.
Now calculate REV with real porous rock data
import glob
import os
from PIL import Image
folder='N:\\PRESALT_FLOW\\FROM CTLAB\\PRESALT FLOW\FF4402H_EXPT_1_FEB_2018\\H2O_Saturated\converted\\NLM'
list_of_imgs=sorted(glob.glob(os.path.join(folder,'*.tif')))
def import_img(file):
img=Image.open(file)
img=np.asarray(img)[196:825,170:850]
return img
img0=import_img(list_of_imgs[0])
plt.imshow(img0,cmap='gray')
<matplotlib.image.AxesImage at 0x4279608>
Apply a simple thresholding to segment the porosity and rock
porespace=np.zeros(img0.shape)
porespace[img0>38]=1
plt.imshow(porespace,cmap='gray')
<matplotlib.image.AxesImage at 0x42fd2c8>
Load the image data into a 3D np.array
sample_3d=[]
#n_slice=len(list_of_imgs)
for i in tqdm(list_of_imgs):
img= import_img(i)
porespace=np.zeros(img.shape)
porespace[img>38]=1
sample_3d.append(porespace)
sample_3darray=np.asarray(sample_3d)
print ('generating a 3D volume of size',sample_3darray.shape)
print ('total porosity:',porosity(sample_3darray))
100%|███████████████████████████████████████| 887/887 [00:04<00:00, 183.83it/s]
generating a 3D volume of size (887, 629, 680)
Now run the rev function and show the rev plot
rev_f4402h=rev(sample_3darray,step=10)
0%| | 0/10 [00:00<?, ?it/s]�[A
10%|████▎ | 1/10 [00:00<00:01, 6.58it/s]�[A
20%|████████▌ | 2/10 [00:00<00:02, 3.20it/s]�[A
30%|████████████▉ | 3/10 [00:02<00:05, 1.33it/s]�[A
40%|█████████████████▏ | 4/10 [00:06<00:09, 1.56s/it]�[A
50%|█████████████████████▌ | 5/10 [00:12<00:15, 3.04s/it]�[A
60%|█████████████████████████▊ | 6/10 [00:23<00:21, 5.49s/it]�[A
70%|██████████████████████████████ | 7/10 [00:40<00:26, 8.78s/it]�[A
80%|██████████████████████████████████▍ | 8/10 [01:04<00:26, 13.37s/it]�[A
90%|██████████████████████████████████████▋ | 9/10 [01:36<00:19, 19.00s/it]�[A
100%|██████████████████████████████████████████| 10/10 [02:21<00:00, 14.11s/it]�[A
The subsample porosity measurements converged and stabilised around the sample total porosity (0.1) from around 0.5e8 voxels, therefore it is the REV of this sample.