Advances in critical temperature shift modeling of confined pure fluids using the Kihara potential function

Abstract

A multitude of research works have been conducted in the past decade to better predict the change of critical properties of fluids confined in nanopores, known as critical shift, due to its great impact upon calculations of fluid properties in tight reservoirs. Modeling of this phenomenon commenced with developing equations of state (EOS) and has been continuing with correlations, all based on the two-parameter Lennard–Jones (L–J) potential function. Although these approaches have tried to present passable estimations of critical shift, sufficiently accurate predictions of critical shift are still missing in the literature. In this study, the three-parameter Kihara potential, as a more physically realistic alternative, is used to develop the van der Waal (vdW) EOS, and accordingly, a fluid-dependent expression is derived to calculate the critical temperature of confined fluids, i.e., pore critical temperature (${T}_{\mathrm{cp}}$). Using 50 data points of ${T}_{cp}$ reports for normal alkanes in the literature, the average error of our model is only 2.23%, 6.4% less than that of the L–J model. Furthermore, despite simple correlations of previous works, herein the Kihara parameters are exclusively tuned for each component based on their ${T}_{cp}$ reports, which resulted in an average error of 0.4% for normal alkanes. Finally, the pressure–volume diagrams of vdW and Peng–Robinson EOSs associated with the Kihara potential function are comprehensively discussed. The findings of this study can help researchers with more accurate predictions of the critical temperature of fluids confined in tight porous media, thereby providing more precise calculations of fluid properties and behavior at equilibrium conditions.

Introduction

In recent years, myriads of efforts have been made to understand and evaluate the phase behavior of fluids confined in molecular-scale media, especially hydrocarbons mainly because of the growing demand for production from unconventional reserves, such as shale oil and tight gas. The phenomenon of changes in the critical properties of matters due to confinement is originally known as critical shift and even its applicability has been investigated in several other research areas, such as nanofluidic chipset technologies (Alfi et al. 2016).

Phase behavior and thermodynamics of confined fluids have a significant deviation from their corresponding bulk behavior. Experimental measurements (Morishige and Ito 2002; Morishige and Shikimi 1998; Morishige et al. 2003; Tan et al. 2019) have shown that the thermodynamic and physicochemical properties of fluids, such as critical properties, substantially shift when they are confined in a porous medium with nanometric pore radius. Hence, incorporating the bulk properties (not confined state) of matters in phase behavior calculations undoubtedly leads to poor results.

Various analytical and numerical methods have been examined to predict and model the mentioned theory for pure fluids, e.g., the developed van der Waals (vdW) equation of state (EOS) in the previous study (Zarragoicoechea and Kuz 2002), extended Peng–Robinson and Soave–Redlich–Kwong EOSs (Islam and Sun 2017; Zhang et al. 2019a, b), density functional theory (Evans et al. 1986), and Monte Carlo simulations (Jiang et al. 2004; Jin and Nasrabadi 2016; Lowry and Piri 2018; Singh et al. 2009). Numerical approaches are frequently inconvenient to be applied to a different study scope since they are accompanied by numerous limitations. Having a theoretical and thermodynamic basis, the modified vdW EOS developed by Zarragoicoechea and Kuz has been widely used because of its compatibility with experimental values on the one hand, and its ease of use on the other. Such a convenience would be much preferable to implementing time-consuming and inconvenient huge numerical and molecular simulations for every research case. The EOS approach was developed based on the Lennard–Jones (L–J) potential function with two adjustable parameters ($\sigma_{{{\text{LJ}}}}$, $\varepsilon /k$) whose values can be determined via least-square data analysis of thermodynamic properties, such as experimental reports of viscosity or second virial coefficient (Tee et al. 1966a, b). A three-parameter potential function, nevertheless, is expected to provide more accurate modeling due to its flexibility and more realistic physics.

On the other hand, fluid adsorptions and their huge impacts on phase behavior could have not been ignored (Zhang et al. 2018a, b, c; Zhang et al. 2019a, b). A multi-layer film of fluid is always adsorbed to the pore walls and its contribution to the confined fluid characteristics is great enough to be considered, because in nanoscale (confined) porous media, the wall-molecule interactions are not insignificant anymore as they are in the bulk state. In this regard, analytical (Wang et al. 2021; Zhang et al. 2018a, b, c; Zhang et al. 2019a, b) and numerical (Jiang et al. 2005; Peng et al. 2007; Tan et al. 2017; Wong, et al. 2009) research works have been conducted to include the adsorption phenomenon. Although experimental measurements are always highly reliable and give direct understanding of fluid properties, they are greatly time-consuming and highly expensive setups are needed for confined porous media to be studied. Moreover, it is quite difficult to establish an apparatus that could last under extreme conditions of high pressure and temperature (Zhang et al. 2019a, b).

Over the last few years, various measures have been taken in order to enhance the confinement modeling of pure fluids, though no fundamental and analytical development was presented except the previous original model of L–J (Zarragoicoechea and Kuz 2002). Other fellow researchers have manipulated numbers and coefficients through curve fitting to merely modify every possible aspect of the basic model (Jin et al. 2013; Song et al. 2020; Yang et al. 2019; Yang and Li 2020). In several works, moreover, the mentioned L–J approach has just been repeated with corrections of the theory’s coefficients (Islam et al. 2015; Zhang et al. 2018a, b, c; Zhang et al. 2018a, b, c). These numerical enhancements are strongly related to the available data sets of one specific fluid, and possibly not appropriate for comprehensive applications. Such uncertainty even led to developing exclusively curve-fitted correlations for specific substances, e.g., methane, ethane and carbon dioxide (Tan et al. 2019; Yang et al. 2019).

Kihara potential, despite the Lennard–Jones, is a three-parameter potential function that incorporates the hard-impenetrable core of a molecule into the modeling by its third parameter ${a}_{k}$ (Kihara 1953). In the L–J model, however, there is no limitation for two approaching molecules and they can completely interpenetrate. Kihara successfully fits thermodynamic data, specifically the second virial coefficient, for various types of fluids. L–J potential, in contrast, performs weaker in lower temperatures, and the flexibility of Kihara potential makes it much easier to have a perfect fitting (Prausnitz et al. 1998).

