# Pore-network models and effective medium theory:
A convergence analysis

Jack Edwards<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of Alberta  Peter Berg<EMAIL_ADDRESS>Theoretical Physics
Institute,
University of Alberta

###### Abstract

The convergence between effective medium theory and pore-network modelling is
examined. Electrical conductance on two and three-dimensional cubic resistor
networks is used as an example of transport through composite materials or
porous media. Effective conductance values are calculated for the networks
using effective medium theory and pore-network models. The convergence between
these values is analyzed as a function of network size. Effective medium
theory results are calculated analytically and numerically. Pore-network
results are calculated numerically using Monte Carlo sampling. The reduced
standard deviations of the Monte Carlo sampled pore-network results are
examined as a function of network size. Finally, a “quasi-two-dimensional”
network is investigated to demonstrate the limitations of effective medium
theory when applied to thin porous media. Power law fits are made to these
data to develop simple models governing convergence. These can be used as a
guide for future research that uses both effective medium theory and pore-
network models.

##  I Introduction

EMT is a useful tool for calculating the macroscopic properties of a composite
material based upon the properties of its constituent parts. As defined by
Choy [7], at its most basic, effective medium theory (EMT) “functions by being
able to define averages, which one hopes will be representative of the system
and be connected with experimental measurements.” While the ability to
directly calculate macro- or meso-scale behaviour based upon microscopic
properties has improved with increasing computational power, it remains
unfeasible to do so for all applications. Replacing a composite material with
a uniform effective medium that replicates its larger-scale properties can
greatly simplify calculations involving that material.

PNM techniques are another tool for analyzing the macroscopic behaviour of
materials based upon the physics at the microscopic level. They are used to
model porous media in a variety of fields, most notably Earth sciences, though
it has seen more widespread use in recent years [25]. The underlying logic of
pore-network model (PNM) is to model porous media as a system of
interconnected pores, on which transport equations can be calculated exactly.
The porous medium can then be treated similarly to a network of resistors, and
reduced to a linear system of equations analogous to Kirchhoff’s Circuit Laws
or a more complex, nonlinear system of equations. Ideally, the aggregate
behaviour of the medium can be calculated through this method [12].

The goal of this report is to examine the convergence behaviour between EMT
and PNM models. To help contextualize these techniques, a brief history of EMT
is recounted in Section II.1. The development of PNM is similarly detailed in
Section II.2. In Section II.3, the relationship between the two is explored
and their continued development and use is presented. The central questions
about this relationship that we hope to answer are detailed in Section III.

The methodology used in this analysis is detailed in Section IV. To narrow the
scope, electrical conductance on resistor networks is examined using a
particular EMT derived by Kirkpatrick [15] and PNM techniques. The open-source
PNM Python package OpenPNM is used for the pore-network calculations [14].
Effective conductances are calculated for both two-dimensional and three-
dimensional networks with a variety of resistor conductance distributions.
Through examining networks with different properties, the effect of differing
network properties on the convergence between EMT and PNM techniques can be
determined. Additionally, a “quasi-two-dimensional” network is examined, which
is very large in two dimensions but small in the third. This is to simulate
how a thin porous medium can affect the convergence between these techniques.
The results are presented in Section V and a discussion into their
implications is presented in Section VI.

The ultimate aim of this report is to guide the development of research into
composite and porous materials, especially in regards to transport problems.
Through the examination of the convergence behaviours between the macroscopic
models of EMT and the discretized, microscopic approaches used in PNM,
underlying trends governing their relationship are revealed. This can be used
as a framework for future efforts in these fields, in determining what
problems are likely suitable for EMT treatment and what domain sizes can be
considered sufficient for PNM to adequately model macroscopic properties.

##  II Literature Review

###  II.1 Effective medium theory

Attempts to describe the macroscopic properties of composite materials can be
traced back to the nineteenth century. Early work in the field can be
attributed to the physicists Lorenz [19] and Lorentz [18]. This research
related the refractive indices of a particular medium to its density,
modelling the transport of light through the medium as it pertains to its
average properties. They associated the transport in the dilute limit, where
the obstacles to transport can be considered small in comparison to the
distance between them. Lord Rayleigh built upon this foundation, expanding it
to transport obstacles whose sizes are not negligible in regards to their
spacing [22]. In this treatment, Rayleigh looked at the case of spherical or
cylindrical objects with uniform conductance arranged in a “rectangular or
square” (cubic) order, embedded within a conducting medium. Rayleigh’s
solution for the effective transport coefficient of this two-phase medium is

$\frac{D_{\mathrm{eff}}-D_{2}}{D_{\mathrm{eff}}+mD_{2}}=\frac{\nu_{1}}{\nu_{2}}\frac{D_{1}-D_{2}}{D_{1}+mD_{2}}$
((1))

when $\nu_{1}\ll\nu_{2}$, where $D_{i}$ and $\nu_{i}$ refers to the transport
coefficient and volume fraction of the two phases, respectively (adapted from
Tjaden et al. [23]). The variable $m$ is set to a value of 1 for cylinders and
2 for spheres.

These initial approaches, though often framed from the standpoint of optics
and light refraction, had more general applications to heat or electrical
conduction. These developments were important, but limited in their
applicability. In the cases of Lorenz and Lorentz, their models were valid
only in the dilute limit. For Rayleigh, this model was only valid for objects
arranged in cubic lattices in which the flux of the system is parallel to one
of the lattice constants.

The first traditional EMT was the Maxwell Garnett formula, itself an extension
of the earlier Clausius–Mossotti formula [7, 13]. It provides a method for
calculating an effective (or bulk) dielectric constant, $D_{\mathrm{eff}}$, of
a composite medium. If considering spherical inclusions of material 1 embedded
into material 0, the equation takes the form

$\frac{D_{\mathrm{eff}}-D_{0}}{D_{\mathrm{eff}}+2D_{0}}=\nu_{1}\frac{D_{1}-D_{0}}{D_{1}+2D_{0}}.$
((2))

This is similar to Equation (1). However, the coefficient in front of the
right-hand-side fraction in Equation (2) is only the volume fraction of the
inclusions, as opposed to the ratio of the volume fractions of both materials
in Equation (1). This allows the equation to be generalized to $n$ materials,

$\frac{D_{\mathrm{eff}}-D_{0}}{D_{\mathrm{eff}}+2D_{0}}=\sum_{i=1}^{n}\nu_{i}\frac{D_{i}-D_{0}}{D_{i}+2D_{0}}.$
((3))

However, this EMT was only accurate at small values of $\nu_{i}$. Taken to the
extreme, if $\nu_{1}=1$ in Equation (2) (i.e. only material 1 is present), the
calculated $D_{\mathrm{eff}}$ is still dependent on the value of $D_{0}$
despite material 0 being completely absent from the composite. It also does
not model critical threshold behaviour.

The next major development of EMTs came from Bruggeman [5]. Bruggeman expanded
Rayleigh’s approach from strictly cubic lattices to random structures, and did
so in a way that improved upon the shortcomings present in the Maxwell Garnett
formula [23, 7]. Bruggeman examined a sample with two distinct phases, 1 and
2, embedded within some homogeneous medium. Bruggeman’s model was based upon
three underlying assumptions:

  1. 1.

The two phases are homogeneous and isotropic.

  2. 2.

The particle sizes of the phases are small in comparison to the overall sample
size.

  3. 3.

The phases have been randomly distributed within the sample.

Through an iterative embedding procedure, the entire sample volume is filled
with increasingly greater amounts of phases 1 and 2, completely replacing the
embedding medium. Eventually, as the entire volume is filled with inclusions
of these phases, the following result is found:

$\nu_{1}\frac{D_{1}-D_{\mathrm{eff}}}{D_{1}+mD_{\mathrm{eff}}}+\nu_{2}\frac{D_{2}-D_{\mathrm{eff}}}{D_{2}+mD_{\mathrm{eff}}}=0.$
((4))

