Voronoi tessellation-based algorithm for determining rigorously defined classical and generalized geometric pore size distributions

The geometric pore size distribution (PSD) $P\left(r\right)$ as function of pore radius $r$ is an important characteristic of porous structures, including particle-based systems, because it allows us to analyze adsorption behavior, the strength of materials, etc. Multiple definitions and corresponding algorithms, particularly in the context of computational approaches, exist that aim at calculating a PSD, often without mentioning the employed definition and therefore leading to qualitatively very different and apparently incompatible results. Here, we analyze the differences between the PSDs introduced by Torquato et al. and the more widely accepted one provided by Gelb and Gubbins, here denoted as T-PSD and G-PSD, respectively, and provide rigorous mathematical definitions that allow us to quantify the qualitative differences. We then extend G-PSD to incorporate the ideas of coating, which is significant for nanoparticle-based systems, and of finite probe particles, which is crucial to micro and mesoporous particles. We derive how the extended and classical versions are interrelated and how to calculate them properly. We next analyze various numerical approaches used to calculate classical G-PSDs and may be used to calculate the generalized G-PSD. To this end, we propose a simple yet sufficiently complicated benchmark for which we calculate the different PSDs analytically. This approach allows us to completely rule out a recently proposed algorithm based on radical Voronoi tessellation. Instead, we find and prove that the output of a grid-free classical Voronoi tessellation, namely, the properties of its triangulated faces, can be used to formulate an algorithm, which is capable of calculating the generalized G-PSD for a system of monodisperse spherical particles (or points) to any precision, using analytical expressions. The Voronoi-based algorithm developed and provided here has optimal scaling behavior and outperforms grid-based approaches.

Figure 1

Multidisk setup. Construction, shape and size of the volume $V\left(r\right)=V\left(r\right|0)$ for a 2D system consisting of 20 noncoated circles of radius $r_{\circ }=0.1$ (filled circles) in a periodic box of side length $L=1$ for two selected cases: (a) $r=0.1$, and (b) $r=0.3$. The first column shows the euclidean distance map (EDM, colored) constructed from the particles center positions (dark blue), as well as a single contour line (black) at altitude $r+r_{\circ }$. The second column highlights the region (green) enclosed by the EDM contour line; it carries all possible center coordinates of circular pores of radius $r$. The last column shows the volume $V\left(r\right|0)$ (green).

Figure 2

Enclosing spheres for a material with just two solid 3D spheres. (a) Largest enclosing sphere (LES, black empty circle, actually representing a sphere) of radius $r\left(\mathbf{p}\right)$, centered at $\mathbf{c}$ (cyan), containing the probe point particle centered at $\mathbf{p}$ (green). (b) Largest probe enclosing sphere (LPES, black empty circle) of radius $r(\mathbf{p};r_p)$ centered at another $\mathbf{c}$, containing the probe sphere (radius $r_p>0$) centered at $\mathbf{p}$. (c) LES of radius $r\left(\mathbf{p}\right|r_c)$ within the $r_c$-coated material. (d) A hypothetical material shell of size $r_p$ can be used to calculate $r(\mathbf{p};r_p)-r_p=r\left(\mathbf{p}\right|r_p)$, cf. Eq. (4), and to reduce the LPES problem to the setup shown in panel (c), where $r_c$ is replaced by an effective radius $r_{\text{eff}}=r_c+r_p$, cf. Eq. (1). The distance $R\left(\mathbf{c}\right)$ between $\mathbf{c}$ and the center of the nearest material particle (radius $r_{\circ }$) is different in panels (a)–(c). For sufficiently large probe particles, there are regions in the void space (coating overlap) that cannot host a finitely sized probe particle. Note that the surface of the probe particle must not touch the LES or LPES in general, if there are more than two material spheres. This will become obvious later below.

Figure 3

Multidisk setup. Construction, shape and size of the volume $V\left(r\right|r_{\text{eff}})$ for the system treated already by Fig. 1, here for two selected cases which share their sum $r+r_{\text{eff}}=0.2$: (a) $r=0.2, r_{\text{eff}}=0$, (b) $r=r_{\text{eff}}=r_c+r_p=0.1$ (for example $r_c=0.1$ and $r_p=0$ or $r_c=0$ and $r_p=0.1$). The first column shows the euclidean distance map (EDM, colored) constructed from the particles center positions (dark blue), as well as a single contour line (black) at altitude $r+r_s$. The second column highlights the region (green) enclosed by the EDM contour line; it carries all possible center coordinates of circular pores of radius $r$. The black circles in panel (b) mark the outer surface of $r_{\text{eff}}$-coated particles. The last column shows the volume $V\left(r\right|r_{\text{eff}})$ (green), i.e., $V(0.2,0)$ in panel (a) and $V(0.1,0.1)$ in panel (b). All variants of PSD's are encoded in the two-parametric $V\left(r\right|r_{\text{eff}})$. The numerical result for $V\left(r\right|r_{\text{eff}})$ over the whole range of semipositive $r$ and $r_{\text{eff}}$ for this setup is shown in Fig. 6, along with the corresponding results for the generalized G-PSD.

