Simulating Single Phase Transport#

The point of an OpenPNM simulation is ultimately to compute some transport process. In this notebook we will cover the following subjects:

Defining conductance
Settings boundary conditions

Start by defining a network. We’ll use the Demo class which happens to include all the geometrical pore-scale models already.

import numpy as np
import openpnm as op
op.visualization.set_mpl_style()

pn = op.network.Demo(shape=[5, 5, 1], spacing=5e-5)
print(pn)

══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Demo at 0x7f9730bd8400>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
pore.coords                                                       25 / 25
throat.conns                                                      40 / 40
pore.coordination_number                                          25 / 25
pore.max_size                                                     25 / 25
throat.spacing                                                    40 / 40
pore.seed                                                         25 / 25
pore.diameter                                                     25 / 25
throat.max_size                                                   40 / 40
throat.diameter                                                   40 / 40
throat.cross_sectional_area                                       40 / 40
throat.hydraulic_size_factors                                     40 / 40
throat.diffusive_size_factors                                     40 / 40
throat.lens_volume                                                40 / 40
throat.length                                                     40 / 40
throat.total_volume                                               40 / 40
throat.volume                                                     40 / 40
pore.volume                                                       25 / 25
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
pore.xmin                                                               5
pore.xmax                                                               5
pore.ymin                                                               5
pore.ymax                                                               5
pore.surface                                                           16
throat.surface                                                         16
pore.left                                                               5
pore.right                                                              5
pore.front                                                              5
pore.back                                                               5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Define Phase Viscosity#

To fully illustrate the process of performing transport calculations, we’ll use an empty Phase object and add all the needed properties manually:

water = op.phase.Phase(network=pn)

Let’s assume that we are interested in pressure driven flow. This requires knowing the viscosity of the phase, so let’s add a pore-scale model for computing the viscosity of water:

water.add_model(propname='pore.viscosity',
                model=op.models.phase.viscosity.water_correlation)
print(water)

══════════════════════════════════════════════════════════════════════════════
phase_01 : <openpnm.phase.Phase at 0x7f972e657ba0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.temperature                                                  25 / 25
  3  pore.pressure                                                     25 / 25
  4  pore.viscosity                                                    25 / 25
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                               25
  3  throat.all                                                             40
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

And we can check the individual values to verify they make sense:

print(water['pore.viscosity'])

[0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
00089319]

Basic Conductance Calculation#

Determining the conductance of the conduits between each pair of connected pores is the most important part of performing a simulation. The details of conductance models are covered elsewhere. For this demonstration will assume the very simplest case where all pressure loss occurs in the throats.

Recall the Hagan-Poiseuille equation for fluid flow through a cylindrical tube:

\[ Q = \frac{\pi R^4}{8 \mu L} \Delta P\]

where \(R\) and \(L\) are the radius and length of the throat, and \(\mu\) is the viscosity of the fluid. Together this prefactor can be referred to as the hydraulic conductance, \(g_h\), giving:

\[ Q = g_h \Delta P \]

So the aim is the compute values of \(g_h\) for each throat. We start by doing this manually:

R = pn['throat.diameter']/2
L = pn['throat.length']
mu = water['throat.viscosity']  # See ProTip below
water['throat.hydraulic_conductance'] = np.pi*R**4/(8*mu*L)
print(water['throat.hydraulic_conductance'])

[3.90656713e-15 4.61477128e-15 4.61693012e-15 3.95923078e-15
88089160e-14 3.91085689e-15 3.77158450e-15 3.75567926e-14
17466676e-16 9.87652813e-16 1.94529152e-14 2.08205534e-14
48310587e-14 4.99528793e-15 4.29615905e-15 7.59624050e-15
05489737e-15 3.06522239e-15 1.25043212e-15 6.36912282e-16
92014514e-15 4.66669842e-15 3.86786205e-15 4.50061035e-15
47152255e-15 1.41282395e-14 1.04517963e-15 3.67973864e-15
00191458e-14 4.00479482e-14 1.13306108e-14 1.01532282e-15
84605856e-15 7.10227941e-15 3.90730022e-14 1.57420050e-14
89930283e-15 5.29210640e-15 1.03354292e-15 7.98783750e-16]

