Size Factors and Transport Conductance#

Conductance in pore network models#

Consider an arbitrary pore network model such as the one shown below:

https://user-images.githubusercontent.com/14086031/188179071-01413efc-ed04-4920-b9b1-d99288ced47b.png

For simulating transport in this network, we need to know lots of parameters such as pore sizes, physical properties, etc, but ultimately it comes down to knowing the conductance of the “conduits” in the network. For example, consider simulating Fickian diffusion in a pore network model under steady state conditions. Writing a mass balance around an arbitrary pore $i$ , we have:

\sum_{k}^{N b^{i}} m_{i k} = \sum_{k}^{N b^{i}} g_{i k} (c_{i} - c_{k}) = 0

where $N b^{i}$ is the number of neighboring pores to pore $i$ , $m_{i k}$ is the mass flow rate from pore $i$ to pore $k$ , and $g_{i k}$ is the overall conductance of the conduit containing pores $i$ and $k$ , and the connecting throat $i k$ . In simple terms, it means that any mass that enters/exits pore $i$ must balance such that the net rate is 0. Writing the mass balance for each pore, we get a linear system of $N_{p}$ equations (where $N_{p}$ is the total number of pores in the network), which can be solved for the pore concentrations using a linear solver.

Conductance of a single conduit#

Consider a conduit in a pore network model, as shown by the following figure (we assumed spherical pores and cylindrical throats, but this is not necessary):

https://user-images.githubusercontent.com/14086031/187992821-94dee392-25a1-4fba-bce4-6ee8be8c22f2.png

where L is length, A is the cross-sectional area, and the subscripts $i$ , $k$ , and $i k$ denote pore $i$ , pore $k$ , and throat $i k$ . Now, consider the diffusional mass transfer between pore $i$ and pore $k$ with a mass flow rate of $m_{i k}$ , which can be described by the following equation:

(1)#

m_{i k} = g_{i k} \cdot (c_{i} - c_{k})

where $g_{i k}$ is the overall conductance of the conduit $i k$ .

Note

A conduit is the assembly of a throat and the two pores it connects. So, conduit $i k$ is not the same as throat $i k$ ! The conduit $i k$ begins at the center of pore $i$ and ends at the center of pore $k$ .

The overall conductance is a function of the pore sizes, throat sizes, and physical properties of the fluid. Since elements of a conduit are connected in series, by electrical circuit analogy, the resistance (inverse of the conductance) of each element can be related to the overall resistance by the following equation:

R_{i k} = R_{i}^{p} + R_{k}^{p} + R_{i k}^{t}

where the superscripts $p$ and $t$ denote pore and throat, respectively, and $R_{i k}$ is the resistance of the entire conduit. The equation above can be rewritten in form of conductance as follows:

(2)#

1 / g_{i k} = 1 / g_{i}^{p} + 1 / g_{k}^{p} + 1 / g_{i k}^{t}

For simplicity, let’s assume that the shape of the pores and throats are cylindrical as shown in the figure below:

https://user-images.githubusercontent.com/14086031/188184909-079977ba-d9a8-4760-99f7-4019b4e8f2a9.png

Based on this assumption, the diffusive conductance of each element can be expressed as (more on how to derive these equations comes later):

(3)#

g_{i}^{p} = D \frac{π d_{i}^{2} / 4}{ℓ_{i}}, g_{k}^{p} = D \frac{π d_{k}^{2} / 4}{ℓ_{k}}, g_{i k}^{t} = D \frac{π d_{i k}^{2} / 4}{ℓ_{i k}}

Size factors vs. shape factors#

In reality, however, pores/throats are not cylindrical, and the diffusive conductance needs to be corrected. This correction factor is called the shape factor in the network modeling literature, as defined by the following equation:

S F = \frac{actual conductance}{cylindrical conductance}

Warning

We used to use shape factors in OpenPNM in the past, but for the sake of clarity, we removed them and introduced the concept of “size factor” instead. More on this later.

In OpenPNM, we used to follow this terminology, but because the shape factor depends on both geometrical and physical properties, we coined the term size factor, $S$ to decouple the geometric part of the conductance from its physical properties dependence, such that the diffusional mass flow rate becomes:

m_{i k} = S \cdot D \cdot (c_{i} - c_{k})

where $D$ is the diffusivity. In other words, the size factor is the geometric part of the conductance, and so the conductance can be defined in terms of the size factor as:

g = S \cdot D