Figure 4

Relevant PSD volumes for the triangle and our benchmark setup. Drastic differences between the relevant (gray-shaded) areas used to define the G-PSD $P\left(r\right)$ of Gelb & Gubbins (G) [31] as well as the ${\mathrm{T}}_{\mathrm{c}}$-PSD and ${\mathrm{T}}_{\mathrm{s}}$-PSD of Torquato et al. [30]. While for the ${\mathrm{T}}_{\mathrm{c}}$-PSD (${\mathrm{T}}_{\mathrm{s}}$-PSD) the red circle touches at least one center (one surface) of a material particle, for the G-PSD it touches at least two surfaces of material particles. (a)–(c) Naked triangle serving as material. Its edges can be thought of consisting of infinitely many circles of radius $r_{\circ }=0$. (d)–(f) Our naked benchmark setup, consisting of a single material particle of radius $r_{\circ }$ subject to periodic boundary conditions, i.e., tangent to its periodic image. (a) Shaded area $V\left(r\right|0)$ accessible to a circle of radius $r$, relevant for the G-PSD. (b), (c) Shaded area $V\left(0\right|r-r_{\circ })=V\left(0\right|r)$ enclosed by $r$-coated triangle, relevant for both ${\mathrm{T}}_{\mathrm{c}}$-PSD and ${\mathrm{T}}_{\mathrm{s}}$-PSD, because $r_{\circ }=0$ for the triangle. (d) One quarter of the area $V\left(r\right|0)$ (gray-shaded) accessible to a circle of radius $r$, relevant for the G-PSD. (e) One quarter of the area $V\left(0\right|r-r_{\circ })$ (gray-shaded) created by the material circle coated by a shell of thickness $r-r_{\circ }$ (for this particular example, negatively coated), relevant for the ${\mathrm{T}}_{\mathrm{c}}$-PSD. (f) One quarter of the area $V\left(0\right|r)$ (gray-shaded) created by the material circle coated by a shell of thickness $r$, relevant for the ${\mathrm{T}}_{\mathrm{s}}$-PSD. While the ${\mathrm{T}}_{\mathrm{c}}$- and ${\mathrm{T}}_{\mathrm{s}}$-PSDs are both trivially calculated at minor computational effort, the efficient calculation of the G-PSD still poses a challenge, as these simple examples might already suggest.

Figure 5

The Voronoi-based algorithm locates the center $\mathbf{c}$ (blue point) and radius of the pore (black open circle) for a given point $\mathbf{p}$ (green) either on (a) a vertex or (b) an edge of the Voronoi diagram (black straight lines) corresponding to the material circles (black solid circles). We calculate the position of $\mathbf{c}$ analytically, $t\in \{0,1\}$ for vertex-centered pores, and $t\in (0,1)$ for edge-centered pores in Sec. 3d.

Figure 6

Multidisk setup. Results from the presented algorithm for the configuration inspected in Figs. 1 and 3. The former figures provided visual information on how to obtain 4 selected points (marked by red bullets) in the present quantitative plots (a) $V\left(r\right|r_{\text{eff}})$, (b) the cumulative G-PSD $P_{\mathrm{cum}}\left(r\right|r_{\text{eff}})$, and (c) the G-PSD $P\left(r\right|r_{\text{eff}})$ defined in Eq. (8). Beyond a certain $r+r_{\text{eff}}$ (the radius of the LES), there is no contour line anymore in the left column of Fig. 1, and therefore vanishing $V\left(r\right|r_{\text{eff}})=0, P\left(r\right|r_{\text{eff}})=0$, and $P_{\mathrm{cum}}\left(r\right|r_{\text{eff}})=1$. For the uncoated system ($r_c=0$), the classical PSD $P\left(r\right)$ for point-like probe particles ($r_p=0$) is contained in panel (c) on the $x$ axis, i.e., at $r_{\text{eff}}=r_c+r_p=0$.

Figure 7

Naked (uncoated) benchmark setup. The setup consists of identical circles of radius $r_{\circ }=1$ or, equivalently, infinitely long and aligned cylinders, that just touch each other ($L=2r_{\circ }$). The colored filled circles represent two possible largest circles that can be inscribed within the (white) void space enclosed by the large circles. The blue circle with radius $r_{\max }$ (24), represents the circle with largest possible radius that can be inscribed within this setup. This situation is realized in a periodic simulation cell (black square) of side length $L$, carrying a single circle of radius $r_{\circ }$.

Figure 8

Benchmark geometry. First quadrant of the $r_c$-coated benchmark setup (Fig. 7). Here we introduce radius $r$ of the circle centered at $(x,0)$, radius $r_p$ of the probe particle, and some interrelated angles $\alpha , \beta$, and $\gamma$. The shaded area $V\left(r\right|r_{\text{eff}})$ is the area that needs to be calculated as function of $r$, for given effective shell thickness, $r_{\text{eff}}=r_c+r_p$ to obtain the generalized G-PSD. The classical G-PSD is adsorbed by the special case $r_{\text{eff}}=0$.