In the current work, first, the Kihara potential function is employed to develop the energy integral and derive the vdW-EOS with which the hysteresis critical temperature ($T_{{{\text{ch}}}}$) of confined fluids including argon, nitrogen, oxygen, carbon dioxide, ethylene, and xenon is calculated. Second, the Kihara parameters are adjusted with a brand-new and concrete adjustment method on a thermodynamic basis to achieve compatibility with bulk critical properties. Accordingly, the pore critical temperature ($T_{{{\text{cp}}}}$) is calculated for some components comprising nitrogen, carbon dioxide, methane, ethane, propane, n-butane, and n-octane with substantial error reduction. Third, the Kihara parameters are tuned based on the reported $T_{{{\text{cp}}}}$ empirical reports to have the most likely precise predictions. Finally, the Peng–Robinson (PR) EOS is developed to evaluate and discuss the pressure–volume (P–V) diagrams of fluids confined in nanopores with respect to reduction in pore size. In general, 80 data points are examined to evaluate the capabilities of the proposed model, which is one of the most significant collections of critical temperature shift points ever gathered from research articles (Jiang et al. 2004; Jin and Nasrabadi 2016; Lowry and Piri 2018; Morishige and Shikimi 1998; Morishige et al. 2003; Pitakbunkate et al. 2016; Singh et al. 2009; Tan et al. 2019). Previous studies (Jin et al. 2013; Song et al. 2020; Yang et al. 2019; Yang and Li 2020), however, employed much fewer data points for developing their model, among which Yang and Li (2020) had the largest number of data points, i.e., 53. It is noteworthy that the data used in their work belonged, to a large extent, to methane and the diversity of studied fluids was not wide at all.

This study will address the following major research gaps:

• Developing the vdW EOS for confined fluids using a more physically realistic potential function

• Presenting a fluid-dependent expression for critical temperature shift modeling

• Applying a thermodynamic-based adjustment method to acquire a more accurate method to calculate pore critical temperature

Model and theory

At the first step, the conventional vdW EOS must be developed in a modified form to describe the confinement phenomenon. Although this EOS is a simple, it has the qualities and abilities to properly predict the critical and vapor–liquid equilibrium (VLE) properties, i.e., the properties of matters in liquid and gaseous phases at equilibrium. As shown in Fig. 1, we consider a Cartesian pore model consisting of confined particles interacting through any pair potential. For this model, the pressure $\vec{P}$ [${\text{Pa}}$] is a diagonal tensor with components $\vec{P}\left( {p_{i} , i = x,y,z} \right)$. We start with the Gibbs free energy ($G)$ [${\text{J}}$] of the fluids which is given by (Smith et al. 2005),

$$G\left(p,T\right)=U+{P}_{i}V-\mathrm{TS}$$

(1)

where $U$ is the internal energy [${\text{J}}$]; $P$ is the system pressure [${\text{Pa}}$]; $T$ is the system temperature [${\text{K}}$]; $V$ is volume of the system [${\text{m}}^{{3}}$]; and $S$ is the entropy [${\text{J}}{\text{.K}}^{{ - 1}}$]. Besides, the equation for internal energy is as follows (Borgnakke and Sonntag 2012),

Fig. 1

Schematic of a nano-scale pore in Cartesian coordinates

$$\mathrm{d}U=T\mathrm{d}S-{P}_{i}\mathrm{d}V$$

(2)

The pressure expression of vdW EOS can be derived using the derivate of Helmholtz free energy with respect to volume. The Helmholtz free energy of a system including $N$ particles interacting by a pair potential $U\left({r}_{12}\right)$ can be derived from the perturbation theory as (Hansen and McDonald 1990)

$$F=\mathrm{f}\left(T\right)-\frac{kT{N}^{2}}{2{V}^{2}} \iint \left({e}^{-\beta U({r}_{12})}-1\right)\mathrm{d}{V}_{1}\mathrm{d}{V}_{2}$$

(3)

where ${\text{f}}\left( T \right)$ is the free energy of ideal gas [${\text{J}}$]; $k$ is the Boltzmann constant [${\text{J}}{\text{.K}}^{{ - {1}}}$]; and $\beta$ is $1/kT$ [${\text{J}}^{ - 1}$]. The term $e^{{ - \beta U\left( {r_{12} } \right)}} - 1$ is often called the Mayer function. Here, in spite of the former work (Zarragoicoechea and Kuz 2002), we utilized the Kihara potential (Fig. 2) for the interactions between particles. In the Lennard–Jones theory, two molecules can completely interpenetrate provided that they have enough energy (Prausnitz et al. 1998), whereas according to the Kihara model, the repulsion interaction in the rigid spherical core is split into two parts; a rigid core and an outer soft spherical repulsion region (Reed and Gubbins 1973). Therefore, as Fig. 3 shows, we suppose that soft-penetrable electron clouds surround a hard-impenetrable core whose mathematical expression would be as (Prausnitz et al. 1998),