Much like the Maxwell Garnett formula, the terms in Equation (4) are dependent
only on the volume fraction $\nu_{i}$ of a single phase. This allows it to be
generalized to an arbitrary number of phases, similar to Equation (3):

$\sum_{i=1}^{n}\nu_{i}\frac{D_{i}-D_{\mathrm{eff}}}{D_{i}+mD_{\mathrm{eff}}}=0.$
((5))

This is the main result from Bruggeman’s analysis, which has become a widely
used EMT across many disciplines.

Considering how important this result has been for the development of EMTs,
the barriers to accessing the original publication have hindered discussions
of the result. It was published in German and, to our knowledge, has not been
translated into other languages. However, excellent explanations of the theory
are given by both Tjaden et al. [23] and Choy [7].

###  II.2 Pore-network modelling

The techniques of EMT can also be applied to transport through porous media,
analogously to their applications in optics and electrical conductance.
However, transport through porous media has been studied for longer than EMTs
have been available, and there exists a history of distinct techniques in this
field. The earliest empirical studies can be traced to Henry Darcy and the
development of Darcy’s Law [8]. This research, based upon the flow of water
through soil and sand, yielded a simple transport equation for a fluid with a
known viscosity flowing through a medium with a known permeability under the
influence of a known pressure gradient. One limitation of Darcy’s Law is that
it only applies to laminar flows, reducing its applicability in certain
scenarios (for instance, turbulent flows through gravel). Nonetheless, Darcy’s
Law continues to be widely used.

Attempts to describe the transport through porous media in more detail
followed. Models involving systems of packed spheres were used, notably the
Kozeny–Carman equation [12]. Through this method, the permeability of the
medium was related to the porosity and geometry of the sphere packs, with some
proportionality constant. However, the complexity of the pore-geometry within
this sphere-packed model prevented accurate mathematical derivations. This
method, though a step toward more complex models, proved to be of limited use
outside of specific fields (e.g. estimating the surface area of powders) [12].

The next attempts involved modelling porous media as bundles of tubes through
which fluid flows [12]. This geometry allowed for comparatively simpler and
rigorous mathematical solutions. These gains made in mathematical rigour,
however, were offset by the model’s inability to accurately model several key
features of porous media. For example, porous media tend to be isotropic,
whereas the bundle of tubes model is anisotropic. Additionally, these
approaches were incapable of modelling multiphase flow phenomena.

In order to address some of the shortcomings present in these models, Fatt
[12] took elements from both the packed-spheres models and bundle of tubes
models to develop what he referred to as the network model, or network of
tubes. This is now most-commonly referred to as a pore-network model (PNM).
Fatt [12] took the interconnectedness of the packed-spheres models, which
better represents the isotropic properties of most porous media; and combined
it with the cylindrical tubes of the bundle of tubes models, which are much
simpler mathematically and allow for exact flow calculations to be performed.
In this treatment, porous media can essentially be modelled as networks of
resistors. Fatt [12] discusses how manipulating the distribution of cylinder
radii within the network, lattice structure of the network, and other
geometric properties can be used to model the properties of a variety of
different materials. Indeed, it is this aspect of PNMs that makes them useful
to this day.

PNMs were initially mostly used within fields of Earth sciences, modelling the
flow of fluids through porous rock or soil. Since then, their uses and
applications have spread into many other fields, from nuclear waste management
to electrochemical energy systems [25, 11, 3].

###  II.3 Pore-network modelling and effective medium theory

When Fatt [12] first introduced PNMs, it was an important achievement in the
modelling of porous media. The ability to introduce random physical
characteristics to the network, as well as the notions of “coordination
number” to describe the connectivity of individual nodes within the network,
was a pioneering achievement. However, Fatt’s treatments were limited by the
computing power available at the time (the initial calculations were performed
by hand) [17, 12]. The networks considered were necessarily small and results
were dominated by boundary effects. Additionally, randomized characteristics
had to be described through discrete distributions, as opposed to more
physically-accurate continuous distributions. Much like in Bruggeman’s time,
EMTs remained a necessity in describing macroscopic properties of porous or
composite media. Direct calculations through PNM techniques were not feasible;
only effective medium approximations could be practically used.

As computing technology advanced, PNMs were able to grow larger and more
complex. Extensions to Fatt’s models were developed, adding more detail and
increasing the size of the domains examined [21, 10]. Comparisons between the
predictions of PNMs and EMT were made, examining the applicability of each in
various situations [15, 11, 4, 6]. Today, computational power is such that
many problems can be “brute-forced” in PNM and similar models.

Nevertheless, EMTs are still extensively used. The Bruggeman model continues
to be wide-spread, recently seeing greater applications in fields such as fuel
cell and battery research [23, 24, 9, 2]. When applied to porous media,
transport through one of the phases is often negligible (e.g. water will flow
in the spaces between grains of sand, but not through the sand itself). This
allows Equation (4) to be simplified to

$D_{\mathrm{eff}}=D\frac{\nu}{\tau},$ ((6))

where $D$ and $\nu$ are the transport coefficient and volume fraction of the
transporting phase. $\tau$ is the tortuosity factor, given by

$\tau=\begin{cases}\nu^{-1/2},&\text{for spherical inclusions}\\\
\nu^{-1},&\text{for cylindrical inclusions}\end{cases}.$ ((7))

Equation (6) is the EMT that is most often applied to porous media transport
[23].

##  III Problem Statement

In this report, we will examine the applicability of EMT to a given transport
problem within a composite material. Specifically, we will calculate the
effective electrical conductance within a composite material, though the
results will be applicable to other problems which can be modelled with
similar PNM or resistor techniques. This is to answer several related
questions: When calculating effective transport properties in a composite
medium, how large does the system have to be for EMT to model it accurately?
How large of a network is needed for PNM techniques to agree with the
macroscale limit that EMT describes? What is the rate of convergence between
PNM and EMT results as the system size grows? What effect does the specific
random conductance distribution have on this behaviour?

##  IV Methods

###  IV.1 Effective medium theory

To address the problem in Section III, it is necessary to compare EMT
transport solutions to those directly calculated by PNM techniques. To narrow
the focus of this work, we will specifically look at an EMT developed by
Kirkpatrick [15], a formulation based upon Bruggeman’s work. Kirkpatrick’s
formulation considers the case of a network of random resistors, arranged in a
square (2D) or cubic (3D) lattice. An example network is shown in Figure 1.

Within such a network, the electric potential differences across the various
resistors are generated by an external electric field. The potentials can then
be treated as consisting of two contributions: the portion supplied by the
external field, which increases the voltage at a given node by a uniform
amount per row; and the portion from a fluctuating local field, which averages
to zero over a sufficiently large area. The average effects of the random
resistors are represented through an effective medium, where the total field
inside it is equal to the external field. For this to hold, the medium must be
homogeneous, represented as every resistor having an equal conductance
$g_{m}$, as seen in Figure 1. This formulation, used by Kirkpatrick, only
considers identical resistors in nearest-neighbour arrangement on a cubic
lattice 111“For simplicity we shall consider [the network] to be made up of a
set of equal conductances, $g_{m}$, connecting nearest neighbors on the cubic
mesh,” Kirkpatrick [15].. As Kirchhoff’s Laws stipulate that the total current
entering and exiting a node in the network must be zero, this allows the
network to be described as a linear system of equations,

$\sum_{i}g_{ij}\left(V_{i}-V_{j}\right)=0,$ ((8))

where $g_{ij}$ represents conductance, $V_{k}$ is electric potential, and the
subscripts $i$ refer to the nodes adjacent to node $j$ on the cubic network.

Figure 1: Example resistor network. The random network has been replaced with
an effective medium of uniform conductance.