Phase can do automatic interpolation to get throat values

The Phase class has a special ability to interpolate pore values to throats, and vice-versa. In PNM simulations all the balances are solved for each pore, so the thermodynamic properties like temperature, pressure, etc. are all defined on pores. Consequently, the physical properties are also defined in pores, like viscosity. However, as shown above we often want viscosity values in the throats. OpenPNM provides a shortcut for this, such that if you request a throat property that does not exist, it will attempt to fetch the pores values and do a linear interpolation of values to produce an array of throat values. There is also a function for this, water.interpolate_data('throat.viscosity'), but the water['throat.viscosity'] shortcut is very convenient. The automatic interpolation can be disabled in the phase.settings.

water.interpolate_data('throat.viscosity')

array([0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319])

water['throat.viscosity']

array([0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319])

Create Algorithm Object#

OpenPNM contains a variety of Algorithm classes in the openpnm.algorithms module. Let’s initialize a StokesFlow algorithm, since this simulates pressure driven flow through the network.

sf = op.algorithms.StokesFlow(network=pn, phase=water)
print(sf)

══════════════════════════════════════════════════════════════════════════════
stokes_01 : <openpnm.algorithms.StokesFlow at 0x7f972e679620>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.bc.rate                                                       0 / 25
  3  pore.bc.value                                                      0 / 25
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                               25
  3  throat.all                                                             40
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Assign Boundary Conditions#

As can be seen in the print-out above there are predefined 'pore.bc' arrays, but they contain no valid values, meaning they are all nans. Once we set some boundary conditions, this will change. Let’s apply pressure BCs on one side of the network, and rate BCs on the other:

sf.set_value_BC(pores=pn.pores('left'), values=100_000)
sf.set_rate_BC(pores=pn.pores('right'), rates=1e-10)
print(sf)

══════════════════════════════════════════════════════════════════════════════
stokes_01 : <openpnm.algorithms.StokesFlow at 0x7f972e679620>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.bc.rate                                                       5 / 25
  3  pore.bc.value                                                      5 / 25
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                               25
  3  throat.all                                                             40
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

All boundary conditions are preceeded with ‘pore.bc’

All boundary conditions are stored as ‘pore.bc.’, which means that OpenPNM’s dictionary lookup tricks can be used to see all types and values of bc’s using: sf['pore.bc'] which will return a dict. This can be used as shown below:

print(sf['pore.bc'].keys())

dict_keys(['rate', 'value'])

Now we can see there are 5 valid values of each type. The sf algorithm will look for 'throat.hydraulic_conductance' on water by default, so we can just run:

soln = sf.run()

The run method solves the mass balance around each pore and computes the pressure within each pore that is required to sustain the flow defined by the boundary conditions. The soln object that is returned is a dictionary with the key corresponding to the quantity that was solved for.

print(soln)

None

The reason for the dict format is to provide a consistent API for single components and multiphysics, where multiple different quantities might be solved for. However, these 'pore.pressure' values are also written to the dictionary of the algorithm object as well:

print(sf)

══════════════════════════════════════════════════════════════════════════════
stokes_01 : <openpnm.algorithms.StokesFlow at 0x7f972e679620>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.bc.rate                                                       5 / 25
  3  pore.bc.value                                                      5 / 25
  4  pore.pressure                                                     25 / 25
  5  pore.initial_guess                                                25 / 25
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                               25
  3  throat.all                                                             40
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Finally we can look at how much pressure was required in the “right” pores to meet the required flow rate:

sf['pore.pressure'][pn.pores('right')]

array([152894.58128806, 167635.75510035, 170052.05136202, 204084.42360037,
       229894.24284339])

So we can see that 150 kPa was required to accomplish the requested flow.

Rigorous Conductance Calculation#

The above demonstration used a very simplistic conductance calculation. It was also stated that computing conductance is the most important part of doing a PNM simulation. To finish this notebook, we’ll look more closely at this process.

Manual Method#

Let’s print the network object again:

print(pn)

