The Effect of Topology on Phase Behavior under Confinement

Abstract

This work studies how morphology (i.e., the shape of a structure) and topology (i.e., how different structures are connected) influence wall adsorption and capillary condensation under tight confinement. Numerical simulations based on classical density functional theory (cDFT) are run for a wide variety of geometries using both hard-sphere and Lennard-Jones fluids. These cDFT computations are compared to results obtained using the Minkowski functionals. It is found that the Minkowski functionals can provide a good description of the behavior of Lennard-Jones fluids down to small system sizes. In addition, through decomposition of the free energy, the Minkowski functionals provide a good framework to better understand what are the dominant contributions to the phase behavior of a system. Lastly, while studying the phase envelope shift as a function of the Minkowski functionals it is found that topology has a different effect depending on whether the phase transition under consideration is a continuous or a discrete (first-order) transition.

1. Introduction

Under tight confinement, a gas can form a condensed phase at a pressure below the bulk vapor pressure. This phenomenon is known as capillary condensation and has applications in many fields of science and engineering, including the storage of hydrogen carriers [1,2,3], battery technology [4], hydrocarbon extraction from unconventional reservoirs [5], and carbon dioxide sequestration [6]. Capillary condensation can have a large effect on transport properties, including effective diffusion coefficients [7], imbibition [8], and mass flow rates [9], and it is reported in the literature that both morphology (i.e., the shape of a structure) and topology (i.e., how different structures are connected) have a strong effect on the sorption of both sub- and supercritical fluids [2,10,11]. For ordered porous media, the relation between capillary condensation and geometry is well understood [12]; however, in practice, many porous media are disordered rather than ordered. To get a better understanding of capillary condensation in disordered porous media, simple geometries like cylinders, slit pores, ink bottles, and spheres [13,14] have been studied extensively. Unfortunately, topology is generally not considered in these model systems and, because there is no well-defined geometrical parameter space, the results of these studies are difficult to compare with each other. As a consequence, it is still not well understood how phase transitions in porous media are affected by topology [14,15,16,17,18,19].

In this work, we study the effect of topology on capillary condensation and wall adsorption under confinement (i.e., small pores) through the lens of the Minkowski functionals. The aim is to develop a description of phase behavior under confinement that is universally valid, independent of the details of the pore structure under consideration. The Minkowski functionals are a concept from integral geometry which not only characterize the morphology, but also the topology of spatial patterns [20]. They have been applied in a wide array of research areas including astronomy [21,22], statistical physics [23], phase behavior [24], granular materials [25,26], and porous media [27,28]. For a system in D dimensions, there are $D+1$ Minkowski functionals and in the case of a two-dimensional system, these functionals are related to the surface area, circumference, and signed curvature (i.e., the Euler characteristic) of the system [29]. A central component of the Minkowski functionals framework is Hadwiger’s theorem [30]. It states that in many cases, the free energy of a system, or any other functional, can be expressed as a linear combination of Minkowski functionals [30,31], thus decoupling the thermodynamical behavior of a system from its morphology and topology. This has two important consequences: (i) Characterization and classification; when the assumptions of Hadwiger’s theorem are fulfilled, the Minkowski functionals provide a complete spatial description of a physical system. This means that they provide a finite well defined geometrical parameter space that allows for the direct comparison of the results of different studies using different geometries. (ii) Scale translation; as mentioned above, the Minkowski functionals allow for the separation of the thermodynamics of a system from its geometry. This makes it possible to perform a simulation of a small system with known Minkowski functionals, compute the linear coefficients for various thermodynamic properties, and use these to derive the properties of a much larger system with known Minkowski functionals. In addition, this also implies that simple geometries can be used to represent complex geometries.

To study fluids under confinement, the Minkowski functionals have been used in combination with both hard-sphere fluids [31] and more complex force fields in conjunction with molecular dynamics [32]. Building on the work of König et al. [31], this study looks at the efficacy of using the Minkowski functionals to describe the behavior of Lennard-Jones fluids under confinement. A Lennard-Jones potential is more realistic than a hard-sphere potential for many applications. However, the interaction length of a Lennard-Jones fluid is much longer than the interaction length of a hard-sphere fluid and can become of the same order as the pore size. This can result in the breakdown of Hadwiger’s theorem. Minkowski functionals can be used in combination with experiments, theory, or simulations. In this work, classical density functional theory (cDFT) is employed [33] to compute the free energy and adsorption of a Lennard-Jones fluid for a wide variety of geometries. This simulation technique is commonly used to study capillary condensation and wetting [34,35,36,37] and can capture larger length scales than molecular dynamics simulations. The Minkowski functionals are fitted to the cDFT results to study how well they can describe the complex phase behavior of a Lennard-Jones fluid under tight confinement.

We find that Minkowski functionals can provide a good description of Lennard-Jones fluids down to fairly small system sizes. In addition, it is found that topology has a different effect on the phase envelope shift of a Lennard-Jones fluid under confinement depending on whether the phase transition under consideration is a continuous or discrete (first-order) transition. The reason behind this intriguing behavior is still under active investigation.

2. Theory

2.1. Minkowski Functionals

The Minkowski functionals are a concept from integral geometry. These functionals characterize both the morphology and the topology of a spatial pattern [38,39]. For a D dimensional space, there are $D+1$ functionals. Considering a 2D system with a surface, X, and a smooth boundary, $\delta X$, the following functionals can be defined:

$$\begin{array}{cc}{\displaystyle M_0\left(X\right)={\int }_{\delta X}dA}\hfill & =A\left(X\right)\hfill \\ {\displaystyle M_1\left(X\right)=\frac{1}{2}{\int }_{\delta X}dL}\hfill & =C\left(X\right)\hfill \\ {\displaystyle M_2\left(X\right)=\frac{1}{2}{\int }_{\delta X}k\left(X\right)dL}\hfill & =K\left(X\right)\hfill \end{array}$$