$$U\left( {r_{12} } \right) = \left\{ {\begin{array}{*{20}c} \infty \\ {4\varepsilon \left[ {\left( {\frac{{\sigma_{k} - 2a_{{\text{k}}} }}{{r - 2a_{{\text{k}}} }}} \right)^{12} - \left( {\frac{{\sigma_{{\text{k}}} - 2a_{{\text{k}}} }}{{r - 2a_{{\text{k}}} }}} \right)^{6} } \right]} \\ \end{array} } \right. \begin{array}{*{20}c} {\begin{array}{*{20}c} {r < 2a_{{\text{k}}} } \\ \\ \end{array} } \\ {r \ge 2a_{{\text{k}}} } \\ \end{array}$$

(4)

where $a_{{\text{k}}}$ is the radius of the spherical core [${\text{m}}$]; $\varepsilon$ [${\text{J}}$] is the depth of the energy well (we employ $\varepsilon_{k} = \varepsilon /k$ [${\text{K}}$] as the third Kihara parameter instead of $\varepsilon$ in this study); and $\sigma_{k}$ is the collision diameter [$m$], i.e., the distance $r$ [${\text{m}}$] between molecular centers when $U\left( {r_{12} } \right) = 0$. The intermolecular distance goes from $r_{12} = 2a_{{\text{k}}}$ (where the repulsion force is infinity) to $r_{12} = \infty$.

Fig. 2

Kihara potential function with three adjustable parameters

Fig. 3

A specific position of two molecules in which the value of the Kihara potential function approaches to infinity

Given the Eq. (3), the standard vdW equation can be derived by integrating over an infinite volume. For a finite volume, the same procedure can be applied by splitting the integral into two regions, $2a_{{\text{k}}} < r_{12} < \sigma_{{\text{k}}}$ and $r_{12} > \sigma_{k}$:

$$F=f\left(T\right)-\quad \frac{{\mathrm{kTN}}^{2}}{2{V}^{2}} \left[\underset{2a<{r}_{12}<\sigma }{\overset{ }{\iint }}\left({e}^{-\beta U({r}_{12})}-1\right)\mathrm{d}{V}_{1}\mathrm{d}{V}_{2}+ \underset{{r}_{12}>\sigma }{\overset{ }{\iint }}\left({e}^{-\beta U({r}_{12})}-1\right)\mathrm{d}{V}_{1}\mathrm{d}{V}_{2}\right]$$

(5)

Assuming two approximations for each integral, $e^{{ - \beta U\left( {r_{12} } \right)}} \cong 0$ for the former because $U\left( {r_{12} } \right) \to \infty$ and $e^{{ - \beta U\left( {r_{12} } \right)}} \cong 1 - U\left( {r_{12} } \right)/{\text{kT}}$ for the latter, Eq. (5) becomes

$$F = {\text{f}}\left( T \right) + \frac{{{\text{kTN}}^{2} }}{V}b + \frac{{{\text{kTN}}^{2} }}{{2V^{2} }}\iint\limits_{{r_{12} > \sigma }} {\frac{{U\left( {r_{12} } \right)}}{{{\text{kT}}}}}{\text{d}}V_{1} {\text{d}}V_{2}$$

(6)

with $b = \frac{2\pi }{3}\left( {\sigma_{k}^{3} - \left( {2a_{{\text{k}}} } \right)^{3} } \right)$ and $V = L_{x} L_{y} L_{z}$ in which $L_{x} = L_{y} = L$ and $L_{z} = L_{z}$ for a pore model in the Cartesian coordinates.

The expression for $b$ – that is slightly different between the Kihara and L-J potentials – comes from the analytical solution of the first integral in Eq. (5). Since $b$ represents the repulsive term, its form clearly manifests the effect of the hard-core during the repulsion. The integral in Eq. (6) can be semi-analytically calculated as,

$$\frac{1}{V}\underset{{r}_{12}>\sigma }{\overset{ }{\iint }}\frac{U({r}_{12})}{\mathrm{kT}}\mathrm{d}{V}_{1}\mathrm{d}{V}_{2}=\frac{4\varepsilon }{kT}{\sigma }_{\mathrm{k}}^{3}I(A,{a}_{\mathrm{k}}^{*})$$

(7)

$$I\left(A,{a}_{k}^{*}\right)={C}_{0}\left({a}_{k}^{*}\right)+\frac{{C}_{1}\left({a}_{k}^{*}\right)}{\sqrt{A}}+\frac{{C}_{2}\left({a}_{k}^{*}\right)}{A}$$

(8)

where $A = L_{x} L_{y} /\sigma_{{\text{k}}}^{2}$ is the reduced area of the square cross section of a pore, $a_{k}^{*}$ is the reduced Kihara parameter ($2a_{{\text{k}}} /\sigma_{{\text{k}}}$), and $C_{0}$ to $C_{2}$ are coefficients obtained from the numerical integration. Yet, $C_{0}$. has an analytical expression (Eq. (9)) that is achieved from the integral when $A = \infty$. Besides, $C_{1}$ and $C_{2}$ can be obtained from nonlinear fitting for each fluid that will be presented in Sect. "Results and discussion".

$${C}_{0}\left({a}_{\mathrm{k}}^{*}\right)=\frac{2\pi }{495}\left({a}_{\mathrm{k}}^{*}-1\right)\left(220+{a}_{\mathrm{k}}^{*}\left(31{a}_{\mathrm{k}}^{*}-143\right)\right)$$

(9)

The Helmholtz free energy of a confined fluid is formulated by Eq. (10) derived through applying the ideal gas free energy (${\text{f}}\left( T \right) = - {\text{NkT}}\left[ {\ln V/N + c} \right]$ (Hill 1986)), writing $C_{0} = - a_{k} /2\varepsilon \sigma_{k}^{3}$, and considering the limited compressibility of matter via the substitution $\ln V - \left( \frac{N}{V} \right)b\sim \ln \left( {V - Nb} \right)$.

$$F=\mathrm{f}\left(T\right)-\mathrm{NkTln}\left(V-Nb\right)\quad+2\frac{{N}^{2}}{V}\varepsilon {\sigma }_{\mathrm{k}}^{3}\left(-\frac{{a}_{\mathrm{k}}}{2\varepsilon {\sigma }_{k}^{3}}+\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{A}}+\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{A}\right)$$

(10)

In the vdW theory, a molecule’s volume (volume excluded from the system) is represented by variable $b$ that is, in other words, the available space for molecules to overlap. With this definition, other molecules will have an accessible volume of $V - Nb$ to move in (Johnston 2014). This is the backbone of the idea laying behind the substitution for Eq. (10) which is concretely independent of the applied potential function.

As shown in Fig. 1, a cross section confinement occurs in an axially infinite pore for which the axial ($P_{zz}$) and transverse ($P_{xx}$ and $P_{yy}$) components of the pressure tensor $\vec{P}$ are as follows,

$${P}_{xx}={P}_{yy}=-\frac{{\sigma }_{\mathrm{k}}^{2}}{A}\frac{\partial F}{\partial L}$$

(11)

$${P}_{zz}=-\frac{{\sigma }_{k}^{2}}{A}\frac{\partial F}{\partial A}$$

(12)

Using Eqs. (11) and (12), and applying Eq. (10), the axial and transverse components are given by Eq. (13) and (14).

$${P}_{xx}={P}_{yy}=\frac{\mathrm{NkT}}{V-Nb}-\frac{{N}^{2}}{{V}^{2}} \left[a-\varepsilon {\sigma }_{\mathrm{k}}^{3}\left(3\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{A}}+4\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{A}\right)\right]$$

(13)

$$P_{zz} = \frac{{{\text{NkT}}}}{V - Nb} - \frac{{N^{2} }}{{V^{2} }} \left[ {a - 2\varepsilon \sigma_{{\text{k}}}^{3} \left( {\frac{{C_{1} \left( {a_{{\text{k}}}^{*} } \right)}}{\sqrt A } + \frac{{C_{2} \left( {a_{{\text{k}}}^{*} } \right)}}{A}} \right)} \right]$$

(14)

Clearly, when $A \to \infty$, the bulk vdW equation is obtained and $P_{xx} = P_{yy} = P_{zz}$. In reduced coordinates, Eqs. (13) and (14) can be written as,

$${P}_{xx}^{*}={P}_{yy}^{*}=\frac{{T}^{*}}{{v}^{*}-{b}^{*}}-\frac{{a}^{*}-\left(3\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{A}}+4\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{A}\right)}{{{v}^{*}}^{2}} ,$$