══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Demo at 0x7f9730bd8400>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
pore.coords                                                       25 / 25
throat.conns                                                      40 / 40
pore.coordination_number                                          25 / 25
pore.max_size                                                     25 / 25
throat.spacing                                                    40 / 40
pore.seed                                                         25 / 25
pore.diameter                                                     25 / 25
throat.max_size                                                   40 / 40
throat.diameter                                                   40 / 40
throat.cross_sectional_area                                       40 / 40
throat.hydraulic_size_factors                                     40 / 40
throat.diffusive_size_factors                                     40 / 40
throat.lens_volume                                                40 / 40
throat.length                                                     40 / 40
throat.total_volume                                               40 / 40
throat.volume                                                     40 / 40
pore.volume                                                       25 / 25
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
pore.xmin                                                               5
pore.xmax                                                               5
pore.ymin                                                               5
pore.ymax                                                               5
pore.surface                                                           16
throat.surface                                                         16
pore.left                                                               5
pore.right                                                              5
pore.front                                                              5
pore.back                                                               5
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Note the 'throat.hydraulic_size_factors' array. This is computed by a pore-scale model on the Demo network. This computation is more rigorous in the following ways:

The conductance of each half pore and that throat is considered.
The throat length is computed carefully by accounting for the ‘lens’ between the intersection of the spherical pore bodies and the cylindrical throat.
The net cross-sectional area of the pores are computed by integrating between the pore center and the pore-throat intersection point

The conductance of each element in the conduit is returned as an Nt-by-3 array, where columns 1 and 3 contain the hydraulic conductance of the half pore on either end of the throat, and the column 1 contains the throat conductance.

pn['throat.hydraulic_size_factors']

array([[1.10847142e-16, 3.48931040e-18, 9.55917974e-17],
       [9.55917974e-17, 4.12187193e-18, 1.77142522e-16],
       [1.77179362e-16, 4.12380018e-18, 9.56236884e-17],
       [9.56236884e-17, 3.53634909e-18, 1.16037553e-16],
       [6.05768484e-16, 5.25275914e-17, 5.28879136e-16],
       [1.67682687e-16, 3.49314198e-18, 8.42711436e-17],
       [8.42711436e-17, 3.36874514e-18, 1.54647046e-16],
       [4.26436245e-16, 3.35453872e-17, 4.29098937e-16],
       [6.41050434e-17, 8.19472930e-19, 3.26315557e-17],
       [3.26315557e-17, 8.82162552e-19, 8.03762465e-17],
       [3.04467390e-16, 1.73751677e-17, 2.75269824e-16],
       [2.75269824e-16, 1.85967298e-17, 3.48386558e-16],
       [2.22008324e-16, 1.32469674e-17, 2.85778952e-16],
       [1.76937093e-16, 4.46174596e-18, 1.02635642e-16],
       [1.02635642e-16, 3.83729037e-18, 1.14536153e-16],
       [1.38835697e-16, 6.78489325e-18, 2.10964639e-16],
       [1.55508680e-16, 2.72860667e-18, 6.95411101e-17],
       [6.95411101e-17, 2.73782889e-18, 1.56562581e-16],
       [1.04193408e-16, 1.11687465e-18, 3.65137220e-17],
       [3.04885685e-17, 5.68884285e-19, 2.75139571e-17],
       [1.43410974e-16, 7.96739289e-18, 2.68791181e-16],
       [9.55917974e-17, 4.16825278e-18, 1.81415456e-16],
       [1.63744738e-16, 3.45473937e-18, 8.42711436e-17],
       [9.56236884e-17, 4.01990442e-18, 1.67309772e-16],
       [1.62021063e-16, 8.45987820e-18, 2.33702165e-16],
       [3.28484229e-16, 1.26192156e-17, 1.97927213e-16],
       [9.22257146e-17, 9.33544982e-19, 3.26315557e-17],
       [8.42711436e-17, 3.28670924e-18, 1.45600520e-16],
       [3.23035819e-16, 1.78809198e-17, 2.75269824e-16],
       [4.33948438e-16, 3.57704648e-17, 4.63746441e-16],
       [1.97927213e-16, 1.01203990e-17, 2.06256203e-16],
       [3.26315557e-17, 9.06877149e-19, 8.62334635e-17],
       [1.64581311e-16, 4.32845564e-18, 1.02635642e-16],
       [1.78726559e-16, 6.34369167e-18, 1.38835697e-16],
       [4.58746612e-16, 3.48996519e-17, 4.26601608e-16],
       [2.22008324e-16, 1.40606163e-17, 3.20645023e-16],
       [1.38807853e-16, 2.58963104e-18, 6.95411101e-17],
       [1.02635642e-16, 4.72686153e-18, 1.99798897e-16],
       [6.01486062e-17, 9.23151182e-19, 3.65137220e-17],
       [7.65146779e-17, 7.13466414e-19, 2.75139571e-17]])