(1)

where $A\left(X\right)$ is the total surface area, $dA$, a surface element, $C\left(X\right)$, the total circumference, $dL$ a circumference element, $K\left(X\right)$, the total signed curvature, and $k\left(X\right)$ is the local signed curvature [40]. Following the Gauss–Bonnet theorem, in the absence of geodesic curvature, the signed curvature is directly proportional to the Euler characteristic, $\chi$, which is a measure of connectivity/topology. Now consider a functional, $\mathcal{M}\left(X\right)$, which is additive:

$$\mathcal{M}\left(X_1\cup X_2\right)=\mathcal{M}\left(X_1\right)+\mathcal{M}\left(X_2\right)-\mathcal{M}\left(X_1\cap X_2\right),$$

(2)

motion invariant:

$$\mathcal{M}\left(gX\right)=\mathcal{M}\left(X\right),$$

(3)

and continuous:

$$\mathcal{M}\left(X_n\right)\to \mathcal{M}\left(X\right)\mathit{if}:X_n\to X\mathit{for}:n\to \infty .$$

(4)

Then, following Hadwiger’s theorem [30], this functional, $\mathcal{M}\left(X\right)$, can be expressed as a linear combination of Minkowski functionals, $M_{\nu }\left(X\right)$, as follows:

$$\mathcal{M}\left(X\right)=\sum _{\nu =0}^d{c}_{\nu }{M}_{\nu }\left(X\right).$$

(5)

Additivity implies that the functional of the domain consisting of both $X_1$ and $X_2$ is equal to the sum of the functionals of the individual domains $X_1$ and $X_2$, minus the boundary. Motion invariance signifies that the value of the functional is independent of the translation or rotation of the domain X. Lastly, continuity means that as a series of domains $X_n$ converges to domain X, the functional of the series $X_n$ should converge to the functional of X as well. An example of such a functional is the grand potential or Landau free energy, $\Omega \left(\mu ,V,T\right)$, which is commonly used to study capillary condensation [41,42]. For a 2D system, the grand potential can thus be written as [20]:

$$\frac{\Omega \left(X\right)}{L}=-p\left(\mu ,T\right)A\left(X\right)+\sigma \left(\mu ,T\right)C\left(X\right)+\kappa \left(\mu ,T\right)K\left(X\right),$$

(6)

where L is a unit length perpendicular to the 2D system, $p\left(\mu ,T\right)$ is the average local pressure, $\sigma \left(\mu ,T\right)$ is the surface tension between the liquid and the wall, and $\kappa \left(\mu ,T\right)$ is the signed bending rigidity. The above expression for the grand potential demonstrates the importance of Hadwiger’s theorem. The pressure, surface tension, and signed bending rigidity are all only dependent on the chemical potential, $\mu$, and the temperature, T. Thus, the above equation shows how thermodynamics can be separated from morphology and topology [31]. Once the grand potential is known, other thermodynamic properties can be derived. This includes the excess free energy or surface tension [13,43]:

$$\begin{array}{cc}\hfill \gamma =& \frac{\left(\Omega /L+p_{b}A\right)}{C}\hfill \\ \hfill =& \left(p_b\left(\mu ,T\right)-p\left(\mu ,T\right)\right)\frac{A\left(X\right)}{C\left(X\right)}\hfill \\ \hfill & +\sigma \left(\mu ,T\right)+\kappa \left(\mu ,T\right)\frac{K\left(X\right)}{C\left(X\right)},\hfill \end{array}$$

(7)

and, through Gibbs’ theorem, the excess adsorption:

$$\begin{array}{cc}\hfill -{\Gamma }_{\mathrm{ex}}=& \frac{1}{C}\int {\rho }_b-\rho \left(r\right)dA={\left(\frac{\partial \gamma }{\partial \mu }\right)}_{T,V}\hfill \\ \hfill =& \frac{\partial }{\partial \mu }\left(p_b-p\left(\mu ,T\right)\right)\frac{A\left(X\right)}{C\left(X\right)}\hfill \\ \hfill & +\frac{\partial \sigma }{\partial \mu }+\frac{\partial \kappa }{\partial \mu }\frac{K\left(X\right)}{C\left(X\right)}.\hfill \end{array}$$

(8)

In the above equations, $p_b$ is the bulk pressure and ${\rho }_b$ is the bulk density. While it has been suggested in the past that Gibbs’ theorem is not valid for some systems undergoing capillary condensation [34,41], later work found that Gibbs’ theorem is not violated when using an arc length tracking algorithm [13,44].

In addition to the excess adsorption, one can also compute the effect of confinement on the phase envelope [43]:

$$\Delta p\left(\mu ,T\right)={\sigma }_{lg}^{\prime }\left(\mu ,T\right)\frac{C\left(X\right)}{A\left(X\right)}+{\kappa }_{lg}^{\prime }\left(\mu ,T\right)\frac{K\left(X\right)}{A\left(X\right)},$$

(9)

With:

$$\begin{array}{cc}\hfill {\sigma }_{lg}^{\prime }\left(\mu ,T\right)=& {\sigma }_{sg}\left(\mu ,T\right)-{\sigma }_{sl}\left(\mu ,T\right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\kappa }_{lg}^{\prime }\left(\mu ,T\right)=& {\kappa }_{sg}\left(\mu ,T\right)-{\kappa }_{sl}\left(\mu ,T\right),\hfill \end{array}$$

(11)