(15)

$${P}_{zz}^{*}=\frac{{T}^{*}}{{v}^{*}-{b}^{*}}-\frac{{a}^{*}-2\left(\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{A}}+\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{A}\right)}{{{v}^{*}}^{2}} ,$$

(16)

where $P^{*} = P\sigma_{{\text{k}}}^{3} /\varepsilon$; $T^{*} = {\text{kT}}/\varepsilon$; $v^{*} = \left( {V/N} \right)\sigma_{{\text{k}}}^{ - 3}$; $b^{*} = b\sigma_{{\text{k}}}^{ - 3}$; and $a^{*} = a/\varepsilon \sigma_{{\text{k}}}^{3}$.

From the theory of corresponding states and the continuity principle of gaseous and liquid phases, the following conditions are required for vdW type EOSs at the critical point (Danesh 1998):

$${\left(\frac{\partial {P}_{zz}^{*}}{\partial {v}^{*}}\right)}_{{T}^{*}={T}_{\mathrm{c}}^{*}}={\left(\frac{{\partial }^{2}{P}_{zz}^{*}}{\partial {{v}^{*}}^{2}}\right)}_{{T}^{*}={T}_{\mathrm{c}}^{*}}=0$$

(17)

Applying the requirements, we are furnished with the critical properties’ equations.

$${T}_{c}^{*}=\frac{8}{27{b}^{*}}\left({a}^{*}-2\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{A}}-2\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{A}\right)$$

(18)

$${P}_{\mathrm{c}}^{*}=\frac{{a}^{*}-2\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{A}}-2\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{A}}{27{{b}^{*}}^{2}}$$

(19)

$$v_{{\text{c}}}^{*} = 3b^{*}$$

(20)

Correspondingly, the above modeling can be applied to nano-scale pores in cylindrical coordinates (Fig. 4) in which pores are infinite over the $Z$ direction. This is what has been considered for all calculations in this work, similar to the previous work (Zarragoicoechea and Kuz 2002) in which the square section area in Eq. (18) is equated to the area corresponding to the experimental cylindrical pore, $A=\pi {\left({r}_{\mathrm{p}}/{\sigma }_{\mathrm{k}}\right)}^{2}$. Equations (18) and (19) reveal a shift in the critical temperature and pressure in comparison with those of the bulk vdW equation.

Fig. 4

Schematic of a nano-scale pore in cylindrical coordinates

Since molecule-wall interactions always lead to forming of a layer onto the pore walls (adsorption thickness, $\delta_{{\text{p}}}$ [${\text{nm}}$]), the available volume for free fluid molecules to move will shrink. To consider adsorption thickness in our modelling, we use some measured data for various fluids from the literature, i.e., adsorption thickness at different pore sizes. Since these reported adsorption thicknesses are not necessarily measured at our desired pore radii, we attempt to find a correlation, exclusively for each component, that would have the least possible deviation from the measured data.

Mathematically, there will be a slight change in Eq. (18) and the areal cross section would decrease due to the adsorption thickness. Therefore, the critical temperature and pressure shift equations turn into:

$${T}_{\mathrm{c}}^{*}=\frac{8}{27{b}^{*}}\left({a}^{*}-2\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{\pi }\left(\frac{{r}_{\mathrm{p}}-{\delta }_{\mathrm{p}}}{{\sigma }_{\mathrm{k}}}\right)}-2\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{\pi {\left(\frac{{r}_{p}-{\delta }_{p}}{{\sigma }_{k}}\right)}^{2}}\right)$$

(21)

$${P}_{\mathrm{c}}^{*}=\frac{{a}^{*}-2\frac{{C}_{1}\left({a}_{\mathrm{k}}^{*}\right)}{\sqrt{\pi }\left(\frac{{r}_{\mathrm{p}}-{\delta }_{\mathrm{p}}}{{\sigma }_{\mathrm{k}}}\right)}-2\frac{{C}_{2}\left({a}_{\mathrm{k}}^{*}\right)}{\pi {\left(\frac{{r}_{\mathrm{p}}-{\delta }_{\mathrm{p}}}{{\sigma }_{\mathrm{k}}}\right)}^{2}}}{27{{b}^{*}}^{2}}$$

(22)

where $\delta_{{\text{p}}}$ [${\text{nm}}$] is the adsorption thickness for which here we employ the values Zhang et al. (2019a; b) have reported in their work. In this study, instead of a logarithmic function proposed in their paper, a rational model as Eq. (23) is used for least square fitting that results in less deviation:

$${\delta }_{p}=\frac{a+b{\times r}_{p}}{1+c\times {r}_{p}+d\times {r}_{p}^{2}}$$

(23)

in which $r_{{\text{p}}}$ and $\delta_{{\text{p}}}$ are both in ${\text{nm}}$. Moreover, $a$, $b$, $c$, and $d$ are coefficients whose values are listed in Table 1 for N₂, CO₂, and normal alkanes from CH₄ to n-C₁₀H₂₂. Figure 5 shows the measured (Zhang et al. 2019a, b) and fitted values of $\delta_{{\text{p}}}$. Clearly, there is a decreasing trend as pore radius increases, implying that the effect of adsorption thickness is less at larger pores.

Table 1 Values of coefficients in the adsorption thickness equation utilized in this work

Fig. 5

Measured and fitted values of adsorption thickness

Furthermore, the following expressions may be used for normal alkanes based on the values presented in Table 1:

$$a = 0.12{ }\ln \omega + 1.0495\quad R^{2} = \, 0.9235$$

(24)

$$b = 0.0225{ }\ln \omega + 0.1332\quad R^{2} = \, 0.8675$$

(25)

$$c = - 1.0654 \omega^{2} + 0.1643 \omega + 0.6555\quad R^{2} = \, 0.7315$$

(26)

$$d = 6.3464 \times 10^{ - 6} \ln \omega - 5.9508 \times 10^{ - 6} \quad R^{2} = \, 0.8454$$

(27)

Considering the adsorption effect, the critical shift phenomena is now modeled better and more realistic than the previous work (Zarragoicoechea and Kuz 2002) which was developed for neutral pore walls.

Results and discussion

Hysteresis critical temperature, ${{\varvec{T}}}_{{\varvec{c}}{\varvec{h}}}$