This data is called the size factor because it is purely the geometrical information required for the computation of the hydraulic conductance. So the Hagan-Poisseiulle equation for each element is written as:

\[ Q = \frac{F_h}{\mu} \Delta P = g_h \Delta P\]

Note that both the \(\frac{\pi R^4}{8}\) term and \(L\) have been rolled into the \(F_h\) value.

The total conductance of the pore-throat-pore conduit can be found as the sum of three resistors in series. Since we have conductance values, we add the inverses, and invert again. The full expression for the hydraulic conductance between pores i and j, through throat k, is:

\[ Q = \bigg( \frac{\mu}{F_{h, i}} + \frac{\mu}{F_{h, k}} + \frac{\mu}{F_{h, j}} \bigg) ^ {-1} \Delta P \]

This can be computed by hand:

F_h = water['throat.hydraulic_size_factors']
water['throat.hydraulic_conductance'] = (mu * (1/F_h).sum(axis=1))**(-1)

water['throat.hydraulic_conductance']

array([3.65790024e-15, 4.32747785e-15, 4.32945427e-15, 3.70902778e-15,
       4.95846353e-14, 3.68155811e-15, 3.55220556e-15, 3.24649611e-14,
       8.83967677e-16, 9.51487307e-16, 1.73657606e-14, 1.85742277e-14,
       1.34093622e-14, 4.67422361e-15, 4.01176361e-15, 7.02684648e-15,
       2.89075003e-15, 2.90031880e-15, 1.20082939e-15, 6.12807414e-16,
       8.21983236e-15, 4.37538119e-15, 3.64173286e-15, 4.22170149e-15,
       8.70213158e-15, 1.28185275e-14, 1.00620813e-15, 3.46630092e-15,
       1.78692828e-14, 3.45370778e-14, 1.02986917e-14, 9.77862894e-16,
       4.53550005e-15, 6.56897143e-15, 3.37451601e-14, 1.42180469e-14,
       2.74582474e-15, 4.94722138e-15, 9.93189486e-16, 7.71581161e-16])

sf = op.algorithms.StokesFlow(network=pn, phase=water)
sf.set_value_BC(pores=pn.pores('left'), values=100_000)
sf.set_rate_BC(pores=pn.pores('right'), rates=1e-10)
soln = sf.run()
sf['pore.pressure'][pn.pores('right')]

array([157627.11848238, 172977.15697151, 175460.26969393, 210314.98342933,
       236560.40550604])

As can be seen the numbers are about the same as with the simple case, but should be somewhat more correct. In fact, these above pressures are a bit higher, which is because the total conductance of the conduit is lower due to the inclusion of the pore body lengths into the total length, compared to above where only the throat length was included.

Pore-Scale Model Method#

Instead of computing the hydraulic conductance manually as done above, there is a pore-scale model available:

water.add_model(propname='throat.hydraulic_conductance',
                model=op.models.physics.hydraulic_conductance.generic_hydraulic)

sf = op.algorithms.StokesFlow(network=pn, phase=water)
sf.set_value_BC(pores=pn.pores('left'), values=100_000)
sf.set_rate_BC(pores=pn.pores('right'), rates=1e-10)
soln = sf.run()
sf['pore.pressure'][pn.pores('right')]

array([157627.11848238, 172977.15697151, 175460.26969393, 210314.98342933,
       236560.40550604])

Which gives exactly the same result, without have to manually deal with the conductances-in-series calculation.