where ${\sigma }_{sg}\left(\mu ,T\right)$ and ${\sigma }_{sl}\left(\mu ,T\right)$ are the solid-gas and solid-liquid surface energies, respectively, and ${\kappa }_{sg}\left(\mu ,T\right)$ and ${\kappa }_{sl}\left(\mu ,T\right)$ are the solid-gas and solid-liquid bending rigidities. The above equation is a generalization of the Kelvin equation, and an equation of the same form can be derived for the temperature shift [45]. Two points should be noted about the above equation: (i) due to diverging density fluctuations, a mean field approach is not expected to fully capture the correct scaling at the critical point [46], and (ii) the correlation length that measures the range of density fluctuations at the critical point also diverges, resulting in a potential violation of Hadwiger’s theorem. When the characteristic interaction length scale becomes of the same order as the system size it is no longer additive. However, while these points need to be investigated further, the above equation should give a good first approximation of how phase behavior is affected by topology.

3. Methods

3.1. Density Functional Theory

The Minkowski functionals can be used with either experiments, theory, or simulations. In this work, classical density functional theory (cDFT) is used to compute the coefficients in front of the Minkowski functionals. DFT is a mean field approach that was first developed for quantum mechanics [47] but was later adapted to describe classical mechanical systems [33] as well. This mean field approach has the advantage of giving a description of the physics at the nanoscopic molecular level while scaling up to the mesoscopic level at which capillary condensation occurs.

The two basic assumptions of density functional theory are (i) the Hohenberg-Kohn variational principle, which states that there is a functional of the ground state free energy which can be fully recovered from the ground-state one-particle density distribution, and (ii) the Gibbs’ inequality, which states that any particle density distribution that is not the ground state will have a higher free energy than the ground state free energy [47]. Formulated in the grand canonical ($\mu$,V,T) ensemble, at the most basic level this means that classical DFT solves the following minimization problem:

$$\frac{\delta \Omega }{\delta \rho \left(\mathit{r}\right)}=0,$$

(12)

where $\Omega$ is the grand potential or Landau free energy, $\rho$ is the density, and $\delta$ is the Fréchet (functional) derivative [48]. To solve the above equation, we use the cDFT solver Tramonto, which is developed at the Sandia National Laboratories [13,44,48,49,50]. This code uses perturbation theory where the grand potential is split up as:

$$\Omega =F_{\mathrm{id}}+F_{\mathrm{hs}}+F_{\mathrm{p}}-\int d\mathit{r}\rho \left(\mathit{r}\right)\left[V\left(\mathit{r}\right)-\mu \right],$$

(13)

where $F_{\mathrm{id}}$ is the ideal contribution, $F_{\mathrm{hs}}$ is the hard-sphere contribution, and $F_{\mathrm{p}}$ is the perturbation contribution. $V\left(\mathit{r}\right)$ is the external potential resulting from (pore) walls acting on the fluid. The individual contributions are given by the following integrals:

$$\begin{array}{cc}\hfill F_{\mathrm{id}}& ={\beta }^{-1}\int d\mathit{r}\rho \left(\mathit{r}\right)\left\{ln\left[{\Lambda }^3\rho \left(\mathit{r}\right)\right]-1\right\}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill F_{\mathrm{hs}}& =\int d\mathit{r}\Phi \left\{{\overline{\rho }}_{\gamma }\left(\mathit{r}\right)\right\}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill F_{\mathrm{p}}& =\frac{1}{2}\int d\mathit{r}\int {\mathit{r}}^{\prime }\rho \left(\mathit{r}\right)\rho \left({\mathit{r}}^{\prime }\right)U_p\left(|\mathit{r}-{\mathit{r}}^{\prime }|\right).\hfill \end{array}$$

(16)

In the above equations, ${\beta }^{-1}=k_{B}T$, with $k_B$ the Boltzmann constant, and T the temperature, $\Lambda$ is the thermal de Broglie wavelength, $\Phi$ is the excess free energy density which is a function of ${\overline{\rho }}_{\gamma }$, a set of weighted non-local densities, and $U_p$ is an interaction potential. This potential is based on the Weeks-Chandler-Anderson approach [51], which splits an interaction potential as $U_p\left(r\right)=u\left(r_{\min }\right)$ for $r<r_{\min }$ and $U_p\left(r\right)=u\left(r\right)$ for $r\ge r_{\min }$. The potential $u\left(r\right)$ is a cut and shifted Lennard-Jones potential with $u\left(r\right)=u_{\mathrm{LJ}}\left(r\right)-u_{\mathrm{LJ}}\left(r_c\right)$, where:

$$u_{\mathrm{LJ}}\left(r\right)=4{ϵ}_{\mathrm{ff}}\left[{\left(\frac{{\sigma }_{\mathrm{ff}}}{r}\right)}^{12}-{\left(\frac{{\sigma }_{\mathrm{ff}}}{r}\right)}^6\right],$$

(17)

and $r_c={\sigma }_{\mathrm{ff}}$. Here ${ϵ}_{\mathrm{ff}}$ is the depth of the potential well and ${\sigma }_{\mathrm{ff}}$ is the finite distance at which the potential is zero. In this work, the Fundamental Measure Theory (FMT) is used with the White Bear functional [52]. The weighted non-local densities are:

$${\overline{\rho }}_{\gamma }\left(\mathit{r}\right)=\int d{\mathit{r}}^{\prime }\rho \left({\mathit{r}}^{\prime }\right)w^{\left(\gamma \right)}\left(|\mathit{r}-{\mathit{r}}^{\prime }|\right),$$

(18)

with the weight functions:

$$\begin{array}{cc}\hfill w^{\left(3\right)}\left(r\right)& =\theta \left(r-R\right),\hfill \end{array}$$

(19)

$$\begin{array}{cccccc}\hfill w^{\left(2\right)}\left(r\right)& =4\pi R{w}^{\left(1\right)}\left(r\right)\hfill & & =4\pi R^2{w}^{\left(0\right)}\left(r\right)\hfill & & \hfill \\ & & & =\delta \left(r-R\right),\mathrm{and}\hfill & & \end{array}$$