for diffusive transport, and

g = S \cdot \frac{1}{μ}

for fluid flow, where $μ$ is the viscosity.

Built-in size factor models in OpenPNM#

OpenPNM has a number of built-in functions to calculate the size factor, which can be found under openpnm.models.geometry. Let’s inspect the available diffusive and hydraulic size factor models:

import openpnm as op
op.visualization.set_mpl_style()
op.Workspace().settings.loglevel = 50
diffusive_size_factors = op.models.geometry.diffusive_size_factors
print("\n".join(diffusive_size_factors._funcs.__all__))

spheres_and_cylinders
circles_and_rectangles
cones_and_cylinders
intersecting_cones
hybrid_cones_and_cylinders
trapezoids_and_rectangles
hybrid_trapezoids_and_rectangles
intersecting_trapezoids
pyramids_and_cuboids
intersecting_pyramids
hybrid_pyramids_and_cuboids
cubes_and_cuboids
squares_and_rectangles
intersecting_trapezoids
ncylinders_in_series

import openpnm as op
hydraulic_size_factors = op.models.geometry.hydraulic_size_factors
print("\n".join(hydraulic_size_factors._funcs.__all__))

spheres_and_cylinders
circles_and_rectangles
cones_and_cylinders
intersecting_cones
hybrid_cones_and_cylinders
trapezoids_and_rectangles
intersecting_trapezoids
hybrid_trapezoids_and_rectangles
pyramids_and_cuboids
intersecting_pyramids
hybrid_pyramids_and_cuboids
cubes_and_cuboids
squares_and_rectangles
ncylinders_in_series

Note

The diffusive size factor models can be used for any physics that involves molecular diffusion such as heat conduction, electron transport, etc, and therefore, it is not limited to the diffusive mass transport. The hydraulic size factor models, however, are only applicable to fluid flow algorithms such as StokesFlow.

Deriving the size factor for arbitrary-shaped pores and throats#

In this section, we will derive the equation for the size factor for arbitrary-shaped pores and throats, then we will apply it to the case of cylindrical pores and throats to retrieve the equations we used in the previous section.

Consider the following conduit, which consists of pore $i$ , throat $i k$ , and pore $k$ :

Note that the choice of spheres and cylinders in the above figure is arbitrary and will not be taken into account in the derivation.

Diffusive size factor#

From the Fick’s law of diffusion, the diffusive mass flow rate between pore $i$ and pore $k$ is given by (assuming that the transport begins at the center of pore $i$ and ends at the center of pore $k$ ):

m_{i k} = - D \cdot A (x) \cdot \frac{d c}{d x}

where $A (x)$ is the cross-sectional area of the conduit at position $x$ , and $d c / d x$ is the concentration gradient. Rearranging the equation above, we get:

m_{i k} \cdot \frac{d x}{A (x)} = - D \cdot d c

Integrating the above equation over the entire conduit, we get:

m_{i k} \int_{i \to k} \frac{d x}{A (x)} = - D \int_{i \to k} d c

Note that $m_{i k}$ and $D$ are not a function of $x$ , so they can be taken out of the integral. Since the conduit elements are connected in series, the above equation can be rewritten as:

m_{i k} [(\int \frac{d x}{A (x)})_{pore i} + (\int \frac{d x}{A (x)})_{throat i k} + (\int \frac{d x}{A (x)})_{pore k}] = D \cdot (c_{i} - c_{k})

and after rearranging, we get:

m_{i k} = [(\int \frac{d x}{A (x)})_{pore i} + (\int \frac{d x}{A (x)})_{throat i k} + (\int \frac{d x}{A (x)})_{pore k}]^{- 1} \cdot D \cdot (c_{i} - c_{k})

By comparing the above equation with equations (1), and (2), we can get:

\frac{1}{g_{i k}} = \frac{1}{g_{i}^{p}} + \frac{1}{g_{i k}^{t}} + \frac{1}{g_{k}^{p}} = [(\int \frac{d x}{A (x)})_{pore i} + (\int \frac{d x}{A (x)})_{throat i k} + (\int \frac{d x}{A (x)})_{pore k}] \cdot \frac{1}{D}

Combining the above equation with the definition of the diffusive size factor, $g = S \cdot D$ , we get:

\frac{1}{D} \cdot \frac{1}{S_{i k}} = \frac{1}{D} \cdot [\frac{1}{S_{i}^{p}} + \frac{1}{S_{i k}^{t}} + \frac{1}{S_{k}^{p}}] = [(\int \frac{d x}{A (x)})_{pore i} + (\int \frac{d x}{A (x)})_{throat i k} + (\int \frac{d x}{A (x)})_{pore k}] \cdot \frac{1}{D}

Canceling the $1 / D$ terms, we can finally get the size factors for each conduit element:

S_{i}^{p} = [\int_{ℓ_{i}} \frac{d x}{A_{i} (x)}]^{- 1}, S_{i k}^{t} = [\int_{ℓ_{i k}} \frac{d x}{A_{i k} (x)}]^{- 1}, S_{k}^{p} = [\int_{ℓ_{k}} \frac{d x}{A_{k} (x)}]^{- 1}

where $A_{i} (x)$ , $A_{k} (x)$ , and $A_{i k} (x)$ are the cross-sectional areas of pores $i$ , $k$ , and throat $i k$ , respectively. Note that the above equations are valid for any arbitrary shape of the pores and throats.

As a test case, let’s consider cylindrical pores and throats. In this case, the cross-sectional area of the each element is given by:

(4)#

A_{i} = π d_{i}^{2} / 4, A_{k} = π d_{k}^{2} / 4, A_{i k} = π d_{i k}^{2} / 4

where $d_{i}$ , $d_{k}$ , and $d_{i k}$ are the diameters of pores $i$ , $k$ , and throat $i k$ , respectively. Substituting the above equations into the above equations, we get:

S_{i}^{p} = \frac{π d_{i}^{2}}{4 \cdot ℓ_{i}}, S_{i k}^{t} = \frac{π d_{i k}^{2}}{4 \cdot ℓ_{i k}}, S_{k}^{p} = \frac{π d_{k}^{2}}{4 \cdot ℓ_{k}}

which matches what we declared earlier in equation (3).

Hydraulic size factor#

For brevity, we will not go into the nitty-gritty details of the derivation of the hydraulic size factors. Instead, we will give you some pointers on where to start. The derivation of the hydraulic size factors is slightly more involved than that of the diffusive size factors. The main reason is that fluid flow equations need to satisfy two sets of boundary conditions, pressure drop along the domain ( $p_{i}$ and $p_{k}$ at the two ends of the conduit), and no slip at the lateral boundaries (walls), whereas for diffusive transport, only the concentration boundary conditions need to be satisfied ( $c_{i}$ and $c_{k}$ at the two ends of the conduit).

Assuming that the cross-section of the conduit is gradually changing, the pressure drop along the conduit can be approximated by the following equation [1]:

- \frac{1}{Q} \frac{d p}{d x} = μ [16 π^{2} \cdot \frac{I_{p}^{*}}{A (x)^{2}} - \frac{2 ρ Q}{μ} \cdot \frac{d A (x) / d x}{A (x)^{3}}]

where $Q$ is the volumetric flow rate, $I_{p}^{*}$ is the specific polar moment of inertia [2] (don’t worry if you don’t know what this is, there are tables in the literature that you can use to look up the values for different shapes), $A (x)$ is the cross-sectional area of the conduit, and $μ$ is the viscosity of the fluid. Rearranging and integrating the above equation, we get:

(5)#

- \frac{1}{Q} \int d p = μ [16 π^{2} \int \frac{I_{p}^{*}}{A (x)} d x - \frac{2 ρ Q}{μ} \cdot \int \frac{d A (x) / d x}{A (x)^{3}} d x]

From calculus, we know that:

\frac{d A (x) / d x}{A (x)^{3}} = - \frac{1}{2} \frac{d}{d x} (\frac{1}{A (x)^{2}})

Substituting the above equation into equation (5), we get:

\frac{p_{i} - p_{j}}{Q} = μ [16 π^{2} \int \frac{I_{p}^{*}}{A (x)} d x + \frac{ρ Q}{μ} \cdot \int d (\frac{1}{A (x)^{2}})]

which can further be simplified to:

\frac{p_{i} - p_{j}}{Q} = 16 π^{2} μ \int \frac{I_{p}^{*}}{A (x)^{2}} d x + ρ Q \cdot (\frac{1}{A_{k}^{2}} - \frac{1}{A_{i}^{2}})

where the first term in the right-hand-side is the frictional loss, and the second term is the inertial loss. The frictional loss is a function of the conduit geometry, and the inertial loss is a function of the flow rate and the conduit geometry. The frictional loss is usually much larger than the inertial loss, so we can safely ignore the inertial loss to obtain:

(6)#

\frac{p_{i} - p_{j}}{Q} = 16 π^{2} μ \int \frac{I_{p}^{*}}{A (x)^{2}} d x

Similar to the diffusive transport, the hydraulic conductance is defined as the ratio of the flow rate to the pressure drop as follows:

(7)#

Q = g_{i k} \cdot (p_{i} - p_{k})

where $g_{i k}$ is the hydraulic conductance (we use the same notation as the diffusive conductance to keep it simple in this notebook). Expanding the integral over the conduit on each element (pores $i$ and $k$ , and throat $i k$ ), and by comparing equations (6) and (7), we get:

\frac{1}{g_{i k}} = \frac{1}{g_{i}^{p}} + \frac{1}{g_{i k}^{t}} + \frac{1}{g_{k}^{p}} = 16 π^{2} μ [(\int \frac{I_{p}^{*}}{A_{i} (x)^{2}} d x)_{pore i} + (\int \frac{I_{p}^{*}}{A_{i k} (x)^{2}} d x)_{throat i k} + (\int \frac{I_{p}^{*}}{A_{k} (x)^{2}} d x)_{pore k}]

Combining the above equation with the definition of the hydraulic size factor, $g = S / μ$ , we get:

\frac{μ}{S_{i k}} = μ \cdot [\frac{1}{S_{i}^{p}} + \frac{1}{S_{i k}^{t}} + \frac{1}{S_{k}^{p}}] = 16 π^{2} μ [(\int \frac{I_{p}^{*}}{A_{i} (x)^{2}} d x)_{pore i} + (\int \frac{I_{p}^{*}}{A_{i k} (x)^{2}} d x)_{throat i k} + (\int \frac{I_{p}^{*}}{A_{k} (x)^{2}} d x)_{pore k}]

from which, we can obtain the hydraulic size factor for each conduit element as follows:

(8)#

S_{i}^{p} = [16 π^{2} \int_{ℓ_{i}} \frac{I_{p}^{*}}{A_{i} (x)^{2}} d x]^{- 1}, S_{i k}^{t} = [16 π^{2} \int_{ℓ_{i k}} \frac{I_{p}^{*}}{A_{i k} (x)^{2}} d x]^{- 1}, S_{k}^{p} = [16 π^{2} \int_{ℓ_{k}} \frac{I_{p}^{*}}{A_{k} (x)^{2}} d x]^{- 1}

As a test case, consider cylidnrical pores and throats. The specific polar moment of inertia for a cylinder is:

I_{p}^{*} = \frac{1}{4 π}

Substituting the above equation and the cross-sectional areas from equation (4) into equation (8), we get:

S_{i}^{p} = \frac{π d_{i}^{4}}{128}, S_{i k}^{t} = \frac{π d_{i k}^{4}}{128}, S_{k}^{p} = \frac{π d_{k}^{4}}{128}

To verify the above results, we can inspect the Hagen-Poiseuille equation for a cylindrical conduit [3]:

Δ p = \frac{8 π μ ℓ Q}{A^{2}} \Rightarrow \frac{Q}{Δ p} = \frac{A^{2}}{8 π μ ℓ} = \frac{π d^{4}}{128 μ ℓ}

Since the hydraulic size factor is defined as $g / μ$ , and the hydraulic conductance is defined as $Q / Δ p$ , we can derive the following relationship for the hydraulic size factor based on the Hagen-Poiseuille equation:

S = \frac{π d^{4}}{128 ℓ}

which matches our derivation.

Conduit vs. element size factor#

So far, we have derived individual equations for the size factor of each conduit element as well the entire conduit (based on resistors in series model). Note that eventually, only the size factor of the entire conduit is used by the algorithms, although as we discussed, it is directly calculated from the size factors of the individual elements.

Almost all of the built-in size factor models in OpenPNM return an Nt by 3 array, where Nt is the number of conduits (same as the number of throats!). The first column is the size factor of the pore on the left side of the conduit, the second column is the size factor of the throat, and the third column is the size factor of the pore on the right side of the conduit.

However, we do have one model (at the time of writing this notebook) that uses a machine learning method [4] to predict the diffusive size factor of the entire conduit directly from a voxel image of the porous medium. For this reason, it only returns an Nt by 1 vector. Our built-in conductance models (which take in the size factor and physical properties as input and return the conductance) are designed to work with both types of size factor models.