Now, the issue is to determine the value of $g_{m}$ that replaces the random
conductances of the network such that the predicted effective properties are
accurate for the given network. To do this, it must be ensured that the local
fields in the effective medium average to zero. Kirkpatrick achieves this by
considering a single resistor in the network, oriented along the external
field, with a conductance of $g_{AB}=g_{0}$. This is seen in Figure 1. The
uniform field solution of this network has the voltages increasing by a
constant amount per row, $V_{m}$, caused by the external field. The effects of
a fictitious current $i_{0}$ are added to this system, introduced through node
$A$ and extracted from node $B$. This causes an imbalance in the current
conservation described in Equation (8). The new current equation caused by
this imbalance is

$V_{m}\left(g_{m}-g_{0}\right)=i_{0}.$ ((9))

There is an extra voltage induced between $A$ and $B$, introduced by $i_{0}$.
It can be calculated using the conductance $G^{\prime}_{AB}$, the value for
the equivalent resistor replacing the network with the resistor $g_{0}$
removed. This equivalent network is shown in Figure 2. The extra voltage is
given by

$V_{0}=\frac{i_{0}}{g_{m}+G^{\prime}_{AB}}.$ ((10))

The conductance of the total effective medium, with the missing resistor
reintroduced with $g_{m}$, is given by $G_{AB}=G^{\prime}_{AB}+g_{m}$.

Figure 2: Example resistor network. The random network has been replaced with
an equivalent resistor with conductance $G^{\prime}_{AB}$.

Examining the network as shown in Figure 1, the conductance $g_{0}$ can be
replaced with the effective medium $g_{m}$. The current distribution in this
network can be described as the sum of two components: a current $i_{0}$
introduced at $A$ and extracted at infinity, and an equal current introduced
at infinity and extracted at $B$. In both cases, the current flowing through
the $z$ equivalent bonds at the point the current enters the $AB$ resistor is
$i_{0}/z$. This gives a total current of $2i_{0}/z$ flowing through the $AB$
bond. From this, we know that $G_{AB}=(z/2)g_{m}$, and subsequently that
$G^{\prime}_{AB}=(z/2-1)g_{m}$. Knowing now the current across this resistor,
and knowing its conductance to be $g_{m}$, the voltage across this resistor
can be determined. Substituting Equation (9) into Equation (10), and
substituting for $G^{\prime}_{AB}$, we obtain the expression

$V_{0}=\frac{V_{m}\left(g_{m}-g_{0}\right)}{g_{0}+\left(\frac{z}{2}-1\right)g_{m}}.$
((11))

As the conductances in the resistor networks being examined are random, we can
distribute the conductances with a probability density function (PDF) $f(g)$,
which may be continuous or discrete. We can then define a function $F(g_{m})$,

$F(g_{m})=\int
f(g)\frac{V_{0}(g)}{V_{m}}\mathop{}\\!\mathrm{d}g=\int\frac{f(g)(g_{m}-g)}{g+\left(\frac{z}{2}-1\right)g_{m}}\mathop{}\\!\mathrm{d}g.$
((12))

As stated previously, it is required that $V_{0}$ averages to zero. This
requirement leads to the expression

$F(g_{m})=\int\frac{f(g)(g_{m}-g)}{g+\left(\frac{z}{2}-1\right)g_{m}}\mathop{}\\!\mathrm{d}g=0,$
((13))

which can be used to determine the value of $g_{m}$. In essence, the solution
to Equation (13) determines the effective conductance $g_{m}$ that can replace
the random network with a uniform effective medium. Given any idealized square
or cubic infinite lattice, with a known distribution of conductances,
transport calculations can be performed exactly on it using EMT.

For certain discrete distributions, Equation (13) has analytic solutions. For
example, take a discrete binary distribution

$f(g)=p\delta(g-\alpha)+(1-p)\delta(g-\beta),$ ((14))

where $\delta(x)$ is the Dirac delta function, $p$ is the probability that a
given bond has conductance $\alpha$, and $(1-p)$ is the probability that a
given bond has conductance $\beta$. Substituting Equation (14) into Equation
(13), we obtain the quadratic equation

$\left(\frac{z}{2}-1\right)g_{m}^{2}+\left\\{\beta\left[\frac{z}{2}\left(p-1\right)+1\right]+\alpha\left[1-\frac{z}{2}p\right]\right\\}g_{m}\\\
-\alpha\beta=0.$ ((15))

There are two roots in this equation, positive and negative. The appropriate
value of $g_{m}$ is the positive root. For any discrete distribution of $N$
positive delta functions (e.g. a ternary distribution), Equation (13) will
result in a polynomial equation of order $N$.

For a continuous PDF, numerical techniques can be used to determine the
solution. To show that Equation (12) has exactly one positive root, the
derivative with respect to $g_{m}$ is needed. This is given by

$F^{\prime}(g_{m})=\frac{\mathop{}\\!\mathrm{d}F(g_{m})}{\mathop{}\\!\mathrm{d}g_{m}}=\int\frac{f(g)\frac{z}{2}g}{\left(g+\left(\frac{z}{2}-1\right)g_{m}\right)^{2}}\mathop{}\\!\mathrm{d}g.$
((16))

The numerator of Equation (16) is the product of a PDF, $f(g)$, which is
always non-negative; a positive integer, $z/2$; and a conductance, $g$, which
is always non-negative. The denominator is the square of a positive real
number, which is always positive. As a result, the value of Equation (16) is
never negative, and is, in fact, always positive for any physically meaningful
PDF, $f(g)$. Moreover, we know that $F(0)=-1$ and that
$\raisebox{2.15277pt}{\scalebox{0.8}{$\displaystyle\lim_{g_{m}\rightarrow\infty}\;$}}F(g_{m})>0$.
So, $F(g_{m})$ in Equation (12) always has a positive slope and always has
exactly one intercept with the positive $g_{m}$ axis. This guarantees that the
EMT solution to Equation (13) is unique. Therefore, root finding methods
always converge (e.g. Newton’s method).

To solve for the root, Equation (12) and Equation (16) were written as Python
functions. These functions accept a PDF corresponding to $f(g)$ as an argument
and numerically calculate the integral function utilizing Python’s native
numerical integration capabilities (namely, the scipy.integrate package). The
Python codes implementing $F(g)$, $F^{\prime}(g)$, and the root finding
methods are available in Appendices A.1 and A.2.

It should be noted that Kirkpatrick’s effective-medium approach breaks down
near percolation thresholds but, to the best of our knowledge, this breakdown
has not been investigated in detail in the literature to-date.

###  IV.2 Pore-network modelling

To assess the accuracy of EMT solutions, we need to compare them to some
“empirical” value. To that end, PNM techniques are used to simulate resistor
networks and directly calculate the effective transport property of interest,
namely the effective conductance of the network. OpenPNM, an open-source PNM
package for Python, was used to perform these direct simulations [14]. Using
OpenPNM, square or cubic lattice resistor networks can be simulated, with
random conductances distributed according to specified PDFs. This leads to the
linear system of equations from Equation (8), which can be solved directly.
Technically, OpenPNM solves for an effective conductivity for the network.
This is easily converted into a conductance, however, by multiplying by the
length of the bonds.

For convenience, all simulated networks have a lattice constant of one, so the
effective conductivity and the effective conductance of the network have the
same numerical value. Networks examined were two and three-dimensional, with
an equal number of nodes $n$ in each direction. A third “quasi-two-
dimensional” case was also looked at, in which the network was very large in
two dimensions and small in the third. Specifically, the network had
dimensions of $n\times n\times 5$. This case is to represent thin porous
media, such as ionomer membranes in fuel cells. These membranes can be
millimetres or centimetres wide, but only nanometres or micrometres thick,
with channel lengths on the order of a few nanometres [16, 1]. Visualizations
of example 2D, 3D, and quasi-2D networks can be seen in Figure 3.

2D network.

3D network.

Quasi-2D network.

Figure 3: Examples of PNM networks with network voltages, $n=10$.

###  IV.3 Comparison of EMT and PNM

To compare EMT and PNM results, a simple procedure was followed for a given
network.

  1. 1.

Choose a PDF to describe network conductances.

  2. 2.