(20)

$$\begin{array}{ccccc}\hfill w^{\left(V2\right)}\left(r\right)& =4\pi R{w}^{\left(V1\right)}\left(r\right)\hfill & & =\left(1/r\right)\delta \left(r-R\right).\hfill & \end{array}$$

(21)

The excess free energy density is given by $\Phi ={\Phi }_s+{\Phi }_v$ with:

$$\begin{array}{cc}\hfill {\Phi }_s& =-{\overline{\rho }}_{0}ln\left(1-{\overline{\rho }}_3\right)+\frac{{\overline{\rho }}_1{\overline{\rho }}_2}{1-{\overline{\rho }}_3},\mathrm{and}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\Phi }_v& =-\frac{{\overline{\rho }}_{V1}·{\overline{\rho }}_{V2}}{1-{\overline{\rho }}_3}+\frac{1}{36\pi {\overline{\rho }}_3^2{(1-{\overline{\rho }}_3)}^2}\times \hfill \\ \hfill & {\left({\overline{\rho }}_2-\frac{{\overline{\rho }}_{V2}·{\overline{\rho }}_{V2}}{{\overline{\rho }}_2}\right)}^3\left({\overline{\rho }}_3+{(1-{\overline{\rho }}_3)}^{2}ln\left(1-{\overline{\rho }}_3\right)\right).\hfill \end{array}$$

(23)

The last term that needs to be defined is the external potential, $V\left(r\right)$, which is defined as:

$$V\left(\mathit{r}\right)={\rho }_s\int d{\mathit{r}}_s{v}_{\mathrm{LJ}}\left(|\mathit{r}-{\mathit{r}}_s|\right)-v_{\mathrm{LJ}}\left(r_c\right),$$

(24)

where the integral is taken over all elements assigned to the (pore) wall. The potential $v_{\mathrm{LJ}}\left(r\right)$ is the same as the Lennard-Jones potential defined in Equation (17), but with ${ϵ}_{\mathrm{ff}}$ replaced with ${ϵ}_{\mathrm{sf}}$ and ${\sigma }_{\mathrm{ff}}$ replaced with ${\sigma }_{\mathrm{sf}}$. More details about the discretization of the above equations, their numerical implementation, and how to solve them in parallel can be found in the literature [13,44,48,49,50]. Phase transitions are tracked using the pseudo arc length continuation algorithm of Keller [44,53] which has been implemented in the LOCA software library [54].

3.2. Geometries

As mentioned in the introduction, the Minkowski functionals framework allows for the separation of thermodynamics and geometry. Therefore, one can study a problem using simple geometries and extrapolate the results to complex geometries with known Minkowski functionals. Figure 1 shows these various 2D geometries and topologies. Along the vertical axis, the various shapes show cross-sections of pores with different radii, $r_p$. Along the horizontal axis rods with radius $r_r=1.0$ (≈$0.36\mathrm{nm}$) are placed inside the pores to modify the Euler characteristic $\chi$. While these specific pores may or may not be found in nature or technology, they represent any material with a similar topology. Materials with a negative Euler characteristic are common in nature [55,56]. A negative Euler characteristic indicates that the pore space is multiple times connected. The various 2D Minkowski functionals associated with the surface area, circumference, and signed curvature, respectively, can be computed with the following set of equations:

$$\begin{array}{cccc}\hfill A& =\pi r_p^2\hfill & & -n_r\pi r_r^2,\hfill \end{array}$$

(25)

$$\begin{array}{cccc}\hfill C& =\pi r_p\hfill & & +n_r\pi r_r,\hfill \end{array}$$

(26)

$$\begin{array}{cccc}\hfill K& =\pi \hfill & & -n_r\pi \phantom{r_r},\hfill \end{array}$$

(27)

where $n_r$ is the number of rods inside the pore. As mentioned in Section 2.1, for the Minkowski functionals to accurately capture the physics of capillary condensation, the conditions in Hadwiger’s theorem need to be met. Considering that for very small pores the characteristic interaction length between molecules becomes of the same order as the pore size, the additivity constraint (Equation (2)) is expected to break down first. In the literature it is reported that an error of about $1\%$ is found when the system size becomes of the order of ≈$10{\sigma }_{\mathrm{ff}}$, where ${\sigma }_{\mathrm{ff}}$ is the characteristic length scale of the interaction potential between molecules [31]. For the pores in gray in Figure 1 the minimum distance between the walls of the pore and/or the rods inside the pore is smaller than this distance, and Hadwiger’s theorem is expected to break down. Another potential source of error is that Hadwiger’s theorem applies to convex geometries while the pores are concave.

3.3. Simulation Parameters

Existing literature has focused on the behavior of hard-sphere fluids [31] and TraPPE force fields [32]. In this work, classical density functional theory (cDFT) simulations have been performed for a Lennard-Jones fluid. The combination of a Lennard-Jones potential with cDFT allows for more complex behavior than a hard-sphere potential, while also reducing the need for computational resources compared to the TraPPE force fields combined with molecular dynamics simulations.

Table 1. cDFT parameters of N₂ and SiO₂ [60]. The number density for SiO₂ is ${\rho }_{\mathrm{s}}=66.15{\mathrm{nm}}^{-3}$ [35]. Fluid-fluid interactions are truncated at $5{\sigma }_{\mathrm{ff}}$. The simulations are performed at $77.3\mathrm{K}$.

	${\mathit{ϵ}}_{\mathbf{ff}}/{\mathit{k}}_{\mathit{B}}$	${\mathit{\sigma }}_{\mathbf{ff}}$	${\mathit{d}}_{\mathbf{HS}}$	${\mathit{ϵ}}_{\mathbf{sf}}/{\mathit{k}}_{\mathit{B}}$	${\mathit{\sigma }}_{\mathbf{sf}}$
	$\left[\mathbf{K}\right]$	$\left[\mathbf{nm}\right]$	$\left[\mathbf{nm}\right]$	$\left[\mathbf{K}\right]$	$\left[\mathbf{nm}\right]$