As illustrated in Fig. 6 (Morishige and Shikimi 1998), Hysteresis Critical Temperature ($T_{{{\text{ch}}}}$) [${\text{K}}$] is the temperature above which the hysteresis loop vanishes in adsorption–desorption diagrams. Moreover, the jump in the adsorption curve associated with capillary condensation decreases as the temperature rises and eventually disappears at Pore Critical Temperature ($T_{{{\text{cp}}}}$) [${\text{K}}$] that is higher than $T_{{{\text{ch}}}}$ (Morishige and Shikimi 1998). As it has been followed previously (Zarragoicoechea and Kuz 2002), Eqs. (21) and (22) are supposed to calculate $T_{{{\text{ch}}}}$.

Fig. 6

Adsorption/desorption volume versus relative fugacity for argon at 2.1 nm

Table 2 shows the values of $C_{0}$ to $C_{2}$ (coefficients of Eq. (8)) for argon, nitrogen, oxygen, carbon dioxide, xenon, and ethylene. To make an in-depth clarification, the numerical values and their nonlinear fitted curves are presented in Fig. 7 and compared to those of Lennard–Jones (Zarragoicoechea and Kuz 2002).

Table 2 Values of ${C}_{0}$ to ${C}_{2}$ for argon, nitrogen, oxygen, carbon dioxide, xenon and ethylene

Fig. 7

Numerical solution to Eq. (8) for Ar, O₂, C₂H₄, CO₂ and their fitted curves as well as that of the Lennard–Jones model

Calculated $T_{{{\text{ch}}}}$ values from this study are listed in Table 3 and compared with the corresponding L-J outputs as well as the empirical reports. Employing the Kihara potential function, we observe that results are more satisfactory and much closer to the experimental values, whereas the L-J model has led to far inaccurate predictions. For all components studied in this section, the overall error improvement was almost 9 percent on average. The highest and lowest improvements belonged to CO₂ (19.67%) and Ar (0.18%), respectively. Nevertheless, CO₂, O₂, and Xe still have great deviation from empirical measurements of $T_{{{\text{ch}}}}$, even though our model is more physically realistic than the L-J model. Two factors can account for such a deviation; uncertainties in experimental reports, and incapability of both potential functions to express extreme flexibility in modeling. Regarding the latter, although the Kihara potential is much more flexible in terms of mathematical modeling, both potential functions are merely a representative and approximation of the real potential behavior of two interactive molecules (Prausnitz et al. 1998). Hence, they may lead to erroneous results for some fluids whose true potential behavior is far different from them. On the other hand, $T_{{{\text{ch}}}}$ is a concept to which only a few research works have turned their attention (Morishige and Ito 2002; Morishige and Shikimi 1998; Morishige et al. 2003), and there is no crystal-clear explanation for its true nature.

Table 3 Calculated hysteresis critical temperatures of Ar, N₂, O₂, C₂H₄, CO₂, and Xe plus those of the Lennard–Jones model and the experimental reports

Pore critical temperature,${\varvec{T}}_{{{\text{cp}}}}$

We followed the same procedure of the previous section herein for nitrogen, carbon dioxide, and normal alkanes including methane, ethane, propane, n-butane, and n-octane. Due to having a decent molecular weight disparity, these components are a superb representative of n-alkanes and their behavior to be studied in nano-scale pores. As we discussed earlier, the essential point is that the application of Kihara potential provides specific values of $C$ coefficients for each type of fluid (Table 4) rather than the unique constants of the Lennard–Jones. For convenience, Eqs. (28) and (29) are proposed as a function of ${C}_{0}$ for values in Table 4, with ${R}^{2}=1$.

Table 4 Values of C₀ to C₂ for CH₄, C₂H₆, C₃H₈, n-C₄H₁₀, and n-C₈H₁₈

$${C}_{1}\left({a}_{\mathrm{k}}^{*}\right)=0.158\times {C}_{0}{\left({a}_{\mathrm{k}}^{*}\right)}^{2}-1.3081\times {C}_{0}\left({a}_{\mathrm{k}}^{*}\right)-0.2264$$

(28)

$${C}_{2}\left({a}_{\mathrm{k}}^{*}\right)=-0.1485\times {C}_{0}{\left({a}_{\mathrm{k}}^{*}\right)}^{2}+0.4126\times {C}_{0}\left({a}_{\mathrm{k}}^{*}\right)+0.1902$$

(29)

When the pore cross section goes to infinity and the fluid is not confined anymore, it is totally expected that Eqs. (21) and (22) must result in the bulk critical temperature ($T_{{\text{C}}}$) and pressure ($P_{{\text{c}}}$), or $T_{{{\text{cp}}}} \left( {r = \infty } \right) = T_{{\text{C}}}$ and $P_{{{\text{cp}}}} \left( {r = \infty } \right) = P_{{\text{c}}}$. Nonetheless this will not take place and the approach developed for $T_{ch}$ must be adjusted, accordingly. Although this issue is tackled via a critical shift correlation for the L-J model (Zarragoicoechea and Kuz 2004), here we offer and follow a thermodynamic-based perspective in which $a^{*}$ and $b^{*}$ must be equal to the vdW constants $a_{{{\text{vdW}}}}$ and $b_{{{\text{vdW}}}}$, respectively. This will manipulate the Kihara parameters until the demanded equality between $T_{{{\text{cp}}}}$ ($P_{{{\text{cp}}}}$) at $r = \infty$ and $T_{{\text{c}}}$ $\left( {P_{{\text{c}}} } \right)$ is achieved.

Being a three-parameters potential function, the Kihara model also needs a third equation for the adjustment procedure, $a_{{\text{k}}}^{*} = {\text{f}}\left( \omega \right)$, where $\omega$ is the acentric factor. Furthermore, it seems that the value of $a_{{\text{k}}}^{*}$ remains fixed for every fluid, and it comes from one of the previous works for Kihara parameters (Tee et al. 1966a, b) in which the value of $2{a}_{\mathrm{k}}/({\sigma }_{k}-2{a}_{\mathrm{k}})$ has been determined through numerical curve fittings of $B(T)$ and/or viscosity ($\mu$) data.

Adjusted Kihara parameters of normal alkanes from CH₄ to n-C₈H₁₈ are brought in Table 5.

Table 5 Adjusted Kihara parameters of normal alkanes