Calculate the EMT effective conductance, as described in Section IV.1.

  3. 3.

Calculate the effective conductance with PNM for a given $n$, using Monte
Carlo sampling to average the results.

  4. 4.

Repeat Step 3 on progressively larger networks.

  5. 5.

Calculate the relative difference (RD) between EMT and PNM effective
conductances as a function of network size $n$.

  6. 6.

Calculate the reduced standard deviation (RSD) of the PNM conductances as a
function of network size $n$.

These steps were repeated for a variety of different distributions $f(g)$ in
both two and three dimensions. The quasi-2D case was only examined for a
single PDF. Initially, 100 Monte Carlo samples were taken per distribution at
each $n$ while performing the PNM calculations. In the interest of
computational time required per simulation, this was subsequently reduced to
50 Monte Carlo samples, which proved to be a sufficient sample size.

RD for a given network was calculated according to

$\varepsilon=\frac{\lvert
g_{m}-\overline{g}_{\mathrm{eff}}\rvert}{g_{m}}\times 100\%,$ ((17))

where $g_{m}$ is the effective conductance according to EMT and
$\overline{g}_{\mathrm{eff}}$ is the averaged effective conductance from the
Monte Carlo samplings of the PNM calculations at a given $n$. Relative
standard deviation was calculated as

$\sigma_{\mathrm{RSD}}=\frac{\sigma_{g_{\mathrm{eff}}}}{\overline{g}_{\mathrm{eff}}}\times
100\%,$ ((18))

where $\sigma_{g_{\mathrm{eff}}}$ is the standard deviation of the Monte Carlo
samplings.

The goal of the above procedure is to provide guidance for future
computational research. The relationship between the RD and the RSD, and
network size $n$ provides scaling laws governing the convergence of PNM
results. This can serve as a rule of thumb, providing an estimate for when the
macroscopic limit is approached in PNM simulations, and when EMT can be
applied instead of a numerical solution of the PNM.

##  V Results

The procedure described in Section IV.3 was carried out for nine statistical
distributions of network conductance. These are summarized in Table 1. For
illustration, only the plots from the PDF referred to as “Weibull 2” will be
shown in this section. Plots from every distribution in Table 1 can be seen in
Appendix B.

The PDFs were chosen to examine the convergence behaviour between EMT and PNM
on a diverse array of network conductances. The arcsine distribution is
continuous and bimodal, with two sharp peaks in conductance value. The bimodal
sine distribution is also continuous and bimodal, but with less dramatic
peaks. The binary and trinary distributions were chosen to examine PDFs that
are both discrete and have analytic EMT solutions. The uniform distribution is
to examine when each conductance within a certain range has an equal
probability. The Weibull distributions were chosen to examine unimodal PDFs
with a peak at zero conductance (Weibull 1, equivalent to an exponential
distribution), near zero (Weibull 2), and away from zero (Weibull 3). The PDF
for Weibull 2 is shown in Figure 4. All other PDFs are visible in Figure 8.

Figure 4: PDF of Weibull 2 distribution with shape factor $k=1.5$. Table 1: Summary of probability density functions examined. Distribution |  PDF, $f(g)$
---|---
Arcsine | $\frac{1}{\pi\sqrt{g(1-g)}},~{}g\in[0,1]$
Bimodal sine | $\frac{\pi}{2}\begin{cases}\sin{(2\pi g)},&g\in[0,\frac{1}{2})\\\ -\sin{(2\pi g)},&g\in[\frac{1}{2},1]\\\ 0,&\text{otherwise}\end{cases}$
Binary 1 | $0.1\cdot\delta(g-0.3)+0.9\cdot\delta(g-0.6)$
Binary 2 | $0.5\cdot\delta(g-0.3)+0.5\cdot\delta(g-0.6)$
Trinary | $\frac{1}{3}\delta(g-0.25)+\frac{1}{3}\delta(g-0.5)+\frac{1}{3}\delta(g-0.75)$
Uniform | $\begin{cases}1,&g\in[0,1]\\\ 0,&\text{otherwise}\end{cases}$
Weibull 1 | $e^{-g},~{}g\in[0,\infty)$
Weibull 2 | $\frac{3}{2}g^{\frac{1}{2}}e^{-g^{\frac{3}{2}}},~{}g\in[0,\infty)$
Weibull 3 | $5g^{4}e^{-g^{5}},~{}g\in[0,\infty)$

The EMT solutions to Equation (13) on both 2D and 3D networks are plotted in
Figure 5. The intercept with the $g_{m}$ axis in the plots corresponds to the
effective conductance $g_{m}$ that replaces the infinite network. Solutions
for every PDF can be seen in Figure 10. The solutions for every network are
summarized in Table 2.

Table 2: EMT solutions to Equation (13). Distribution | 2D network | 3D network
---|---|---
Arcsine | 0.333 | 0.400
Bimodal sine | 0.414 | 0.445
Binary 1 | 0.561 | 0.564
Binary 2 | 0.424 | 0.433
Trinary | 0.455 | 0.475
Uniform | 0.398 | 0.437
Weibull 1 | 0.610 | 0.727
Weibull 2 | 0.705 | 0.771
Weibull 3 | 0.892 | 0.901

Figure 5: Equation (13) solutions for Weibull 2 distribution. Table 3: Least-square power law fits for RD and RSD, 2D and 3D networks. Distribution |  2D $a_{\mathrm{RD}}$ |  2D $b_{\mathrm{RD}}$ |  2D $a_{\mathrm{RSD}}$ |  2D $b_{\mathrm{RSD}}$ |  3D $a_{\mathrm{RD}}$ |  3D $b_{\mathrm{RD}}$ |  3D $a_{\mathrm{RSD}}$ |  3D $b_{\mathrm{RSD}}$
---|---|---|---|---|---|---|---|---
Arcsine | $\mathbf{-2.638}$ | $\mathbf{17947.503}$ | $-1.035$ | $154.325$ | $-1.088$ | $136.064$ | $-1.576$ | $132.446$
Bimodal sine | $-1.108$ | $149.310$ | $-0.946$ | $56.447$ | $-1.076$ | $134.850$ | $-1.472$ | $60.106$
Binary 1 | $-0.992$ | $99.080$ | $-1.014$ | $20.701$ | $-1.031$ | $114.829$ | $-1.248$ | $7.825$
Binary 2 | $-1.010$ | $107.622$ | $-0.911$ | $23.143$ | $-1.035$ | $117.497$ | $-1.793$ | $94.133$
Trinary | $-1.016$ | $109.052$ | $-0.997$ | $46.446$ | $\mathbf{-1.330}$ | $\mathbf{233.389}$ | $-1.625$ | $70.957$
Uniform | $-1.092$ | $129.761$ | $-1.020$ | $97.236$ | $-1.014$ | $106.608$ | $-1.571$ | $90.009$
Weibull 1 | $-0.961$ | $74.099$ | $-1.097$ | $172.700$ | $-1.000$ | $108.749$ | $-1.459$ | $89.196$
Weibull 2 | $-1.061$ | $129.554$ | $-1.041$ | $99.519$ | $-1.035$ | $117.506$ | $-1.508$ | $68.338$
Weibull 3 | $-1.038$ | $119.236$ | $-1.046$ | $31.159$ | $-1.043$ | $120.107$ | $-1.621$ | $36.443$

Least-squares fits of the RD between EMT and PNM (Equation (17)) and the PNM
Monte Carlo sampling RSD (Equation (18)) were made using Python’s native
scipy.optimize.curve_fit function. Fits were made to a power law,

$y(n)=bn^{a},$ ((19))

where $n$ is the size of the network, and $b$ and $a$ are the parameters of
the fit. The fits for the Weibull 2 distribution are in Figure 6. All other
fits can be seen in Figure 13. Fit parameter values are summarized in Table 3
for 2D and 3D networks. Errors in the fit parameters are presented in Appendix
C.

PNM and EMT RD.