N2	94.45	0.3575	0.3575	147.3	0.317

shows the parameters used in the cDFT simulations. Because it is a commonly used model system [57,58,59], the Lennard-Jones fluid is parameterized as Nitrogen in Vycor glass. The parameters are the same as those used by Ravikovitch et al. [60] and Ustinov et al. [35] and are very similar to the parameters used by Gelb & Gubbins [61] in their grand canonical Monte Carlo simulations of nitrogen in Vycor glass. Figure 2 shows a comparison between cDFT simulations and experiments [62] for the adsorption isotherm of N₂ in a SiO₂ slit pore with a width of $L/{\sigma }_{\mathrm{ff}}=500$. The results confirm that the parameters listed in Table 1 are a reasonable choice. With more advanced models for the interaction between N₂ and SiO₂, a better match can be obtained between cDFT simulations and experiments [35,60]. However, the choice of the same potential for particle-particle and wall-particle keeps the system simple and the results easier to interpret. The computations are performed in the grand canonical ($\mu ,V,T$) ensemble and the relation between the chemical potential and pressure is obtained from a bulk cDFT simulation. Since the pore space is in equilibrium with the bulk, both the temperature and chemical potential are constant throughout the system [63].

3.4. Minkowski Coefficients

To compute the Minkowski coefficients, $p\left(\mu ,T\right)$, the local pressure, $\sigma \left(\mu ,T\right)$, the surface tension between the liquid and the wall, and $\kappa \left(\mu ,T\right)$ the signed bending rigidity, for a given temperature, T, and chemical potential $\mu$, one needs to know the grand potential and Minkowski functionals for three different geometries. To improve accuracy, cDFT simulations are performed for all geometries shown in Figure 1 and one set of coefficients is computed by performing a least-squares fit to all the simulation data. Any simulation technique or theory that computes free energy can thus be used in combination with the Minkowski functionals. For experiments, the excess adsorption can be measured instead.

4. Results

In this section, the results for the Lennard-Jones fluid cDFT simulations are presented, and how these results can be described using the Minkowski functionals. All the parameters have been made dimensionless with $k_{\mathrm{b}}T$ and ${\sigma }_{\mathrm{ff}}$. An overview of the various dimensionless variables is shown in

Table 2. An overview of the dimensionless variables used in the results section. “$\widehat{}$” denotes the dimensional equivalent of a dimensionless variable. In addition, $k_b$ is the Boltzmann constant, ${\sigma }_{\mathrm{ff}}=0.3575\mathrm{nm}$, $T=77.3\mathrm{K}$, and L is one dimensionless unit length.

Description	Symbol	Dimensionless Definition
Distance	$\mathrm{length}$	$\widehat{\mathrm{length}}/{\sigma }_{\mathrm{ff}}$
Density	$\rho$	$\widehat{\rho }{\sigma }_{\mathrm{ff}}^3$
Grand potential	$\Omega /L$	$\widehat{\Omega }{\sigma }_{\mathrm{ff}}/\widehat{L}{k}_{b}T$
Excess adsorption	$\Gamma /L$	$\widehat{\Gamma }{\sigma }_{\mathrm{ff}}/\widehat{L}$
Chemical potential	$\mu$	$\widehat{\mu }/k_{b}T$
Pressure	p	$\widehat{p}{\sigma }_{\mathrm{ff}}^3/k_{b}T$
Surface tension	$\sigma$	$\widehat{\sigma }{\sigma }_{\mathrm{ff}}^2/k_{b}T$
Bending rigidity	$\kappa$	$\widehat{\kappa }{\sigma }_{\mathrm{ff}}/k_{b}T$
Pseudo pressure	$p^{\prime }$	${\widehat{p}}^{\prime }{\sigma }_{\mathrm{ff}}^2/k_{b}T$
Pseudo surface tension	${\sigma }^{\prime }$	${\widehat{\sigma }}^{\prime }{\sigma }_{\mathrm{ff}}^{7/4}/k_{b}T$

.

Figure 3 shows some typical results of our cDFT simulations; it shows the 2D dimensionless free energy $\Omega /L$ as a function of the dimensionless chemical potential, $\mu$, for $\chi =1$. For a pore with radius $r_p=25$ (≈$9\mathrm{nm}$) three different regimes can be identified. Initially, when the chemical potential is very low, the pore is completely empty. At about $\mu \approx -10$ (≈$-0.067\mathrm{eV}$) a phase transition can be observed, which marks the onset of gas adsorption on the pore wall. Eventually, at about $\mu \approx -6$ (≈$-0.040\mathrm{eV}$) another phase transition can be observed, which is the capillary condensation phase transition. Consistently with the literature, the capillary condensation phase transition in this 2D system is a first-order phase transition [46]. It is marked by a van der Waals loop in the grand potential and a jump in the excess adsorption as function of the chemical potential [45].

The three different regimes are also identified in the literature [13,46]. The two phase transitions are shown in more detail in Figure 3a–d. In Figure 3a the derivative of the free energy changes continuously, suggesting that the adsorption of gas on the pore wall is a continuous phase transition. Figure 3d shows the matching divergence of the derivative of the excess adsorption. In Figure 3b the derivative of the free energy changes abruptly, indicating a discrete (first-order) phase transition. Figure 3d shows the matching jump in the excess adsorption. Consistently with the literature, as the pore radius becomes smaller the phase behavior changes [64]. For pores smaller than $r_p=7.5$ (≈$2.7\mathrm{nm}$) the discrete phase transition can no longer be observed and becomes continuous instead.

