Are there better ways to generate mulitscale porous medias ?

[wangkeupc]

Network modelling of multiscale porous medias (such as shale or carbonate) is hot topic . I know some methods such as subdivide , extend pores , stitch and dual network offer good ideas and ways to do this . And these methods are limmited for complex porous media . So I want to know are there better ways in OpenPNM to generate mulitscale porous medias ?

[jgostick]

Hi @wangkeupc We have been working on a class that generates Vuggy carbonates. It's not quite ready for merging to Openpnm yet as we finalize V2. It's basically a faster way to do merging of pores. I think a better way would be to generate random points in space that are clustered in a certain way, then send them to the Delaunay class. The would naturally create multiscale networks, based on how closely the base points are clustered. I am not sure how to do this, but it's probably not too hard.

[jgostick]

Here is a very rough first attempt to illustrate what I mean...it does not work with openpnm yet:

import scipy as sp
import scipy.ndimage as spim
import matplotlib.pyplot as plt
from skimage.morphology import disk, ball, rectangle
import scipy.spatial as sptl

def correlated_random_field(shape, strel=None):
    im = sp.rand(*shape)
    if strel is None:
        if im.ndim == 3:
            strel = ball(3)
        else:
            strel = disk(3)
    im = spim.convolve(im, strel)
    # Convolution is no longer randomly distributed, so fit a gaussian
    # and find it's seeds
    im = (im - sp.mean(im))/sp.std(im)
    im = 1/2*sp.special.erfc(-im/sp.sqrt(2))
    return im

if __name__ == '__main__':
    s = rectangle(75, 25)
    # s = disk(10)
    f = correlated_random_field(shape=[1200, 800], strel=s)
    plt.imshow(sp.fliplr(spim.rotate(f, angle=-90)), cmap=plt.cm.viridis)
    n = 20000
    m = 0
    pts = sp.zeros([n, f.ndim])
    while m < n:
        x, y = sp.rand(2)
        i, j = sp.array([x*f.shape[0]-1, y*f.shape[1]-1], dtype=int)
        if (f[i, j] > 0.5*sp.rand()):
            pts[m] = sp.around((x, y), decimals=3)
            m += 1
    pts = pts*f.shape
    tri = sptl.Delaunay(pts)
    ax = plt.gca()
    sptl.delaunay_plot_2d(tri, ax=ax)
    plt.axis('off')

[jgostick]

If you run this you'll get a nice 2D picture of hierarchically spaced points. It doesn't work with openpnm yet because there are a lot of duplicate points that need to be removed, and the edges need to be dealt with. Anyway, let me know if this is interesting to you and I can develop it further.