PNM RSD.

PNM and EMT RD.

PNM RSD.

Figure 6: Power law fits to RD and RSD.

In addition to these comparisons between EMT and PNM on 2D and 3D networks,
the “quasi-2D” case was examined for the Weibull 2 distribution. Effective
conductance was calculated for transport along one of the long axes of $n$
nodes length, known as in-plane transport, and along the five-node axis, known
as through-plane transport. Both are scenarios outside of the scope considered
by Kirkpatrick’s EMT, which only considers infinite 2D or 3D networks. The
idealized quasi-2D case can be thought of as a 3D network that is infinite in
two dimensions, but is finite in the third dimension. Consequently, there is
no EMT solution to which the PNM results can be directly compared. However,
the RD between the PNM solution and both the 2D and 3D EMT solutions was still
calculated, in addition to the RSD. These can be seen in Figure 7. The fit
parameters for the in-plane transport are summarized in Table 4. Through-plane
transport reached a constant effective conductivity value after $n\sim 20$, so
fits were not performed on those data. Little information could be gained from
the latter analysis.

Table 4: Quasi-2D network fit parameters, in-plane transport. Distribution |  $a_{\mathrm{RD}}$, 2D EMT |  $b_{\mathrm{RD}}$, 2D EMT |  $a_{\mathrm{RD}}$, 3D EMT |  $b_{\mathrm{RD}}$, 3D EMT | $a_{\mathrm{RSD}}$ | $b_{\mathrm{RSD}}$
---|---|---|---|---|---|---
Weibull 2 | $-0.213$ | $24.328$ | $-1.755$ | $873.108$ | $-1.066$ | $41.651$

RD between 2D and quasi-2D.

RD between 3D and quasi-2D.

RSD of quasi-2D network.

Figure 7: Power law fits to quasi-2d relative differences and RSD, in-plane.

##  VI Discussion

###  VI.1 Relative difference

The problem posited in Section III was to examine the relationship between EMT
and PNM in modelling transport phenomena within composite materials and porous
media. Namely, what is the convergence behaviour between the two techniques as
the size of the system being examined grows? This was investigated for the
conductance distributions in Table 1, to see the effect of different
distributions on this behaviour. Two and three-dimensional networks were
examined. As evidenced from the fit parameters in Table 3, it appears that the
distribution of conductances has little effect on the convergence behaviour of
RD values between EMT and PNM. Indeed, the dimension of the network also seems
to have little effect.

For most of the PDFs examined, on both 2D and 3D cubic networks, the power law
convergence behaviour for RD follows the general relationship

$\varepsilon(n)\propto n^{-1}.$ ((20))

There is much agreement in the value of the exponential fit parameter $a$ in
Equation (19). There is less agreement in the scaling fit parameter $b$,
though the scaling parameter ultimately is less informative as to the
underlying trends being exhibited.

The RD is inversely proportional to the network size. This provides a handy
rule of thumb. If one would like their PNM solutions to exhibit a certain RD
from the theoretically exact EMT value for infinitely large systems, one can
use Equation (20) to determine approximately how large the network must be.
For example, if one desires an RD of $2\,\%$, then for most conductivity
distributions a network size, $n$, of

$n\approx\left(\frac{2\,\%}{100\,\%}\right)^{-1}=50$ ((21))

would be approximately sufficient, where dividing by $100\,\%$ accounts for
the $b$ scaling parameter, which is of order $100$ for almost all simulations.
This rule of thumb is an aid in determining appropriate domain sizes for
direct PNM calculations for continuum macroscale behaviour to emerge. Though
these calculations were specifically performed for electrical conductance,
transport behaviours of a variety of phenomena behave in analogous ways, as
discussed in Section II. It can be expected that similar trends will hold for
other transport problems.

There are two exceptions to this trend in our simulations, whose values are
bolded in Table 3. The most notable is the arcsine distribution on a 2D
network. Its exponent fit parameter is $a=-2.638$ and scaling fit parameter is
$b=17947.504$, very different from every other distribution. This distribution
has a large number of bonds in the network having a conductance at or very
near to zero. It is analogous to a softened binary distribution at zero and
one with equal probability. Such a binary distribution is exactly the
percolation threshold on 2D networks, where transport through the network
ceases entirely. The reduced connectivity in a square lattice (four nearest-
neighbour bonds, as opposed to six in a cubic lattice) places the arcsine
distribution very close to this threshold, which likely leads to atypical
transport behaviour within the network. As mentioned above, Kirkpatrick’s
method breaks down near the percolation method and the arcsine example
supports that notion.

The arcsine distribution is an unusual PDF and is not often used in the study
of composite media. The fact that it does not follow the general trends of
other, more conventional PDFs, demonstrates that great care should be taken
when applying such atypical models.

The other exception is in the trinary distribution on a 3D network. Its
exponent fit parameter is $a=-1.330$. This is not so outlandish as the arcsine
power law, but it remains significantly different from the other distributions
examined. Its scaling fit parameter is $b=233.389$, within an order of
magnitude of the other distributions, but still roughly twice as large. It is
unclear why this particular distribution differs from the overall trend. It is
seldom physically accurate to apply discrete PDFs to the properties of
composite media, so this discrepancy can perhaps be overlooked. Nevertheless,
the behaviour of this particular distribution remains an open question in this
analysis. Further studies into the convergence behaviour between PNM and EMT
on a greater variety of discrete PDFs could provide further insight.

###  VI.2 Reduced standard deviation

Much like the RD convergence discussed in Section VI.1, the exponential fit
parameter of the RSD values of the PNM Monte Carlo sampled conductances
appeared to be independent of the PDF used to distribute the conductances.
However, there was greater discrepancy between the results from 2D and 3D
networks, with 3D networks reducing their RSD values more quickly. This,
however, is to be expected. For a given network size $n$, a 3D network will
contain more bonds and, hence, a larger sample volume. Again, these results
can be applied more generally to other transport phenomena within porous
media. Additionally, there was greater discrepancy in the fit parameter $b$
between different distributions.

On 2D networks, the exponential fit parameter $a$ of the power law governing
the relationship between RSD and network size was similar to Equation (20);
inversely proportional. Consequently, the rule of thumb used to approximate
the network size required to achieve a desired RSD is the same as Equation
(21).

Discrete or narrow distributions (such as Weibull 3) have smaller $b$ fit
parameters, whereas wider distributions have larger $b$ fit parameters. This
indicates that the expected RSD is smaller for narrower distributions or ones
with more limited available conductance values, and it is larger for wide and
continuous distributions.

On 3D networks, the approximate power law is

$\sigma_{\mathrm{RSD}}\propto n^{-1.5},$ ((22))

which provides a good rule of thumb to estimate required network size for a
desired RSD. This simplified relationship neglects the influence of the $b$
parameter. Unless the distribution has been examined ahead of time, the $b$
parameter cannot be easily taken into account. However, the general power-law
behaviour of the networks as described by the $a$ parameter will still hold.

Again, let us take the example of $2\,\%$ as the desired value and a value of
$b$ of order $100$. This would provide an estimate for an appropriate network
size of

$n\approx\left(\frac{2\,\%}{100\,\%}\right)^{-0.667}\sim 14.$ ((23))

This is of the right magnitude but somewhat smaller than results found from
PNM simulations undertaken by Berg and Nadon [3]. There, it was determined
that a 3D network size of $n=24$ was needed to get the RSD values below
$2\,\%$. However, the PNM treatment in that study was significantly more
complex than undertaken here, which would be reflected by a value of $b$ that
deviates from $100$. The general trends found in this study can still provide
a useful starting point when determining the appropriate domain sizes for
simulation.

###  VI.3 Quasi-two-dimensional case

The final case examined is the quasi-2D network. As stated in Section IV.2,
this is defined as a network with dimensions of $n\times n\times 5$, imitating
thin porous media such as ionomer membranes. Only the Weibull 2 distribution
was examined on this network. The RSD behaviour on this network for in-plane
and through-plane transport was very similar to that displayed by true 2D
networks, following an approximately inverse relationship between
$\sigma_{\mathrm{RSD}}$ and $n$.