The continuous and discrete phase transitions are now further analyzed using Equation (9), which describes the phase envelope shift as a function of the Minkowski functionals. By only selecting geometries with $K\left(X\right)=0$, one can evaluate the coefficient ${\sigma }_{lg}^{\prime }$. Figure 4 shows the dimensionless pressure difference between the continuous phase transition associated with adsorption of gas onto a wall in the bulk and inside a pore times the surface area Minkowski functional, $\Delta p_w·A\left(X\right)$, for $K\left(X\right)=0$ [46]. This pressure difference can be predicted well by the function: $f_w\left(C^{3/4}\right)={\Omega }_{\lg }+{\sigma }_{\lg }^{\prime }{C}^{3/4}$, where the coefficients are: ${\Omega }_{\lg }=-8.3\times {10}^{-3}\pm 0.1\times {10}^{-3}$ and ${\sigma }_{\lg }^{\prime }=1.538\times {10}^{-3}\pm 0.008\times {10}^{-3}$. The discrepancy between the power of $3/4$ and the theoretical prediction of 1 shows how, due to the small system size, Hadwiger’s additivity assumption starts to break down. Excluded volume effects keeping molecules away from the wall could be another potential source of error.

To find the value of ${\kappa }_{\lg }^{\prime }$ and to validate how well the function $f_w\left(C^{3/4}\right)$ describes the data, Figure 5 shows the pressure difference, $\Delta p_w$, times the surface area Minkowski functional minus the function $f_w\left(C^{3/4}\right)$, $\Delta p_w·A\left(X\right)-f_w\left(C^{3/4}\right)$, as a function of the Minkowski functional, $K\left(X\right)$. The graph shows a collapse of the data and, as expected, a linear fit with: $g_w\left(K\right)={\kappa }_{\lg }^{\prime }K\left(X\right)$, where ${\kappa }_{\lg }^{\prime }=4.6\times {10}^{-4}\pm 0.2\times {10}^{-4}$. The collapse of the data confirms that the grand potential and thus the continuous pressure shift in the phase envelope times the surface area Minkowski functional, $\Delta p_w·A\left(X\right)$, is proportional to $C^{3/4}\left(X\right)$ and is linearly dependent on the topology of the system $K\left(X\right)$. In addition, the data shows that the sensitivity of the pressure shift, $\Delta p_w$, to topology is about an order of magnitude smaller than the sensitivity to the pseudo surface tension, ${\sigma }_{\lg }^{\prime }$. The outlier, $r_p=2.5$ (≈$0.9\mathrm{nm}$), $\chi =1$, confirms the breakdown of Hadwiger’s theorem for small pore sizes.

A similar analysis can be performed for the phase envelope shift of the discrete capillary condensation phase transition. Figure 6 shows the pressure difference between the capillary condensation pressure and the bulk phase transition pressure times the surface area Minkowski functional, $\Delta p·A\left(X\right)$, as a function of pore radius, $r_p$. Because the rod size is constant across different simulations, the Minkowski functionals $A\left(X\right)$ and $C\left(X\right)$ can be expressed as a function of only the pore radius $r_p$ for constant bending rigidity, $K\left(X\right)$. The Minkowski functional for topology is zero in all the cDFT simulations shown in Figure 6, $K\left(X\right)=0$. The curve fit is equal to: $f(C^{3/4},A^{1/2})={\Omega }_{\lg }+{\sigma }_{\lg }^{\prime }{C}^{3/4}+p_{\lg }^{\prime }{A}^{1/2}$, where the Minkowski functionals $C\left(X\right)$ and $A\left(X\right)$ only depend on $r_p$. The coefficients are: ${\Omega }_{\lg }=-0.4\pm 0.2$, ${\sigma }_{\lg }^{\prime }=0.12\pm 0.03$, and $p_{\lg }^{\prime }=-0.04\pm 0.02$.

To find the value of ${\kappa }_{\lg }^{\prime }$ and to validate the data fit shown in Figure 6, Figure 7 shows the pressure difference, $\Delta p$, times the surface area Minkowski functional minus the function from Figure 6, $\Delta p·A\left(X\right)-f(C^{3/4},A^{1/2})$, as a function of the Minkowski functional, $K\left(X\right)$. The graph shows a collapse of the data and, as predicted by theory, a linear fit with: $g_w\left(K\right)={\kappa }_{\lg }^{\prime }K\left(X\right)$, where ${\kappa }_{\lg }^{\prime }=4.6\times {10}^{-4}\pm 0.2\times {10}^{-4}$. The collapse of the data into one single line, confirms that the pressure difference times the surface area Minkowski functional, $\Delta p·A\left(X\right)$, is proportional to $C^{3/4}\left(X\right)$, $A^{1/2}\left(X\right)$, and $K\left(X\right)$. The fitting parameters ${\sigma }_{\lg }^{\prime }$ and $p_{\lg }^{\prime }$ show that the difference between the capillary condensation pressure and the bulk phase transition pressure is the most sensitive to changes in the pseudo surface tension, ${\sigma }_{\lg }^{\prime }$, and the pseudo pressure, $p_{\lg }^{\prime }$. The sensitivity to topology changes, represented by the parameter ${\kappa }_{\lg }^{\prime }$, is much smaller. The outlier, $r_p=10.0$ (≈$3.6\mathrm{nm}$), $\chi =-3$, is a system where the topology change compared to, $\chi =1$, caused an additional phase transition to occur.

Surprisingly, Figure 5 and Figure 7 show that the phase envelope shift dependence on topology is different depending on whether the phase transition is a discrete first-order or a continuous phase transition. In the case of gas adsorbing onto the pore wall, adding rods to the system decreases the pressure shift in the phase envelope. In the case of the discrete capillary condensation phase transition adding rods to a pore increases the pressure shift. What is behind this behavior and what happens when the discrete capillary condensation phase transition turns into a continuous phase transition at the critical pore radius is a topic for further research. In addition, the physical meaning of the pseudo pressure and the pseudo surface tension needs to be investigated further.