The microscopic and macroscopic theories of the corresponding states could be related by establishing a connection between one theory's parameters and those of the other. In the former theory, there are two independent parameters (energy and distance parameters), and in the latter, there appear to be three ($v_{{\text{c}}}$, $T_{{\text{c}}}$, and $P_{{\text{c}}}$). Yet, in fact, only two of these are independent since the compressibility factor at the critical point $\left( {z_{{\text{c}}} } \right)$ is the same for all fluids (Prausnitz et al. 1998). For a fluid at a characteristic state, where the liquid and the gaseous phases become identical, the critical temperature is a measure of kinetic energy; Hence, a relation between $\varepsilon _{{\text{k}}}$ and $T_{{\text{c}}}$ seems quite reasonable. Identically, the critical molar volume ($v_{{\text{c}}}$) represents the molecular size, and proportionality between $\sigma _{{\text{k}}}^{3}$ or $a_{{\text{k}}}^{3}$ and $v_{{\text{c}}}$ also seems to be reasonable (Prausnitz et al. 1998). Thus, Eqs. (30) to (32) are developed to be used as alternative correlations for the values listed in Table 5, with the average ${R}^{2}$ of 0.993,

$$\frac{{a}_{\mathrm{k}}}{\sqrt[3]{{v}_{\mathrm{c}}}}=-26.834{\omega }^{3}+20.26{\omega }^{2}+0.1882\omega +0.743$$

(30)

$$\frac{{\sigma }_{\mathrm{k}}}{\sqrt[3]{{v}_{\mathrm{c}}}}=6.3389{\omega }^{3}+1.1671{\omega }^{2}+0.6967\omega +7.0352$$

(31)

$$\frac{{\varepsilon _{{\text{k}}} }}{{T_{{\text{c}}} }} = - 28.888\omega ^{3} + 24.225\omega ^{2} - 0.336\omega + 1.8308$$

(32)

where ${T}_{c}$ and ${v}_{c}$ are the bulk critical temperature [$\mathrm{K}$] and volume [${\mathrm{m}}^{3}.{\mathrm{Kgmol}}^{-1}$], respectively, and $a_{{\text{k}}}$ and $\sigma _{{\text{k}}}$ both are in angstrom.

Tables 6, 7, 8 present the calculated pore critical temperatures of N₂, CO₂, CH₄, C₂H₆, C₃H₈, n-C₄H₁₀, and n-C₈H₁₈, and compares them to those of the Lennard–Jones as well as the literature reported values. Therein, this study has ameliorated the average overall error of normal alkanes by 6.4%. The overall mean error of this work was less than 2.5%, whereas that of the Lennard–Jones was 8.11%. Note that the highest and the lowest improvements belong to heavier normal alkanes (n-C₄H₁₀ and n-C₈H₁₈) and non-hydrocarbons, respectively.

Table 6 Comparison among the calculated ${\mathrm{T}}_{\mathrm{cp}}$ values of CH₄, n-C₄H₁₀, and n-C₈H₁₈ based on their adjusted Kihara parameters as well as those of the adjusted L-J model and the reported data

Table 7 Comparison among calculated ${\mathrm{T}}_{\mathrm{cp}}$ of C₂H₆ from this study, adjusted L-J, and the reported values

Table 8 Comparison among calculated ${\mathrm{T}}_{\mathrm{cp}}$ of N₂ from this study, the adjusted L-J, and the reported values

Overall, it is noticeable that in smaller pore radii, ${T}_{\mathrm{cp}}$ values calculated by the L–J model decrease dramatically, while our model has a more rational reduction. Such a reasonable behavior accounts for our lower deviation in smaller pore sizes for most of the fluids.

For CH₄, C₂H₆, n-C₄H₁₀, and n-C₈H₁₈, ${T}_{\mathrm{cp}}$ reports were available from various research studies. Regarding the values from Singh et al. (2009), generally, our calculations were closer to $T_{{{\text{cp}}}}$ values of mica porous media. Regarding CH₄, however, the least deviation was achieved when we compared our calculations to the reports from the work of Jin and Nasrabadi (2016), wherein the porous media was, interestingly, made out of graphite. Although the systems studied in different works have different conditions, the results for CH₄ denotes that $T_{{{\text{cp}}}}$ values generated by our model is closer to graphite porous media in larger pore radii, and to those of mica in smaller ones.

According to the ${T}_{\mathrm{cp}}$ reports of C₂H₆ for which the largest number of data points were available, our calculations were closer to measurements except for one pore radius, 3 nm, reported by Pitakbunkate et al. (2016) for graphite porous media. Considering all the data points from their study, the value of pore critical temperature witnessed a drastic reduction from 4 to 3 nm, which was a behavior similar to that of the L–J model in smaller pore sizes.

As a general note, both models underestimate a fluid’s pore critical temperature. Yet, this underestimation is much more severe in the L–J model and the approach proposed in this study succeeded in addressing this issue by employing a more physically realistic potential function and a better adjustment method to calculate ${T}_{\mathrm{cp}}$.

Figures 8 and 9 reveal for the assessed fluids that how much this study generated better results, compared with the adjusted L-J model. Besides, the adjusted Kihara parameters of components in Sect. "Results and discussion" are presented in Table 9. Equations (33) to (35) may also be alternatively employed with an average ${R}^{2}$ of 0.9983 for the values of this table, excluding xenon.

Fig. 8

Cross plots of calculated ${T}_{\mathrm{cp}}$ versus measured ${T}_{\mathrm{cp}}$ for CH₄, C₂H₆, n-C₄H₁₀, and n-C₈H₁₈. Filled symbols represent this work

Fig. 9

Cross plots of calculated ${T}_{\mathrm{cp}}$ versus measured ${T}_{\mathrm{cp}}$ for N₂, CO₂, and C₃H₈. Filled symbols represent this work

Table 9 Adjusted Kihara parameters of Ar, N₂, O₂, CO₂, C₂H₄ and Xe

$$\frac{{a}_{\mathrm{k}}}{\sqrt[3]{{v}_{\mathrm{c}}}}=10.7365-\frac{12.854}{1+\omega }+\frac{2.8788}{1+9.4694\omega -0.00075{\omega }^{2}}$$

(33)

$$\frac{{\sigma }_{\mathrm{k}}}{\sqrt[3]{{v}_{\mathrm{c}}}}=8.67913-\frac{1.7357}{1+\omega }+\frac{1.9238}{1+1933.51\omega +33.84{9\omega }^{2}}$$

(34)

$$\frac{{\varepsilon }_{\mathrm{k}}}{{T}_{\mathrm{c}}}=9.2272-\frac{8.3714}{1+\omega }+\frac{2.99046}{1+111.332\omega +0.497007{\omega }^{2}}$$

(35)

Tuning versus curve fitting