However, the results from PNM for this network did not conform with the
expected EMT results for either 2D or 3D lattices, as can be seen in Figures 7
and 7 for in-plane transport. Instead, as the network size grew, the
$\overline{g}_{\mathrm{eff}}$ approached a value distinct from either solution
(approximately 0.767 for in-plane and 0.995 for through-plane). This
underscores a very important detail when it comes to applying EMT. These
theories are developed to describe idealized effective media. They only
strictly apply in the continuum limit, when microscopic properties are
negligible and the domains are infinite. To use these results on real
materials, it must be assumed that the domains being examined are large enough
to be effectively infinite. Indeed, this is a reasonable approximation for
many fields of research, and EMT can be applied reliably.

However, as the methods of EMT and PNM begin to be applied to model phenomena
on micro- and nanoscales, these approximations can no longer be reasonably
held. Porous thin films used in nanotechnologies are three-dimensional.
However, they are significantly smaller in one dimension than in the other
two, but not to the extent that they can be treated as simply two-dimensional.
The utmost care must be used when attempting to apply EMT techniques in such
fields. Often, EMT will not be suitable.

##  VII Conclusion

The goal of this project is to provide insight into the relationship between
the macroscopic approaches of EMT and the discretized approaches of PNM in
regards to transport through composite and porous media. The specific case
examined in this treatment is electrical conductance, though the trends
discovered here can be generalized to other transport phenomena.

To that end, comparisons are made between predicted effective conductances
from EMT and directly calculated effective conductances from PNM on 2D and 3D
networks. Ultimately, from the analysis of the relative difference (RD)
between these values on progressively larger networks, useful rules of thumb
were developed governing this relationship. From the power law fits made to
the resulting data, it was determined that the RD is approximately inversely
proportional to the size of the network being examined independent of the
chosen PDF for conductance. This relationship held for both 2D and 3D
networks. This can be used as a guide for future research into these areas, to
determine what domain size is sufficient in order for macroscopic EMT models
to accurately replace PNM techniques. It can also be used to determine how
large of a network is required in order for PNM to model bulk transport
properties. All that is necessary is to decide what is an adequate threshold
value for RD for the system under investigation.

There are two caveats, however, to this rule of thumb. It breaks down when
examining the arcsine distribution on 2D networks and – to a lesser extent –
the discrete trinary distribution on 3D networks. This underscores that,
ultimately, the relationship is most useful as a guide but does not guarantee
a particular result. Caution is still required, especially when applying
unorthodox distributions to a particular medium. Further insight could be
gained by examining discrete distributions in more detail. Additionally, as
the EMT being examined was only applicable to square and cubic lattices, these
results will not necessarily apply to more complex network geometries. Also,
nonlinear transport properties (i.e. conductances) may have an impact on the
scaling law, which could be a potential subject of future research.
Nonetheless, the underlying convergence behaviour remains a useful tool.

The reduced standard deviation (RSD) of the Monte Carlo sampled, effective
conductance values from PNM were also examined on 2D and 3D networks, as a
function of network size. Again, it was found that a power law fit described
the convergence behaviour, though dimensionality did affect the result in this
case. On 2D networks, RSD was found to be inversely proportional to network
size. On 3D networks, the power law was closer to $-1.5$; 3D networks converge
more quickly as network size increases. The faster convergence behaviour of 3D
networks was not unexpected, and the results were in-line with RSD values from
another PNM study of transport within porous media [3].

When these techniques were applied to the quasi-2D case, however, the clean
convergence behaviour between EMT and PNM broke down. EMT applies to a
continuous, infinite medium, whether that medium be in two or three
dimensions. The quasi-2D case, however, is analogous to a three-dimensional
network that is infinite in two dimensions, but finite in the third. This can
be thought of as a proxy for thin porous media. As the network size grew
larger, the effective conductance did not converge to either the 2D or 3D
predictions from EMT. Assumptions that continuum methods can be applied to
such a material do not necessarily hold. This underscores that much caution
must be used if one wishes to utilize EMT techniques in analyzing the in-plane
and through-plane behaviour of thin porous media. A closer analysis of the
convergence scaling laws for quasi-2D systems forms another subject for future
research.

## References

  * Allen et al. [2015] F. I. Allen, L. R. Comolli, A. Kusoglu, M. A. Modestino, A. M. Minor, and A. Z. Weber.  Morphology of hydrated as-cast Nafion revealed through cryo electron tomography.  _ACS Macro Letters_ , 4(1):1–5, 2015.  doi: 10.1021/mz500606h.
  * Baschuk and Li [2000] J. J. Baschuk and X. Li.  Modelling of polymer electrolyte membrane fuel cells with variable degrees of water flooding.  _Journal of Power Sources_ , 86(1):181–196, 2000\.  doi: 10.1016/S0378-7753(99)00426-7.
  * Berg and Nadon [2021] P. Berg and P. Nadon.  Random pore-network model for polymer electrolyte membranes.  _Soft Matter_ , 17:5907–5920, 2021.  doi: 10.1039/D0SM02212H.
  * Bernasconi and Wiesmann [1976] J. Bernasconi and H. J. Wiesmann.  Effective-medium theories for site-disordered resistance networks.  _Physical Review B_ , 13(3):1131–1139, 1976.  doi: 10.1103/PhysRevB.13.1131.
  * Bruggeman [1935] D. A. G. Bruggeman.  Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen.  _Annalen der Physik_ , 416(7):636–664, 1935.  doi: 10.1002/andp.19354160705.
  * Burganos and Sotirchos [1987] V. N. Burganos and S. V. Sotirchos.  Diffusion in pore networks: Effective medium theory and smooth field approximation.  _AIChE Journal_ , 33(10):1678–1689, 1987.  doi: 10.1002/aic.690331011.
  * Choy [2015] T. C. Choy.  _Effective Medium Theory_ , Chapter 1. Essentials, pages 1–26.  Oxford University Press, Oxford, 2 edition, 2015.  doi: 10.1093/acprof:oso/9780198705093.003.0001.
  * Darcy [1856] H. Darcy.  _Les fontaines publiques de la ville de Dijon_.  Victor Dalmont, 1856.  Google-Books-ID: DOwbgyt_MzQC.
  * Das et al. [2010] P. K. Das, X. Li, and Z.-S. Liu.  Effective transport coefficients in PEM fuel cell catalyst and gas diffusion layers: Beyond Bruggeman approximation.  _Applied Energy_ , 87(9):2785–2796, 2010.  doi: 10.1016/j.apenergy.2009.05.006.
  * Dullien [1975] F. A. L. Dullien.  New network permeability model of porous media.  _AIChE Journal_ , 21(2):299–307, 1975.  doi: 10.1002/aic.690210211.
  * Eikerling et al. [1997] M. Eikerling, A. A. Kornyshev, and U. Stimming.  Electrophysical properties of polymer electrolyte membranes: A random network model.  _The Journal of Physical Chemistry B_ , 101(50):10807–10820, 1997.  doi: 10.1021/jp972288t.
  * Fatt [1956] I. Fatt.  The network model of porous media.  _Transactions of the AIME_ , 207(01):144–181, 1956.  doi: 10.2118/574-G.
  * Garnett [1904] J. C. M. Garnett.  XII. Colours in metal glasses and in metallic films.  _Philosophical Transactions of the Royal Society A_ , 203:385–420, 1904.  doi: 10.1098/rsta.1904.0024.
  * Gostick et al. [2016] J. Gostick, M. Aghighi, J. Hinebaugh, T. Tranter, M. A. Hoeh, H. Day, B. Spellacy, M. H. Sharqawy, A. Bazylak, A. Burns, W. Lehnert, and A. Putz.  OpenPNM: A pore network modeling package.  _Computing in Science & Engineering_, 18(4):60–74, 2016.  doi: 10.1109/MCSE.2016.49.
  * Kirkpatrick [1973] S. Kirkpatrick.  Percolation and conduction.  _Reviews of Modern Physics_ , 45(4):574–588, 1973\.  doi: 10.1103/RevModPhys.45.574.
  * Kusoglu and Weber [2017] A. Kusoglu and A. Z. Weber.  New insights into perfluorinated sulfonic-acid ionomers.  _Chemical Reviews_ , 117(3):987–1104, 2017.  doi: 10.1021/acs.chemrev.6b00159.
  * Larson et al. [1981] R. Larson, L. Scriven, and H. Davis.  Percolation theory of two phase flow in porous media.  _Chemical Engineering Science_ , 36(1):57–73, 1981.  doi: 10.1016/0009-2509(81)80048-6.
  * Lorentz [1880] H. A. Lorentz.  Ueber die Beziehung zwischen der Fortpflanzungsgeschwindigkeit des Lichtes und der Körperdichte.  _Annalen der Physik_ , 245(4):641–665, 1880.  doi: https://doi.org/10.1002/andp.18802450406.
  * Lorenz [1880] L. Lorenz.  Ueber die Refractionsconstante.  _Annalen der Physik_ , 247(9):70–103, 1880.  doi: https://doi.org/10.1002/andp.18802470905.
  * Note [1] Note1.  “For simplicity we shall consider [the network] to be made up of a set of equal conductances, $g_{m}$, connecting nearest neighbors on the cubic mesh,” Kirkpatrick [15].
  * Payatakes et al. [1973] A. C. Payatakes, C. Tien, and R. M. Turian.  A new model for granular porous media: Part I. Model formulation.  _AIChE Journal_ , 19(1):58–67, 1973.  doi: 10.1002/aic.690190110.
  * Rayleigh [1892] L. Rayleigh.  LVI. On the influence of obstacles arranged in rectangular order upon the properties of a medium.  _The Philosophical Magazine_ , 34(211):481–502, 1892.  doi: 10.1080/14786449208620364.
  * Tjaden et al. [2016] B. Tjaden, S. J. Cooper, D. J. Brett, D. Kramer, and P. R. Shearing.  On the origin and application of the Bruggeman correlation for analysing transport phenomena in electrochemical systems.  _Current Opinion in Chemical Engineering_ , 12:44–51, 2016\.  doi: 10.1016/j.coche.2016.02.006.
  * Vadakkepatt et al. [2015] A. Vadakkepatt, B. Trembacki, S. R. Mathur, and J. Y. Murthy.  Bruggeman’s exponents for effective thermal conductivity of lithium-ion battery electrodes.  _Journal of The Electrochemical Society_ , 163(2):A119, 2015.  doi: 10.1149/2.0151602jes.
  * Xiong et al. [2016] Q. Xiong, T. G. Baychev, and A. P. Jivkov.  Review of pore network modelling of porous media: Experimental characterisations, network constructions and applications to reactive transport.  _Journal of Contaminant Hydrology_ , 192:101–117, 2016\.  doi: 10.1016/j.jconhyd.2016.07.002.