Like Figure 3, Figure 8 shows the dimensionless 2D grand potential, $\Omega /L$, as a function of the dimensionless chemical potential, $\mu$. Only the simulation results with Euler characteristic, $\chi =1$, are shown. The different lines show the results of the cDFT simulations while the symbols show the grand potential reconstructed from the Minkowski functionals and one set of Minkowski functional coefficients: pressure, $p\left(\mu ,T\right)$, surface tension, $\sigma \left(\mu ,T\right)$, bending rigidity, $\kappa \left(\mu ,T\right)$, and the pseudo pressure and surface tension terms $p^{\prime }\left(\mu ,T\right)$ and ${\sigma }^{\prime }\left(\mu ,T\right)$. The terms $p^{\prime }\left(\mu ,T\right)A^{1/2}\left(X\right)$ and ${\sigma }^{\prime }\left(\mu ,T\right)C\left(X\right)$ are added to the expansion of the grand potential based on the analysis of the phase behavior above. Both these terms are in effect circumference (surface) terms. While for a pore with Euler characteristic $\chi =1$, $A^{1/2}\left(X\right)\propto C\left(X\right)$, these two contributions diverge for $\chi \ne 1$. Adding the term $p^{\prime }\left(\mu ,T\right)A^{1/2}\left(X\right)$ to the expression for the grand potential significantly reduced the error in the Minkowski functional fit compared to only having a contribution from the circumference, $C\left(X\right)$.

The corresponding dimensionless Minkowski functional coefficients as a function of the dimensionless chemical potential, $\mu$, are shown in Figure 9a. It can be observed that the pressure coefficient is very similar to the bulk pressure, $p_b$. All coefficients show a large peak around $\mu \approx -6$ (≈$-0.040\mathrm{eV}$) and the range of these peaks extends from about $-150$ to 150. Before these peaks occur, the surface tension and bending rigidity terms show very similar behavior. To get a better understanding of the behavior of a Lennard-Jones fluid under confinement, one can analyze how much individual Minkowski functionals contribute to the grand potential. In Figure 9b, the value of the Minkowski functional coefficients times their corresponding Minkowski functionals is shown for a pore with radius, $r_p=25$ (≈$9\mathrm{nm}$), and an Euler characteristic of, $\chi =-5$. This analysis suggests that the adsorption of gas onto the pore wall is dominated by the surface tension and the pseudo pressure. When capillary condensation occurs, the surface tension and pseudo pressure contributions both become discontinuous and show large increases. Also, the contribution from the pseudo surface tension contribution becomes significant. After the bulk phase transition, the system becomes increasingly dominated by the pressure. The contribution from topology is the third largest, after the surface tension and the pseudo pressure.

Figure 10 shows the average absolute relative error, $<|{ϵ}_{\Omega }|>$, as a function of the square root of the minimal characteristic length scale of the system, $\sqrt{l_{\min }}$. The average absolute relative error is defined as:

$$〈|{ϵ}_{\Omega }|〉={〈|\frac{{\Omega }_{\mathrm{cDFT}}-{\Omega }_{\mathrm{Mink}}}{{\Omega }_{\mathrm{cDFT}}}|〉}_{\mu },$$

(28)

where ${\Omega }_{\mathrm{cDFT}}$ is the grand potential computed using cDFT, and ${\Omega }_{\mathrm{Mink}}$ is the reconstruction of the grand potential using the Minkowski functionals. The error is averaged with respect to the chemical potential, $\mu$. Due to the longer interaction length of the Lennard-Jones potential compared to the hard-sphere potential and the more complex phase behavior, the observed error is larger than the error reported in the literature where an error of $1\%$ is found when the system size becomes of the order of ≈$10{\sigma }_{\mathrm{ff}}$ [31]. While the cut-off length of the interaction potential is equal to $5\sigma$ in the cDFT simulations the Lennard-Jones potential is almost zero at $2\sigma$. One could therefore expect the error to be significantly less than $1\%$ at a system size of ≈$50{\sigma }_{\mathrm{ff}}$, which is indeed the case. This might be in part due to the addition of the pseudo pressure and surface tension terms to the Minkowski functional expression for the grand potential. For pores without rods, it can be observed that $log(<|{ϵ}_{\Omega }|>)$ scales almost linearly with $\sqrt{l_{\min }}$. However, this scaling does not hold for pores with rods. Changing the topology of the system with rods does not have a large effect on the error, $<|{ϵ}_{\Omega }|>$, especially for smaller pores. The fact that, as shown above, the topology has a smaller contribution to the grand potential than the various surface contributions could partly explain this observation. The figure confirms the observation in Figure 8 that for larger pores the fit of the Minkowski functional reconstruction of the grand potential is much better than for smaller pores.