Figure 9

Benchmark (analytic results). Mean pore radius $\langle r\rangle$ for the benchmark configuration, obtained analytically. (a) $\langle r\rangle$ versus $r_c$ for various probe particle radii $r_p$. (b) $\langle r\rangle$ versus $r_c$ for various coating thickness values $r_c$. The sum $r_{\text{eff}}=r_c+r_p$ is limited by $r_{\max }$ for $r_{\text{eff}}=0$, which is $(\sqrt{2}-1)r_{\circ }\approx 0.41r_{\circ }$.

Figure 10

Benchmark. Analytic distribution functions for the benchmark, Eq. (27), all exactly reproduced with the numerical scheme of Sec. 3d. (a) $P_{\mathrm{cum}}(r;r_p|0)=P_{\mathrm{cum}}(r-r_p|r_p)$ for various probe radii $r_p$ in the naked setup. (b) G-PSD $P(r;r_p|0)$ corresponding to the cumulative G-PSD shown in panel (a). (c) $P_{\mathrm{cum}}(r;0|r_c)=P_{\mathrm{cum}}\left(r\right|r_c)$ for point-like probe particle at various coating thickness values $r_c$. (d) $P_{\mathrm{cum}}(r;r_p|r_c)=P_{\mathrm{cum}}(r-r_p|r_p+r_c)$ for probe particles of radius $r_p=0.1$ at various coating thickness values $r_c$.

Figure 11

Benchmark. The largest pore radius $r_{\max }$ is covering a significant portion $V\left(r_{\max }\right|r_{\text{eff}})$ of the total pore volume $V\left(0\right|r_{\text{eff}})$ for the $r_{\text{eff}}$-coated benchmark setup, leading to a jump in the cumulative G-PSD $P_{\mathrm{cum}}\left(r\right|r_{\text{eff}})$ from the value $P_{\mathrm{cum}}\left(r_{\max }\right|r_{\text{eff}})$ to unity for $r>r_{\max }$. The effect of $r_{\text{eff}}$ on this value is displayed here because it is an important characteristics of the whole benchmark G-PSD. Because $r_{\text{eff}}$ can only take values between 0 and $(\sqrt{2}-1)r_{\circ }$ for the benchmark (Fig. 8), we here plotted versus the scaled a dimensionless $r_{\text{eff}}$ that resides on the interval [0,1].

Figure 12

Triangle. For the equilateral triangle with unit side length and a straight edges ($r_{\circ }=0$), the volume $V\left(r\right|r_{\text{eff}})$ can be calculated analytically, cf. Eq. (A4) and leads to the G-PSD quantities (a) $P_{\mathrm{cum}}\left(r\right)=\left[4(9-\sqrt{3}\pi )r_1^2\right]/3$, (b) $P\left(r\right)=\left[8(9-\sqrt{3}\pi )r_1\right]/3$, while both T-PSDs are identical in that case, (a) ${\mathcal{P}}_{\mathrm{cum}}^{\mathrm{c}}\left(r\right)={\mathcal{P}}_{\mathrm{cum}}^{\mathrm{s}}\left(r\right)=4(\sqrt{3}-3r_1)r_1$, and (b) $\mathcal{P}\left(r\right)=4(\sqrt{3}-6r_1)$. This figure demonstrates the huge qualitative difference between G-PSD and T-PSDs. The same is true for our benchmark, cf. Eq. (31) versus Eq. (27) for $r_s=r_{\circ }$.

Figure 13

Benchmark. Comparison (quality and speed) between methods to calculate the cumulative G-PSD $P_{\mathrm{cum}}\left(r\right|r_{\text{eff}})$ for the benchmark setup with $r_{\text{eff}}=0$: (i) the Voronoi-based global analytical optimization developed in this work (Sec. 3d) with $M$ shots, (ii) regular grid-based brute force scan (Sec. 3a) with $M$ nodes. The deviation is a standard deviation calculated against the analytical solution (27). The computational effort is strictly linear in $M$ for the Voronoi-based version, and generally independent of $N$, while the grid-based approach is $\propto M^{5/3}$ for this case, and also independent of $N$. The local constrained optimization approach is almost linear in $M$, but in addition linear in $N$ except for systems with very small pores, where a neighbor list can remove the dependency on $N$ (Sec. 3b).

Figure 14

Multidisk. Comparison of convergence behavior for Voronoi-based and grid-based methods during the calculation of the cumulative G-PSD $P_{\mathrm{cum}}\left(r\right|r_{\text{eff}})$ for the multidisk $N=1000$ setup shown on the left. The right panel displays $P_{\mathrm{cum}}\left(r\right|0)$ versus $r$ for various choices of $M$, for both methods. While the grid-based result slowly converges against the Voronoi-based result, the latter remains basically unaffected by $M$ for sufficiently large $M\gg N$. The underlying reason is the absence of a grid for the Voronoi-based method. At the same time it runs more efficiently and does not suffer from memory limitations. The present result just confirms the trends already available from our benchmark (Fig. 13).