##  Appendix A Code

###  A.1 Kirkpatrick EMT and derivative

The following code was written to calculate the value of the EMT function
$F(g_{m})$ (Equation (12)) and its derivative (Equation (16)), for numerical
root-finding.

⬇

import numpy as np

import scipy.integrate as integrate

\pardef kirkpatrick_integral(f, gm, z):

"""

Numerically calculates the integral described in Kirkpatrick (1973)

\parParameters

———-

f : FUNCTION

Probability density function describing the distribution of

conductances within the network.

gm : FLOAT

The conductivity of the effective medium being considered.

z : INT

Describes the dimensions of the network. 4 corresponds to a 2D square

lattice, each node having 4 nearest neighbours. 6 corresponds to a

3D cubic lattice, each node having 6 nearest neighbours.

\parReturns

——-

F : FLOAT

The numerically integrated value returned by the integral.

\par"""

integrand = lambda g: (f(g)*(gm - g))/(g + (0.5*z - 1)*gm)

\parF, err = integrate.quad(integrand, 0, np.inf)

return F

\pardef kirkpatrick_derivative(f, gm, z):

"""

Numerically calculates the derivative of the integral equation described

in Kirkpatrick (1973) wrt to d/d(gm)

\parParameters

———-

f : FUNCTION

Probability density function describing the distribution of

conductances within the network.

gm : FLOAT

The conductivity of the effective medium being considered.

z : INT

Describes the dimensions of the network. 4 corresponds to a 2D square

lattice, each node having 4 nearest neighbours. 6 corresponds to a

3D cubic lattice, each node having 6 nearest neighbours.

\parReturns

——-

Fprime : FLOAT

The numerically integrated value of the derivative.

\par"""

integrand = lambda g: (f(g)*0.5*z*g)/(g + (0.5*z - 1)*gm)**2

\parFprime, err = integrate.quad(integrand, 0, np.inf)

return Fprime

###  A.2 Root finders

The following functions were used to numerically calculate the root of the EMT
function $F(g_{m})$, using bisection or Newton methods.

⬇

def bisection(f,a,b,N):

”’Approximate solution of f(x)=0 on interval [a,b] by bisection method.

\parParameters

———-

f : function

The function for which we are trying to approximate a solution f(x)=0.

a,b : numbers

The interval in which to search for a solution. The function returns

None if f(a)*f(b) >= 0 since a solution is not guaranteed.

N : (positive) integer

The number of iterations to implement.

\parReturns

——-

x_N : number

The midpoint of the Nth interval computed by the bisection method. The

initial interval [a_0,b_0] is given by [a,b]. If f(m_n) == 0 for some

midpoint m_n = (a_n + b_n)/2, then the function returns this solution.

If all signs of values f(a_n), f(b_n) and f(m_n) are the same at any

iteration, the bisection method fails and return None.

\parExamples

——–

>>> f = lambda x: x**2 - x - 1

>>> bisection(f,1,2,25)

1.618033990263939

>>> f = lambda x: (2*x - 1)*(x - 3)

>>> bisection(f,0,1,10)

0.5

”’

if f(a)*f(b) >= 0:

print("Bisection method fails.")

return None

a_n = a

b_n = b

for n in range(1,N+1):

m_n = (a_n + b_n)/2

f_m_n = f(m_n)

if f(a_n)*f_m_n < 0:

a_n = a_n

b_n = m_n

elif f(b_n)*f_m_n < 0:

a_n = m_n

b_n = b_n

elif f_m_n == 0:

print("Found exact solution.")

return m_n

else:

print("Bisection method fails.")

return None

return (a_n + b_n)/2

\pardef newton(f,Df,x0,epsilon,max_iter):

”’Approximate solution of f(x)=0 by Newton’s method.

\parParameters

———-

f : function

Function for which we are searching for a solution f(x)=0.

Df : function

Derivative of f(x).

x0 : number

Initial guess for a solution f(x)=0.

epsilon : number

Stopping criteria is abs(f(x)) < epsilon.

max_iter : integer

Maximum number of iterations of Newton’s method.

\parReturns

——-

xn : number

Implement Newton’s method: compute the linear approximation

of f(x) at xn and find x intercept by the formula

x = xn - f(xn)/Df(xn)

Continue until abs(f(xn)) < epsilon and return xn.

If Df(xn) == 0, return None. If the number of iterations

exceeds max_iter, then return None.

”’

xn = x0

for n in range(0,max_iter):

fxn = f(xn)

if abs(fxn) < epsilon:

print(’Found solution after’,n,’iterations.’)

return xn

Dfxn = Df(xn)

if Dfxn == 0:

print(’Zero derivative. No solution found.’)

return None

xn = xn - fxn/Dfxn

print(’Exceeded maximum iterations. No solution found.’)

return None

###  A.3 PNM Monte Carlo sampler

