Stokes Flow and Permeability Tensor#

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.

pn =[8, 8, 8], spacing=30e-6)
pn2 =[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.

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.

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

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

fig = op.topotools.plot_connections(pn)

Create other objects needed for simulation#

Generate one Geometry for each layer#

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:

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#

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:

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#

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)
{'pore.pressure': SteadyStateSolution([1., 1., 1., ..., 0., 0., 0.])}

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

fig = op.topotools.plot_coordinates(pn, color=sf_x['pore.pressure'],
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)
{'pore.pressure': SteadyStateSolution([1.        , 1.        , 1.        , ..., 0.13355113,
                      0.02637116, 0.        ])}

Again, let’s visualize the pressure field:

fig = op.topotools.plot_coordinates(pn, color=sf_z['pore.pressure'],

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#

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]
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]

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