# Stokes Flow and Permeability Tensor#

[1]:

import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import numpy as np
import matplotlib.pyplot as plt


## Generate a Layered Network#

In order to see a network with a noticable difference in permeability values in each principle direction, we’ll make a network consisting of two distinct layers then stitch them together. The process of stitching is explained in more detail here.

[2]:

pn = op.network.Cubic(shape=[8, 8, 8], spacing=30e-6)
pn2 = op.network.Cubic(shape=[24, 24, 8], spacing=10e-6)


### Apply Labels to Each Network#

In order to make it easier to find pore indices after stitching, it’s helpful to apply unique and descriptive labels to the networks first. Labels are discussed detail here.

[3]:

pn['pore.coarse'] = True
pn['throat.coarse'] = True
pn2['pore.fine'] = True
pn2['throat.fine'] = True


Now we can stitch the two network together, but first we must translate the fine network along the z-axis so the networks are not overlapping. Stiching is described in more detail here.

[4]:

pn2['pore.coords'] += [0, 0, 8*30e-6]
op.topotools.stitch(network=pn, donor=pn2,
P_network=pn.pores('top'),
P_donor=pn2.pores('bottom'))


Let’s visualize the network to see if it stitched as expected. Quick plotting is described in more detail here.

[5]:

fig = op.topotools.plot_connections(pn)


## Create other objects needed for simulation#

### Generate one Geometry for each layer#

[6]:

np.random.seed(0)
Ps = pn.pores('fine')
Ts = pn.throats('fine')
geo2 = op.geometry.SpheresAndCylinders(network=pn, pores=Ps, throats=Ts)
Ps = pn.pores('coarse')
Ts = pn.throats(['coarse', 'stitched'])
geo1 = op.geometry.SpheresAndCylinders(network=pn, pores=Ps, throats=Ts)


Let’s visualize the pore sizes using histograms:

[7]:

fig = plt.hist(geo1['pore.diameter'], bins=25,
density=True, edgecolor='k', alpha=0.5)
fig = plt.hist(geo2['pore.diameter'], bins=25,
density=True, edgecolor='k', alpha=0.5)


### Create a Phase and 2 Physics objects#

[8]:

air = op.phase.Air(network=pn)
phys1 = op.physics.Standard(network=pn, phase=air, geometry=geo1)
phys2 = op.physics.Standard(network=pn, phase=air, geometry=geo2)


Let’s visualize the hydraulic conductances using histograms. Note that the values are logged, to better visualize the distributions:

[9]:

fig = plt.hist(np.log10(phys1['throat.hydraulic_conductance']),
bins=25, density=True, edgecolor='k', alpha=0.5)
fig = plt.hist(np.log10(phys2['throat.hydraulic_conductance']),
bins=25, density=True, edgecolor='k', alpha=0.5)


## Run a StokesFlow simulation in perpendicular directions#

[10]:

sf_x = op.algorithms.StokesFlow(network=pn, phase=air)
Pin = 1.0
Pout = 0.0
sf_x.set_value_BC(pores=pn.pores('left'), values=1.0)
sf_x.set_value_BC(pores=pn.pores('right'), values=0.0)
sf_x.run()

[10]:

{'pore.pressure': SteadyStateSolution([1., 1., 1., ..., 0., 0., 0.])}


Let’s make sure the pressure field is applied in the correct direction:

[11]:

fig = op.topotools.plot_coordinates(pn, color=sf_x['pore.pressure'],
size_by=pn['pore.diameter'],
markersize=50)

[12]:

sf_z = op.algorithms.StokesFlow(network=pn, phase=air)
Pin = 1.0
Pout = 0.0
sf_z.set_value_BC(pores=pn.pores(['top', 'coarse'], mode='and'), values=Pin)
sf_z.set_value_BC(pores=pn.pores(['top', 'fine'], mode='and'), values=Pout)
sf_z.run()

[12]:

{'pore.pressure': SteadyStateSolution([1.        , 1.        , 1.        , ..., 0.13355113,
0.02637116, 0.        ])}


Again, let’s visualize the pressure field:

[13]:

fig = op.topotools.plot_coordinates(pn, color=sf_z['pore.pressure'],
size_by=pn['pore.diameter'],
markersize=50)


The presure in the coarse network is essentially constant since it’s much more permeable than the fine layer.

### Determine the Permeability Coefficient for each direction#

[14]:

L = 8 * 30e-6
A = (8 * 30e-6)**2 + (8 * 24)*(10e-6)**2
mu = air['pore.viscosity'][0]
Q = sf_x.rate(pores=pn.pores('left'), mode='group')
Kx = Q*L*mu/(A*(Pin - Pout))
print('The permeability coefficient is:', Kx)

The permeability coefficient is: [4.23480546e-14]

[15]:

A = (8 * 30e-6)**2
L = (8 * 30e-6) + (8 * 10e-6)
mu = air['pore.viscosity'][0]
Q = sf_z.rate(pores=pn.pores(['top', 'coarse'], mode='and'), mode='group')
Kz = Q*L*mu/(A*(Pin - Pout))
print('The permeability coefficient is:', Kz)

The permeability coefficient is: [2.4235904e-14]

[16]:

Kx/Kz

[16]:

array([1.74732722])


The permeability coefficient for the direction through the ‘fine’ layer is about 2x lower than the direction parallel to it.