Phase behavior models based on the aforementioned equations may not lead to accurate predictions, even for a well characterized (adjusted) model proposed in the previous section. Hence, in order to have precise results, we will calibrate or tune the model’s parameters against experimental data generated at pertinent conditions for specific studies. Here we manually tune the Kihara parameters for n-alkanes studied in this research, predicting the reported pore critical temperatures with the maximum precision. In favor of three adjustable parameters, the Kihara potential function is entirely flexible to be manipulated for almost all fluids (Tee et al. 1966a, b). Nevertheless, this approach affects the bulk critical temperature calculations when the pore size goes to infinity (details in Supporting Information). Table 10 shows the tuned Kihara parameters for N₂, CO₂, CH₄, C₂H₆, n-C₄H₁₀, and n-C₈H₁₈. Noteworthy it is that the Kihara parameters of CH₄, n-C₄H₁₀, and n-C₈H₁₈ are tuned based on their Mica type ${T}_{cp}$ measurements (Singh et al. 2009) (because of more sensible trend), and those of C₂H₆, CO₂, and N₂ based on the experimental measurements (Morishige and Shikimi 1998; Morishige et al. 2003; Tan et al. 2019).

Table 10 Tuned Kihara parameters of CH₄, C₂H₆, n-C₄H₁₀, and n-C₈H₁₈

With the tuned parameters, the pore critical temperature was calculated once again. From results shown in Tables 11 and 12, the overall error of our approach for n-alkanes was near 0.4%.

Table 11 Calculated pore critical temperatures of CH₄, C₂H₆, n-C₄H₁₀, and n-C₈H₁₈ based on their tuned Kihara parameters

Table 12 Calculated pore critical temperatures of N₂, CO₂, and C₂H₆ based on their tuned Kihara parameters

In recent years, various measures have been taken in order to optimize the Lennard–Jones critical temperature shift equation via curve fitting of experimental values in the literature. The most renowned ones are as follows:

$$1 - \frac{{T_{{{\text{cp}}}} }}{{T_{{\text{c}}} }} = 0.8493{\text{~}}\left( {\frac{{r_{{\text{p}}} }}{{\sigma _{{{\text{LJ}}}} }}} \right)^{{ - 1.241}} + 0.015\quad ({\text{Jin}}\;{\text{et}}\;{\text{al}}.\;2013)$$

(36)

$$1 - \frac{{T_{{{\text{cp}}}} }}{{T_{{\text{c}}} }} = 3.374 \times \left( {\frac{{r_{{\text{p}}} }}{{\sigma _{{{\text{LJ}}}} }}} \right)^{{ - 1.637}} \quad ({\text{Yang}}\;{\text{et}}\;{\text{al}}.\;2019)$$

(37)

$$1 - \frac{{T_{{{\text{cp}}}} }}{{T_{{\text{c}}} }} = 1.7391 \times \left( {\frac{{r_{{\text{p}}} }}{{\sigma _{{{\text{LJ}}}} }}} \right)^{{ - 1.379}} \quad ({\text{Yang and Li 2}}0{\text{2}}0)$$

(38)

$$1 - \frac{{T_{{{\text{cp}}}} }}{{T_{{\text{c}}} }} = 0.6794 \times \left( {\frac{{r_{{\text{p}}} }}{{\sigma _{{{\text{LJ}}}} }}} \right)^{{ - 0.7878}} \quad ({\text{Song}}\;{\text{et}}\;{\text{al}}.\;2020)$$

(39)

The flexibility of Kihara potential provided us with the most precise predictions in favor of our tuning approach. Despite numerical curve fittings in which only the constants in the critical shift equation vary, here for the very first time, we dedicated a specific set of values for each component (Table 10), trying to make the results perfectly matched to the real pore critical temperatures (Figs. 10 and 11). Yet, ${a}_{\mathrm{k}}^{*}$ values has been kept equal to its reported ones (Tee et al. 1966a, b) in order to achieve rational outcomes. Although our adjusted Kihara method is only based on thermodynamic adjustments, it is surprisingly able to generate better results rather than the aforementioned L–J correlations (details in the Supporting Information).

Fig. 10

Cross plots of calculated ${T}_{\mathrm{cp}}$ versus measured ${T}_{\mathrm{cp}}$ for CH₄, C₂H₆, n-C₄H₁₀, and n-C₈H₁₈

Fig. 11

Cross plots of calculated ${T}_{\mathrm{cp}}$ versus measured ${T}_{\mathrm{cp}}$ for N₂ and CO₂

Although the correlation presented by Jin et al. (2013) is the oldest model included in this section, it performed better than the newer models of Yang and Li (2020) and Song et al. (2020) almost for all fluids except for nitrogen. This can be attributed to the fact that the diversity of fluids studied in the work of Jin et al. was much wider than other ones. In fact, the recent models were mostly developed based on ${T}_{\mathrm{cp}}$ reports of methane. That is why the four correlations in this section perform roughly the same for CH₄. This is much more conspicuous for the correlation proposed by Song et al., when applied for other fluids than CH₄.

Capillary condensation–Maxwell construction

In thermodynamic equilibrium, a necessary condition for stability is that pressure does not increase with volume. This basic consistency requirement is sometimes violated in analytic models for first order phase transitions, such as that of vdW-type EOSs (Danesh 1998). The occurrence of this violation in P–V diagrams of vdW-type EOSs is known as capillary condensation (Zarragoicoechea and Kuz 2002).

The vdW P–V isotherm diagram of the Lennard–Jones model has been evaluated in the literature (Zarragoicoechea and Kuz 2002). In this section, we investigate whether capillary condensation will form in vdW and PR EOSs’ P–V diagrams when the EOS is associated with the Kihara potential function. As illustrated in Fig. 12 for vdW, the Maxwell construction is utilized to acquire the gas–liquid coexistence dashed line in the two-phase region where ${{P}_{zz}}_{L}={{P}_{zz}}_{V}$. It is worth mentioning that the critical shift is applied in the EOS through substituting T_c= T_cp and P_c= P_cp .

Fig. 12

Capillary condensation of vdW EOS for methane in axial section (${P}_{zz}$) at ${T}_{\mathrm{r}}=0.8$ and ${r}_{p}=2$ nm

Nonetheless, the controversial subject here is the P–V diagrams of PR EOS as the most widely applied equation in the industry (Danesh 1998). It is has been reported that a confined fluid will be in its supercritical state when PR is coupled with L–J (Islam et al. 2015), whereas the vdW EOS gives capillary condensation at the same bulk reduced temperature (${T}_{\mathrm{r}}=T/{T}_{\mathrm{c}}$); therefore, the predicted state of the fluid is not supercritical.