The following is an example code used to carry out Monte Carlo sampling using
PNM techniques. It is currently configured for examining 3D networks while
taking samples from a binary distribution. OpenPNM contains several native
statistical distributions. The binary distribution is an example of a custom-
written distribution.

⬇

import numpy as np

import openpnm as op

import openpnm.models as mods

import custom_models as cm

\par# Array of network sizes

network_sizes = np.array((5, 10, 15, 20, 25, 30, 35, 40, 43, 45))

# Arrays to store the mean and rsd values of geff

mean_effective_conductivity = np.zeros(len(network_sizes))

rsd_effective_conductivity = np.zeros(len(network_sizes))

\par# Number of Monte Carlo samples to take

M = 50

# Connectivity of the cubic lattice

connectivity = 6

\parfor j in range(0, len(network_sizes)):

# The size of the network

n = network_sizes[j]

x = n

y = n

# 2D or 3D network; 1 if 2D, n if 3D

z = n

\par# Lattice constant of the network; set to 1, then conductivity and

# conductance have the same numerical value

Lc = 1

\par# Array to store the effective conductivities of calculated on the random

# networks

effective_conductivity = np.zeros(M)

\par# Setup the network

pn = op.network.Cubic(shape=[x,y,z], spacing=Lc, connectivity=connectivity)

Ps = pn.pores()

Ts = pn.throats()

geo = op.geometry.GenericGeometry(network=pn, pores=Ps, throats=Ts)

\parfor i in range(0, M):

geo.add_model(propname=’throat.seed’,

model=mods.geometry.throat_seed.random)

geo.add_model(propname=’pore.diameter’,

model=mods.misc.constant,

value=0)

geo.add_model(propname=’throat.diameter’,

model=mods.misc.constant,

value=Lc)

geo.add_model(propname=’pore.area’,

model=mods.geometry.pore_area.cube)

geo.add_model(propname=’throat.area’,

model=mods.geometry.throat_area.cuboid)

geo.add_model(propname=’throat.endpoints’,

model=mods.geometry.throat_endpoints.cubic_pores)

geo.add_model(propname=’throat.conduit_lengths’,

model=mods.geometry.throat_length.conduit_lengths)

\par# Create phase

phase = op.phases.GenericPhase(network=pn)

phase[’throat.seed’] = geo[’throat.seed’]

\par# Assign conductivities according to statistical distribution

phase.add_model(propname=’throat.electrical_conductivity’,

model = cm.binary_distribution, seeds=’throat.seed’,

p=0.1, alpha=0.3, beta=0.6)

# Pores need a value, though it does not affect the geff calculation

phase[’pore.electrical_conductivity’] = np.zeros(len(Ps))

\par# Add physics objects

phys = op.physics.GenericPhysics(network=pn, phase=phase, geometry=geo)

phys.add_model(propname=’throat.electrical_conductance’,

model=mods.physics.electrical_conductance.series_resistors)

\par# Add algorithm for calculated net conductance

ohm = op.algorithms.OhmicConduction(network=pn)

ohm.setup(phase=phase)

ohm.set_value_BC(pores=pn.pores(’left’), values=0)

ohm.set_value_BC(pores=pn.pores(’right’), values=100)

ohm.run()

geff = ohm.calc_effective_conductivity(inlets=pn.pores(’right’),

outlets=pn.pores(’left’),

domain_area=x*z,

domain_length=y)

effective_conductivity[i] = geff

\parmean_geff = np.mean(effective_conductivity)

mean_effective_conductivity[j] = mean_geff

std_geff = np.std(effective_conductivity)

rsd_effective_conductivity[j] = std_geff*100/mean_geff

\parnp.savetxt(’binary_p0p1_3d_mean_g.csv’,

np.array((network_sizes, mean_effective_conductivity)),

delimiter=’,’)

np.savetxt(’binary_p0p1_3d_rsd_g.csv’,

np.array((network_sizes, rsd_effective_conductivity)),

delimiter=’,’)

##  Appendix B Additional plots

###  B.1 Probability density functions

Extreme bimodal values. Continuous bimodal values. 0.1–0.9 split between
values. Even split between values. Even split between values. Uniform
distribution.

Figure 8: Probability density functions (PDFs), $f(g)$, examined in our
analysis.

Unimodal peak near zero-conductance. Unimodal peak away from zero.

Figure 9: Probability density functions (PDFs), $f(g)$, examined in our
analysis.

###  B.2 EMT plots

Figure 10: Solutions of Equation (13) for PDFs listed in Table 1.

Figure 11: Solutions of Equation (13) for PDFs listed in Table 1.

Figure 12: Solutions of Equation (13) for PDFs listed in Table 1.

###  B.3 Power law fit plots: RD and RSD

Figure 13: Power law fits to relative difference (RD) and reduced standard
deviation (RSD) of PNM solution, as a function of $n$.

Figure 14: Power law fits to relative difference (RD) and reduced standard
deviation (RSD) of PNM solution, as a function of $n$.

Figure 15: Power law fits to relative difference (RD) and reduced standard
deviation (RSD) of PNM solution, as a function of $n$.

Figure 16: Power law fits to relative difference (RD) and reduced standard
deviation (RSD) of PNM solution, as a function of $n$.

Figure 17: Power law fits to relative difference (RD) and reduced standard
deviation (RSD) of PNM solution, as a function of $n$.

Figure 18: Power law fits to relative difference (RD) and reduced standard
deviation (RSD) of PNM solution, as a function of $n$.

##  Appendix C Least-square fit parameter errors

Table 5: Least-square fit errors (%) for fit parameters in Table 3. Distribution |  2D $a_{\mathrm{RD}}$ |  2D $b_{\mathrm{RD}}$ |  2D $a_{\mathrm{RSD}}$ |  2D $b_{\mathrm{RSD}}$ |  3D $a_{\mathrm{RD}}$ |  3D $b_{\mathrm{RD}}$ |  3D $a_{\mathrm{RSD}}$ |  3D $b_{\mathrm{RSD}}$
---|---|---|---|---|---|---|---|---
Arcsine | 6.922 | 60.533 | 1.808 | 7.049 | 1.360 | 4.656 | 4.543 | 21.478
Bimodal sine | 2.745 | 11.329 | 3.126 | 11.298 | 1.080 | 3.685 | 2.245 | 10.053
Binary 1 | 1.391 | 5.232 | 4.361 | 16.710 | 0.466 | 1.522 | 5.346 | 20.638
Binary 2 | 2.931 | 11.192 | 6.241 | 21.839 | 0.798 | 2.612 | 1.361 | 7.196
Trinary | 3.665 | 14.067 | 5.583 | 21.084 | 1.757 | 7.171 | 3.640 | 17.674
Uniform | 5.495 | 22.402 | 4.421 | 17.024 | 0.961 | 3.090 | 4.090 | 19.290
Weibull 1 | 8.423 | 30.852 | 5.394 | 22.076 | 1.469 | 4.700 | 1.267 | 5.628
Weibull 2 | 5.738 | 22.832 | 2.734 | 10.710 | 2.981 | 9.765 | 7.632 | 34.744
Weibull 3 | 1.102 | 4.307 | 1.941 | 7.633 | 0.364 | 1.199 | 3.211 | 15.564
Table 6: Least-square fit errors (%) for fit parameters in Table 4. Distribution |  $a_{\mathrm{RD}}$, 2D EMT |  $b_{\mathrm{RD}}$, 2D EMT |  $a_{\mathrm{RD}}$, 3D EMT |  $b_{\mathrm{RD}}$, 3D EMT | $a_{\mathrm{RSD}}$ | $b_{\mathrm{RSD}}$
---|---|---|---|---|---|---
Weibull 2 | 5.436 | 4.961 | 10.100 | 61.340 | 3.087 | 12.332
  *[EMT]: effective medium theory
  *[PNM]: pore-network model
  *[PDF]: probability density function
  *[RD]: relative difference
  *[RSD]: reduced standard deviation