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)
../../../_images/examples_simulations_steady_state_stokes_flow_and_permeability_tensor_9_0.svg

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)
../../../_images/examples_simulations_steady_state_stokes_flow_and_permeability_tensor_13_0.svg

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)
../../../_images/examples_simulations_steady_state_stokes_flow_and_permeability_tensor_17_0.svg

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)
../../../_images/examples_simulations_steady_state_stokes_flow_and_permeability_tensor_21_0.svg
[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)
../../../_images/examples_simulations_steady_state_stokes_flow_and_permeability_tensor_24_0.svg

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.