To address this for the Kihara potential function, first, the critical-shift-associated form of PR EOS must be driven, which is explained as follows (Peng and Robinson 1976):

$$P=\frac{\mathrm{RT}}{v-{b}_{\mathrm{PR}}}-\frac{{a}_{\mathrm{PR}}}{v\left(v+{b}_{\mathrm{PR}}\right)+{b}_{\mathrm{PR}}(v-{b}_{\mathrm{PR}})}$$

(40)

where

$${a}_{\mathrm{PR}}=0.45724\frac{{\left({\mathrm{RT}}_{\mathrm{c}}\right)}^{2}}{{P}_{\mathrm{c}}}\alpha (T)$$

(41)

$$b = 0.0778\frac{{RT_{{\text{C}}} }}{{P_{{\text{c}}} }}$$

(42)

$$\alpha \left( T \right) = \left( {1 + m\left( {1 - \sqrt {T_{{\text{r}}} } } \right)} \right)^{2}$$

(43)

$$m = 0.37464 + 1.54226\omega - 0.26992\omega ^{2}$$

(44)

where $a_{{{\text{PR}}}}$ and $b_{{{\text{PR}}}}$ are EOS constants; v is molar volume [m³. Kgmol⁻¹]; T_c [K] and P_c [MPa] are critical temperature and pressure; T_r is the bulk reduced temperature; w is the acentric factor; andR is the universal gas constant [m³.MPa.K^-1.Kg mol^-1]. To apply the critical shift effect, similar to vdW, T_c and P_c must be recalculated through Eqs. (21) and (22), based on the adjusted Kihara parameters, and C coefficients.

Afterward, we comprehended that the PR EOS, when coupled with the Kihara potential, can demonstrate capillary condensation as the temperature declines. As Fig. 13 reveals, for methane and carbon dioxide at a given reduced temperature and pore size, the confined fluid is homogenous with respect to various axial and transverse pressures ${P}_{zz}$< ${P}_{xx}$.

Fig. 13

Volume dependence of the axial (${P}_{zz}$) and transverse (${P}_{xx}$) pressure tensor components for a methane at ${T}_{\mathrm{r}}=0.94$ and b carbon dioxide at ${T}_{\mathrm{r}}=0.96$ in different pore sizes

Elaborating further, we are well provided with the commencement of the confinement effect on the phase behavior of fluids, designated by the purple-dashed lines in Fig. 13. For instance, at almost the same $T_{r}$, carbon dioxide is under confinement at larger pores (10 to 20 nm) than methane, which feels the critical shift in pore sizes of less than 5 to 10 nm. Such a concept could be directly related to the higher $a_{k}$ value of CO₂, forbidding the molecules to excessively approaching each other rather than the CH₄ ones, which are further allowed in interpenetration. At this step, our results are similar to those of Islam et al. (2015), but further evaluation reveals the differences.

Decreasing the temperature and maintaining the transverse pore area constant where T is less than $T_{{{\text{cp}}}}$ at 2 nm, we observe a loop appearing for the axial pressure $P_{{zz}}$ (Fig. 14). If we take the temperature even lower, the transverse pressure also begins to presents a loop along with the axial part (Fig. 15), which indicates that in contrast to the reports of Islam et al. (2015), both vdW and PR EOSs manifest capillary condensation in the axial part at the same reduced temperature.

Fig. 14

Capillary condensation of Peng-Robinson for a methane at ${T}_{\mathrm{r}}=0.82$ and b carbon dioxide at ${T}_{\mathrm{r}}=0.84$ in the axial section (${P}_{zz}$)

Fig. 15

Capillary condensation of Peng-Robinson for a methane at ${T}_{\mathrm{r}}=0.80$ and b carbon dioxide at ${T}_{\mathrm{r}}=0.82$ in the axial section (${P}_{zz}$) and its tendency to act the same in the transverse one (${P}_{xx}$)

To provide a superior perspective of the phase behavior, Fig. 16 is presented to show and compare the 3D Pressure–Volume-Temperature diagrams of CH₄ and CO₂ for the axial pressure generated by PR EOS, in which the red and purple colors represent the lowest and highest values of axial pressure, respectively.

Fig. 16

P_zz-V-T 3D diagram of methane and carbon dioxide at ${r}_{p}=2$ nm generated by PR-EOS

Conclusions

Critical shift and phase transition of confined fluids have continually been a focal point of interest and many works have been done with the help of molecular thermodynamics science. The model developed in this study benefits from the Kihara potential function, which has never been employed in similar studies by other researchers due to its more complex formulation and high computational cost. This potential is physically more realistic than that of previous studies, i.e., the simple Lennard–Jones potential. Our model led to a more sophisticated as well as accurate form of vdW EOS for confined matters. This, along with the new proposed adjustment method that in contrast to previous works in the literature is thermodynamically reliable, enabled us to have more accurate predictions of critical temperature of fluids in tight porous media. The 80 datapoint used in this work is, to the best of our knowledge, the largest collection of critical shift data points ever gathered, thereby enabling us to more reliably evaluate the performance of our proposed model.

The following major conclusions can be made from this work:

• Comparing with the previous studies, the hysteresis critical temperatures (${T}_{\mathrm{ch}}$) of fluids were predicted with 9 percent less error in this work.

• The Kihara parameters (${a}_{\mathrm{k}}$, ${\sigma }_{\mathrm{k}}$, and ${\varepsilon }_{\mathrm{k}}$) were adjusted in order to accomplish the bulk state requirements of the developed cubic EOS. This new adjustment method was able to calculate the pore critical temperature (${T}_{\mathrm{cp}}$) of fluids with 6.4% less error than previous studies.

• Based on the available T_cp reports of some fluids, their Kihara parameters were tuned. This tuning approach had merely 0.4 percent deviation for n-alkanes.

• Regarding pressure–volume diagrams, the Kihara model demonstrated capillary condensation for both vdW and PR EOSs; however, it was reported that the L-J model does so only for the vdW EOS.

Limitations and suggestions for future works

Since the model in this study is developed particularly for calculating the shift of critical temperature of confined fluids, it is not unexpected to have erroneous predictions of critical pressure shift (${P}_{\mathrm{cp}}$). Hence, it is suggested to consider this issue for future works to achieve a balance between more accurate predictions of ${T}_{\mathrm{cp}}$ as well as sufficiently accurate estimations of $P_{{{\text{cp}}}}$ by using experimental/simulation measurements, if any.