The dimensionless 2D excess adsorption, $\Gamma /L$, as a function of the dimensionless chemical potential, $\mu$, for a Lennard-Jones fluid is shown in Figure 11. Again, only the simulation results when $\chi =1$ are shown. The different lines show the results of the cDFT simulations while the symbols show the grand potential as reconstructed from the Minkowski functionals and one set of Minkowski functional coefficients: the derivatives of pressure, $\partial (p-p_b)/\partial \mu$, surface tension, $\partial \sigma /\partial \mu$, bending rigidity, $\partial \kappa /\partial \mu$, and the pseudo pressure and surface tension terms $\partial p^{\prime }/\partial \mu$ and $\partial {\sigma }^{\prime }/\partial \mu$ with respect to the chemical potential. As is the case for the grand potential in Figure 8, different regimes can be identified in Figure 11. For the pore size, $r_p=25$ (≈$9\mathrm{nm}$), the following regimes can be observed: (i) at low chemical potential the pores are completely empty, (ii) starting at about $\mu \approx -10$ (≈$-0.067\mathrm{eV}$), a continuous phase transition can be observed and gas starts adsorbing on the wall [46], (iii) at about $\mu \approx -6$ capillary condensation can be observed, and iv) around $\mu \approx -5$ (≈$-0.033\mathrm{eV}$) the bulk phase transition occurs and the excess adsorption shows a significant drop. The inset shows that for larger pores there is a very good match between the cDFT computations and the Minkowski functional reconstructions. However, for the smallest pores, the match is quite poor and neither the film adsorption stage nor capillary condensation is captured well.

The derivatives of the Minkowski functional coefficient with respect to the chemical potential, $\mu$, used to reconstruct the excess adsorption, can be observed in Figure 12a. This plot confirms that the pressure term is very similar to the bulk pressure. The peaks at about $\mu \approx -10$ (≈$-0.067\mathrm{eV}$) and $\mu \approx -7$ (≈$-0.047\mathrm{eV}$) are the locations of the continuous phase transitions associated with layers of gas molecules adsorbing onto the pore wall. To be able to analyze how the different terms contribute to the excess adsorption, in Figure 12b the value of the Minkowski functional coefficients times their corresponding Minkowski functionals is shown for a pore with radius $r_p=25$ (≈$9\mathrm{nm}$) and Euler characteristic of $\chi =-5$. Since these coefficients are derivatives of the coefficients used to reconstruct the grand potential, they can be both positive and negative. The amount of gas adsorbed onto the wall is dominated by the pseudo pressure and a negative contribution from the surface tension. The pseudo surface tension only contributes during the continuous phase transitions at about $\mu \approx -10$ (≈$-0.067\mathrm{eV}$) and $\mu \approx -7$ (≈$-0.047\mathrm{eV}$). Capillary condensation is characterized by a much larger contribution of the pseudo surface tension and many discontinuities in all the different terms to accommodate the discontinuity of a first-order phase transition. At this value of the Euler characteristic, topology does not seem to contribute significantly to the excess adsorption.

Figure 13 shows the average absolute relative error $<|{ϵ}_{\Gamma }|>$ as a function of the minimal characteristic length scale of the system, $l_{\min }$. Because excess adsorption is a derivative of the grand potential, the observed error is larger than in Figure 10. For pores without rods, it can be observed that $log(<|{ϵ}_{\Gamma }|>)$ scales almost linearly with $l_{\min }$. Again, this scaling does not hold for pores with rods. This is similar to what was observed in Figure 10 which is probably at least partly caused by the limited contribution of topology to the excess adsorption.

5. Discussion & Conclusions

The effect of morphology and topology on capillary condensation is studied in a systematic manner using a Minkowski functional framework. Compared to a hard-sphere fluid, a Lennard-Jones fluid has a much longer interaction length and phase behavior becomes much more complex. This results in larger errors in both the grand potential and the excess adsorption, which can be partly overcome by introducing additional terms in the expansion of the grand potential. Analyzing the contributions of the various Minkowski functional coefficients to the grand potential shows that wall adsorption is dominated by a surface-tension and a pseudo-pressure term, which is proportional to $A^{1/2}$. During capillary condensation, a pseudo surface tension term proportional to $C^{3/4}$ also gains in importance and this regime is characterized by large discontinuities in the coefficients which match the discontinuities caused by capillary condensation. A similar analysis of the coefficients contributing to the excess adsorption shows the pseudo pressure as a positively contributing term to the wall adsorption and the surface tension as a negatively contributing term. During capillary condensation, the pseudo surface tension term becomes significantly more important. For both the grand potential and the excess adsorption the effect of topology is modest.

Last but not least, the effect of confinement on phase behavior is investigated. It is found that pressure shift in the continuous phase transition describing adsorption on the pore wall is proportional to $C^{3/4}\left(X\right)$ and $K\left(X\right)$, which is close to what is predicted based on Hadwiger’s theorem. However, the pressure for the discrete first-order phase transition describing capillary condensation is proportional to $C^{3/4}\left(X\right)$, $A^{1/2}\left(X\right)$ and $K\left(X\right)$, which is a deviation from theory. In addition, the discrete first-order phase transition and the continuous phase transition have the opposite dependence on topology. In the case of gas adsorbing onto the pore wall, adding rods to the system decreases the pressure shift in the phase envelope. In the case of the discrete capillary condensation phase transition adding rods to a pore increases the pressure shift. Whether this is a finding which holds in general for continuous and discrete phase transitions has to be investigated further.

The Minkowski functionals provide a useful framework to study capillary condensation. The separation of geometry and thermodynamics allows for a method to systematically study the effect of surface area, circumference, and the Euler characteristic on phase behavior. This provides many opportunities for future research. One of the many open questions is whether the different effect that topology has on continuous versus discrete phase transitions is also found for Lennard-Jones fluids in 3D instead of 2D systems, or even for phase transitions in general. Another question is whether it is possible to use the Minkowski functionals for upscaling. The usage of higher order terms breaks the additivity assumption. However, it might be possible to overcome this limitation by using the Minkowski functions, the Minkowski functionals as a function of a threshold value, e.g., the pore size [29]. The idea is to perform a number of simulations on different small geometries with known Minkowski functions and to compute the Minkowski function coefficients for the excess adsorption. These coefficients are then used to predict the excess adsorption for a much larger experimental disordered porous medium with known Minkowski functions. Other questions that remain to be answered are whether sorption hysteresis can be described using the Minkowski functionals [65], and to study whether the Minkowski functionals can be used for higher molecular weight molecules.