Numerical implementation of the multicomponent potential theory of adsorption in Python using the NIST Refprop database

TL;DR

This paper details a Python-based numerical implementation of the multicomponent potential theory of adsorption, utilizing the NIST Refprop database to accurately model gas mixture adsorption in various regimes.

Contribution

It introduces a comprehensive Python implementation of the multicomponent potential theory of adsorption using NIST Refprop, addressing practical challenges and demonstrating with CH4/CO2 mixtures.

Findings

Successful modeling of CH4/CO2 adsorption isotherms

Identification of limitations in the current model

Application across subcritical and supercritical regimes

Abstract

In this paper, we present a detailed numerical implementation of the multicomponent potential theory of adsorption which is among the most accurate gas mixtures adsorption models. The implementation uses the NIST Refprop database to describe fluid properties and applies to pure gases and mixtures in both subcritical and supercritical regimes. The limitations of the model and the issues encountered with its implementation are discussed. The adsorption isotherms of CH4 / CO2 mixture are modeled and parameterized as implementation examples.

Tables

Table 1. Table 1: Variables definition.

Variables	Type	Definition
${\mu }_B$, chem()	function	Chemical potential of the bulk phase (J/mol)
${\mu }_{Ad}$, chem()	function	Chemical potential of the adsorbed phase (J/mol)
$\epsilon$, DRA()	function	Adsorbent surface potential (J/mol)
${\rho }_B$, db	float	Fluid density in the bulk phase (mol/L)
${\rho }_{Ad}$, d_z	float	Fluid density in the adsorbed phase (mol/L)
${\epsilon }_0$, eps0	float, array	Characteristic energy of adsorption (J/mol)
$z$, z	float	Microporous volume (cm3/g)
$z_0$, z0	float	Limiting micropore volume (cm3/g)
$\beta$	float	Heterogeneity parameter (fix to 2 in this work)
$N_{ex}$, N_ex_	function	Excess (Gibbs) adsorption (mol/Kg)
dataP	array	Experimental values of pressure (KPa)
dataD	array	Corresponding values of fluid density (mol/L)
dataAd	array	Experimental values of excess adsorption (mol/Kg)
$x_B$, xB	array	Bulk phase molar fraction array
$x_{Ad}$, x_z	array	Adsorbed phase molar fraction array
x_z_A , x_z_B	float	Adsorbed phase molar fraction of components
d_max	function	Maximal fluid density allowed by the Refprop (mol/L)
d_vap	float, array	Fluid vapor density at dew point (mol/L)
d_liq	float, array	Fluid liquid density at dew point (mol/L)
X	array	Array containing pure gases arrays
d0	float	Initial guess for fluid density (mol/L)
x0	array	Initial guess for mixture molar fraction array
T	float	Fluid temperature (K)

Table 2. Table 2: Initial and fitted values of the parameters for pure gases.

Parameter	Initial value	fitted value	fit error $\%$
$z_{0,C{H}_4}$ ($c{m}^3/g$)	1	0.269	$0.38\%$
$z_{0,N_2}$ ($c{m}^3/g$)	1	0.256	$0.72\%$
$z_{0,C{O}_2}$ ($c{m}^3/g$)	1	0.306	$2.16\%$
${\epsilon }_{0,C{H}_4}$ ($J/mol$)	$\mathrm{10\hspace{0.17em}000}$	8568	$0.51\%$
${\epsilon }_{0,N_2}$ ($J/mol$)	$\mathrm{10\hspace{0.17em}000}$	7020	$0.72\%$
${\epsilon }_{0,C{O}_2}$ ($J/mol$)	$\mathrm{10\hspace{0.17em}000}$	7541	$3.96\%$

Table 3. Table 3: Precision of the fits for pure gases on Calgon F-400

Gas	Mean error (ADDn)	Experimental uncertainty [23]
$C{H}_4$	0.36$\%$	1.8$\%$
$N_2$	0.59$\%$	2.3$\%$
$C{O}_2$	4.71$\%$	6.4$\%$
Mean	1.89$\%$	3.5$\%$

Table 4. Table 4: Values of the parameters for CH 4  / CO 2 mixture on Calgon fitted from pure gas isotherms.

Parameter	Initial value	Fitted value	Fit error $\%$
$z_{0,common}$ ($c{m}^3/g$)	1	0.300	$1.79$
${\epsilon }_{0,C{H}_4}$ ($J/mol$)	10 000	7475	$3.08$
${\epsilon }_{0,C{O}_2}$ ($J/mol$)	10 000	7767	$3.45$

Table 5. Table 5: Mean error between experimental data and the MPTA model for a 80 % percent 80 80\%  CH 4  /  20 % percent 20 20\%  CO 2 mixture on Calgon.

	Model mean error (ADDn)	Experimental uncertainty [23]
CH4	4.42$\%$	4.0$\%$
CO2	4.41$\%$	4.3$\%$
Mixture	3.67$\%$	3.9$\%$

Equations

$${\mu }_B\left(T,{\rho }_B\right)={\mu }_{Ad}\left(T,{\rho }_{Ad}\right)-\epsilon ,$$

$$\epsilon (z)={\epsilon }_0{\left(\ln \frac{z_0}{z}\right)}^{\frac{1}{\beta }}.$$

$$N_{ex}({\rho }_B)={\int }_0^{z_0}\left({\rho }_{Ad}(z)-{\rho }_B\right)dz.$$

$$f_{Ad}\left(T,{\rho }_{Ad}(z)\right)=f_B\left(T,{\rho }_B\right)\exp \left(\frac{\epsilon (z)}{RT}\right).$$

$${\epsilon }^i(z)={\epsilon }_0^i{\left(\ln \frac{z_0}{z}\right)}^{\frac{1}{\beta }},$$

$${\mu }_B^i({\rho }_B,x_B^i)+{\epsilon }^i(z)-{\mu }_{Ad}^i({\rho }_{Ad}(z),x_{Ad}^i(z))$$

$$N_{ex}^i({\rho }_B)={\int }_0^{z_0}\left({\rho }_{Ad}(z)x_{Ad}^i(z)-{\rho }_B{x}_B^i\right)dz.$$

$$N_{ex}({\rho }_B)=\sum _{i=1}^M{N}_{ex}^i({\rho }_B).$$

$$\displaystyle ADD(N)=\frac{100}{k}\sum_{j=1}^{k}\bigg{|}\frac{N_{ex}^{exp}(\rho_{B}^{j})-N_{ex}^{cal}(\rho_{B}^{j})}{N_{ex}^{exp}(\rho_{B}^{j})}\bigg{|}\%.$$

Full text

Keywords: adsorption; mixture adsorption; multicomponent adsorption; potential theory of adsorption; MPTA; density functional theory

1 Introduction

Gas adsorption in porous material is relevant to a wide variety of industrial and scientific applications, ranging from gas purification and separation to gas storage in stationary or mobile applications. While adsorption experiments of pure gases can be performed readily, experiments on gas mixtures can be challenging, expensive and time-consuming [1]. This situation motivates the utilization of theoretical models to predict the behavior of gas mixtures in the presence of adsorbent.

Several approaches have been proposed to model the adsorption isotherms of both pure gases and mixtures [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. In this paper, we will concentrate on the potential theory of adsorption (PTA), a two-parameter thermodynamic model developed by Shapiro and Stenby [16], based on the pore filling approach of Polanyi’s theory of adsorption [17]. For gas mixtures adsorption, the PTA model generalized into the multicomponent potential theory of adsorption (MPTA), which can deliver great performance, especially when considering disordered adsorbents such as activated carbon (an overview of mixture adsorption models is presented in Ref. [18]). For a $N$ components gas mixture, the simplest MPTA model required the adjustment of only $N+1$ parameters, obtained from pure components adsorption isotherms.

The objective of this paper is to present a detailed numerical implementation of the MPTA model. In this spirit, we will write down explicit lines of code, written in Python3, an increasingly popular, easy to follow, non-compiled and non-proprietary programming language. We will also point out pitfalls and subtleties encounter, and show how to identify and deal with phase transitions.

Definitions and symbols used in this paper are described in Table 1.

2 The multicomponent potential theory of adsorption model

When dealing with adsorption, it is useful to define two regions for the adsorbate, the bulk phase and the adsorbed phase. The bulk phase defines a region far from the adsorbent surface where the fluid is not significantly affected by the presence of the adsorbent material. Conversely, the adsorbed phase defines a region near the adsorbent surface where the fluid is significantly affected by the adsorbent material. The thermodynamic properties of the bulk phase are described by the equation of state (EoS) of the fluid.

2.1 Pure gas

The main assumption of the MPTA model is to suppose that the fluid–surface interaction is completely described by a local potential field $\varepsilon$, generated by the surface. The potential $\varepsilon$ is attractive near the surface of the adsorbent and tends to vanish far from the surface. The constitutive equation of the MPTA model is given by Refs. [15, 16]

[TABLE]

where $\mu_{B}$ and $\rho_{B}$ are the chemical potential and the fluid density in the bulk phase, while $\mu_{Ad}$ and $\rho_{Ad}$ are the locals chemical potential and fluid density in the adsorbed phase. The bulk phase properties are assumed to be constant while the adsorbed phase properties vary with position [16]. Hence, the density $\rho_{Ad}$ is usually very high inside pores, close to the adsorbent, and decrease to reach $\rho_{B}$ far from the adsorbent. Given Eq. (1), the local thermodynamic properties in the adsorbed phase are uniquely determined from properties of the bulk phase and the interaction potential $\varepsilon$. Since we consider only equilibrium situations, the temperature is constant.

Along with Eq. (1), the additional information needed is an explicit form for the potential $\varepsilon$. We use the Dubinin–Radushkevich–Astakhov (DRA) equation [19], a two-parameter ($\varepsilon_{0}$ and $z_{0}$), semi-empirical potential originates from the theory of volume filling micropores [20, 21] and which can represent a wide range of energy distributions

[TABLE]

Here, $\varepsilon_{0}$ is the characteristic energy of adsorption of the adsorbate on the adsorbent, a parameter representing the magnitude of the fluid-solid interaction. The variable $z$ has units of a volume and $z_{0}$ is the limiting micropore volume of the adsorbent which is the volume of the pore that significantly contributes to the adsorption. The ratio $z/z_{0}$ represents the fraction of the microporous volume associated to an energy $\varepsilon(z)$. $\beta$ is a measure of the heterogeneity of the adsorbent and is usually set to 2 for activated carbon, but can eventually become a fitting parameter (in this paper, we will use $\beta=2$). Under certain conditions, the variable $z$ can be seen as the distance from the surface of the adsorbent. This is the case when the microporous structure of the adsorbent exhibit a planar symmetry as in slit-pore adsorbents for example. In this paper, we will use this more intuitive approach and interpret $z$ as the distance from the surface of the adsorbent as in Ref. [16].

The presence of the divergence in Eq. (2) as $z\to 0$ is handled, as shown bellow, by setting a cut-off value for the chemical potential in the numerical implementation.

Eq. (1) now need to be inverted to obtain the local density $\rho_{Ad}(z)$ from the chemical potentials. Then, the (Gibbs) excess adsorption $N_{ex}$ (which is what’s experimentally measured), is calculated by

[TABLE]

We should emphasize a crucial point here: even if $\varepsilon_{0}$ and $z_{0}$ have well-defined physical meaning, these two parameters are better used as fitting parameters in the model. Values for the characteristic energy of adsorption ($\varepsilon_{0}$) for physisorption on microporous adsorbents range from 5–15 kJ/mol while typical limiting micropore volume ($z_{0}$) range from 0.1 to 2 cm3/g.

To fit the model to experimental datasets, initial guess of $\varepsilon_{0}$ and $z_{0}$ are made. The differences between the calculated excess adsorption ($N_{ex}(\rho_{B})$) and the experimental data at various pressures are then simultaneously minimized. The result of the minimization between calculated and experimental values of excess adsorption translates into optimal values for $\varepsilon_{0}$ and $z_{0}$.

In the literature, many authors used the fugacity $f$ instead of the chemical potential. In that case, Eq. (1) takes the form

[TABLE]

Fugacity being essentially the exponential of the chemical potential, we choose to work with the later. Also, it is worth noting that when working with the chemical potential, the logarithmic singularity of Eq. (2) is integrable, which is not the case when working with fugacity.

2.2 Mixtures

For gas mixtures, the situation is more challenging, particularly for the description of the local chemical potential in the adsorbed phase. The molar fraction of the component $i$ will be noted $x^{i}$ such that $\sum_{i}^{M}x^{i}=1$ where $M$ is the number of pure components in the mixture. While the molar fraction of the bulk phase can be well known from experimental measurements (gas chromatography for instance), the molar fraction in the adsorbed phase is unknown and may vary substantially from the bulk phase. In gas mixtures, each gas component interacts with the surface through is own potential

[TABLE]

where $i=1\dots M$. Due to the selectivity of adsorbent material, the adsorption of one component will be favored, affecting the local molar fraction of the mixture in the adsorbed phase.

The introduction of the molar fractions adds unknown variables to the problem. Instead of simply inverting Eq. (1) to find $\rho_{Ad}(z)$, we need to solve simultaneously the following system of equations [22]

[TABLE]

for $x_{Ad}^{i}(z)$ and $\rho_{Ad}(z)$, where the superscript $i$ refer to the gas component of the mixture. From now on, the chemical potential also depend on molar fraction $x^{i}$ of gas components ($x_{B}^{i}$ is the molar fraction of the component $i$ of the bulk phase at equilibrium, i.e. when the adsorption process is complete). In Eq. (6), $\rho_{B}$, $x^{i}_{B}$ and $\varepsilon^{i}(z)$ are known and $\rho_{Ad}(z)$ and $x^{i}_{Ad}(z)$ are the unknowns. The dimension of the system of equation to solve is $M$. The set of $x_{Ad}^{i}(z)$ represent $M-1$ independent variables (because $\sum_{i}^{M}x_{Ad}^{i}=1$ decrease by 1 the number of independent variables) and there is one more independent variable $\rho_{Ad}(z)$ which, of course, is the same for all gas components.

Once the $x_{Ad}^{i}(z)$ and $\rho_{Ad}(z)$ are determined, the next step is to calculate the excess adsorption for each component:

[TABLE]

Finally, the total amount of excess gas adsorbed is simply given by

[TABLE]

As for the pure gases, to adjust the model to experience, initial values of $\varepsilon_{0}^{i}$ and a common $z_{0}$ must be postulated. Because $z_{0}$ is considered to be an adsorbent property (the limiting micropore volume), we use a single value for all gas components. The excess adsorption is calculated for each gas component from the postulated initial values. Then, the differences between the calculated and experimental pure excess adsorption is minimized simultaneously for all experimental bulk phase densities and all gas components. This finally gives the optimal values of parameters $\varepsilon_{0}^{i}$ and a common $z_{0}$.

A key point here is that even for mixtures, the adjustment of the model to experimental data is made using pure gas adsorption isotherms, not mixture data [22]. Once the parameters are adjusted from pure gas isotherms, then gas mixture adsorption can be simulated. It is important, however, to understand that $x_{B}$ should be the bulk phase molar fraction at equilibrium which may differ greatly from the initial or feeding molar fraction. For example, in Ref. [23], the feeding and equilibrium composition can easily differ by more than 25%. Then, using initial instead of equilibrium molar fractions could introduce a significant error in the model predictions. Consequently, one should be careful when using the MPTA model on data with unknown bulk phase equilibrium molar fractions.

2.3 Equation of state

Up to now, there is one thing that we left aside: the way we will evaluate the chemical potentials $\mu^{i}_{B}$ and $\mu^{i}_{Ad}$. In order to accurately describe the thermodynamic properties of the fluid, one should use an accurate equation of state (EoS). A well-known EoS commonly used for this kind of purpose is the Soave–Redlich–Kwong (SRK) EoS [24, 25]. However, in this paper, we choose to use the NIST–REFPROP (Refprop) database [26] which give presumably the most accurate results. For the purpose of this paper, the underlying references used by the Refprop are Refs. [27, 28, 29]. Credit goes to Refs. [30, 31] for the Python integration of Refprop.

2.4 Performance analysis

To quantify the precision of the model, we will use the absolute average deviations (AAD) between experimental and fitted data for excess adsorption:

[TABLE]

In this equation, $\rho_{B}^{j}$ is one of the experimental bulk phase density in the set of $k$ values $\{\rho_{B}^{1},\rho_{B}^{2},\dots,\rho_{B}^{k}\}$. $N_{ex}^{exp}(\rho_{B}^{j})$ and $N_{ex}^{cal}(\rho_{B}^{j})$ are the experimental and calculated gas adsorbed density at bulk phase density $\rho_{B}^{j}$, respectively.

3 Numerical implementation

The proposed numerical integration being not very computationally expensive, Python (version 3) [32, 33] might be the best choice of language to obtain a simple, concise and open-source implementation. Indeed, with all the embedded functions and algorithms (principally from SciPy [34] and NumPy [35] packages), most of the numerical operations are easier to do in Python than in other languages. Moreover, the Jupyter (IPython) notebook [36] and matplotlib [37] package provide great and convenient tools for graphical visualization and analysis.

For readers unfamiliar with Python, the language is simple enough to serve as pseudocode for other languages as well. Also, Section 3.1 give a short introduction to Python. Note that some programming convention will be overlooked for space consideration and that the proposed implementation is oriented on simplicity and concision rather than performance.

3.1 Minimal introduction to Python

We present here some of Python concepts needed for the rest of this paper. Firstly, the required package needs to be loaded. In our case, those are the standard packages numpy, scipy, matplotlib, decimal and the non-standard package lmfit [38]. The packages are loaded with

⬇

import matplotlib.pyplot

% matplotlib inline

The last two lines are required for graphical visualization directly inside the notebook.

Python support object oriented programing, which means that properties or sub-functions can be call using dotted notation object.property(). For instance, all the Refprop functions will be called using the syntax refprop.function_name().

Arrays can be created explicitly like x=[1,2,3,4,5], or can be defined generically and values added when required:

⬇

y.append(5)

Values are extracted with $\cdot$[i], where i is the position in the array, starting from 0. In our example, x[1]=2 and y[0]=5. The length of the array x is given by len(x), and the array can by printed out explicitly with print(x).

A function $f(x,y,z)=x^{2}+y^{2}-\text{e}^{z}$ can be define by

⬇

function = x**2 + y**2 - numpy.exp(z)

return function

Note that the indentation here is crucial: Python use indentation to structure the code instead of braces like other languages. The same code with no indentation will produce an error.

To shorten this paper, variables assignment will sometimes be written on a single line

On this single line, the value 5 is assign to the variable x, the string ’yes’ is assign to y and the array [1,2,3] is assign to z. Let us write

⬇

y = x**2

The first line creates an array x containing the values from 0 to 49 while the second line creates an array y where each value is the square of the ones in the x array. In Python, there is a distinction between a list and an array, but we will overlook it since it is basically a matter of performance. The term array will then be used for both lists and arrays. Finally, a simple point plot of y as a function of x can be made with

⬇

matplotlib.pyplot.show(fig)

Refprop Python functions usually return results in the form of a dictionary data structure, which associate a key (a string) to a value. Dictionaries are in fact simply arrays where elements are extracted using a key instead of an index. Then, to extract the numerical value of a variable in a dictionary data structure, we need to add the key associated to that variable in brackets after the function. For example, if the function f() return the dictionary

the command f[’P’] will return the numerical value 5.

3.2 Simplest case: Supercritical pure gas

In this section, we implement the simplest case of a supercritical pure gas far from its critical point. The possibles complications will be added successively in the following sections.

3.2.1 Evaluation of the excess adsorption.

First, we need to import the Refprop package, define the installation path and initialize the fluid that we want to work with (methane in our example):

⬇

refprop.setpath(path=’/usr/local/lib/refprop’)

refprop.setup(u’def’, u’METHANE’)

Next, we need to import experimental dataset of bulk phase pressure and excess adsorption. In the following, we will use experimental data from Ref. [23, run 1]. We can use the function numpy.loadtxt() to import data from a text file. Assuming that the first column in data.txt contain the experimental pressure values of the bulk phase (in KPa) and the second column (separated by white spaces) contain the excess gas adsorbed (in mol/Kg), then the code

⬇

dataAd = numpy.loadtxt(’data.txt’, dtype=’float’, usecols=[1])

creates the two arrays dataP and dataAd for pressure and excess adsorption, respectively.

Refprop functions usually required gas density instead of pressure. The gas density can be obtained using refprop.press(T,d,x), with T the temperature (in K), d the density (in mol/L) and x the gas molar fraction array (for pure gases, x=[1], for a 25% / 75% mixture, x=[0.25, 0.75], etc). The gas molar fraction array will be noted xB and x_z for the bulk phase and the adsorbed phase (at position z) respectively. Of course, for pure gases, xB=x_z=x=[1]. It can be noted that we used d for the density $\rho$ and eps0 for $\varepsilon_{0}$ for convenience. To transform the pressure array dataP into a density array dataD, we first define a function pressure that depend on density as

⬇

return refprop.press(T, float(d), x)[’p’] - p

Then, for each value of pressure in dataP[i] ($i=0\dots N-1$, where $N$ is the number of data), we solve the equation pressure(d,T,xB,p)=0 to find the associated density d. We use the function scipy.optimize.fsolve(f,x0,args=(...)) for that purpose (which is based on MINPACK’s hybrd and hybrj algorithms), but any root finding algorithm would work. Here, f is the function to be solved (assuming the form f=0) and x0 is an initial guess for f. If the function f is multidimensional (as in our case), the parameter args=(...) enable the specification of all the variables but the first in f. For example, if we want to solve f(u,v,w)=0 with respect to u, using the initial guess u=1 and for the specific values v0 and w0, we use scipy.optimize.fsolve(f,1,args=(v0,w0)).

Because there is a bijection between pressure and density, the specification of T, xB and p in pressure() function generate a unique root and then, any reasonable guess for the density will work at this stage, so we take d0=1. Explicitly, this give

⬇

for i in range(0, len(dataP)):

d = scipy.optimize.fsolve(pressure, 1, args=(T, xB, dataP[i]))

dataD.append(float(d[0]))

Simply state, this code solve the equation pressure(d,T,xB,p)=0 with respect to d for known values of T, xB and p, given by args=(T,xB,dataP[i]), starting from an initial guess d0=1, and with $i\in(0,\dots,N-1)$. The results are then put in the array dataD in the last line. Note that fsolve() return an array containing numpy.float64 data type here. To get the numerical value only (meaning a standard float data type), we have to use float(d[0]), which take the first element of the array d and ensures a float data type. Python float data type correspond to the usual double precision (binary64) data type [39].

Next, we need to know the maximum density allowed by the Refprop for the gas considered. We will define a function d_max() to lighten the code. This function will take the temperature T and the molar fraction array x of the mixture (x=[1] everywhere for pure gas), pass it to the function refprop.limitx() and extract the maximum density allowed:

⬇

return refprop.limitx(x, htype=’EOS’, t=T, D=0, p=0)[’Dmax’]

The Refprop also have an embedded function for the chemical potential chempot(T,d,x)[’u’]. This function returns an array containing the chemical potential of all the component of a gas mixture with molar fraction array given by x. For simplicity, let us define a function chem as

⬇

return refprop.chempot(T, float(d), x)[’u’]

In the case of pure gas (or the first component of a gas mixture), chem(T,d,x)[0] will return the numerical value of the chemical potential. This notation will be quite useful for mixtures below.

We can also implement the DRA potential Eq. (2) (with $\beta=2$) as

⬇

return eps0*numpy.log(z0/z)**(1.0/2.0)

From the previous functions, the chemical potential of the fluid at a distance z from the surface of the adsorbent can be written as

In this equation, the unknowns variables are z and d_z (the density in the adsorbed phase at z). For the moment, we supposed that the values of z0 and eps0 are known. Note that we will use dB, d_z, xB and x_z to make it clear that we talk about the bulk or adsorbed phase (remember that xB=[1]=x_z for pure gas, but we will leave xB and x_z here for generality). To extract the density d_z, we need to provide a value for z in the range $0\leq z\leq z0$, and solve the following equation

Here again, we use fsolve() for that purpose, along with an initial guess d0 for the density (let say d0=dB). Explicitly, we can build a density function d_pure() as

⬇

def f(d_z):

    y = chem(T, dB, xB)[0] + DRA(z0, eps0, z) - chem(T, d_z, x_z)[0]

    return y

if chem(T, dB, xB)[0] + DRA(Z0, eps0, z) >= chem(T, d_max(xB), xB)[0]:

    return d_max(xB)

else:

    return scipy.optimize.fsolve(f, dB)[0]

This function return the density of the gas for z in the range $0\leq z\leq z0$. On the fifth line, we added a test condition to ensure that the algorithm will not exceed the maximum density allowed by the Refprop. By putting a cut-off on the chemical potential, we in fact assume that the chemical potential is a monotonous increasing function of the density, which means that $\mu_{max}$ correspond to $\mu(\rho_{max})$. It is this cut-off that takes care of the discontinuity in the DRA potential in Eq. (2) as $z\to 0$. The error introduce by this cut-off is negligible over the whole range $0\leq z\leq z_{0}$. In fact, instead of using the maximal density value d_max, the results of the implementation remain unchanged if we use a drastic cut-off value $\rho=0$.

Figure 1 is an example of the density profile given by d_pure() as a function of the ratio z/z0 for pure CH4, N2 and CO2 on activated carbon Calgon Filtrasorb 400 at 318.2 K with a bulk gas density of 0.45 mol/L (corresponding to bulk pressure arround 1.2 MPa). Values of $\varepsilon_{0}$ and $z_{0}$ used are the one given in table 2. As said before, the experimental data are taken from Ref. [23, run 1].

Finally, we can integrate the function d_pure on z from 0 to z0 to obtain the adsorbed gas quantity. We use the function scipy.integrate.quad() for that purpose:

⬇

y = scipy.integrate.quad(d_pure, 0, z0, args=(T, z0, eps0, dB, xB, xB))

return y[0] - dB*z0

The output of this function is the (Gibbs) excess adsorption (in mol/Kg) considering the temperature, the parameters z0 and eps0, the bulk gas density dB (in mol/L) and the molar fraction xB=[1]. Note that in the last line, we subtract the gas quantity that would be present in the absence of adsorbent, namely dB*z0, in order to obtain the excess adsorption and not the absolute adsorption.

3.2.2 Building and fitting the model.

We have shown how the excess adsorption could be evaluated numerically, provided that z0 and eps0 are known. In the following, we will show how to obtain those values by fitting the MPTA model to experimental dataset for pure gases.

First of all, minimization algorithms are built to minimize a function. Then, we have to build a function, based on our calculated excess adsorption, which can be minimized with respect to parameters z0 and eps0. To do that, we can define the function pure_fit where arguments will be pass by a dictionary data structure this time (as required by the chosen minimization algorithm):

⬇

value = params.valuesdict()

z0 = value[’z0’]

eps0 = value[’eps0’]

difference = [ ]

for i in range(0, len(dataD)):

    difference.append(N_ex_pure(T, z0, eps0, dataD[i], xB) - dataAd[i])

return difference

In words, this function accept the argument params which is a dictionary data structure containing z0 and eps0. For each experimental value of density contain in dataD, the function takes the difference between the excess adsorption given by N_ex_pure(), evaluated with the values z0 and eps0, and the experimental adsorption contain in dataAd. Those differences between calculated and experimental adsorption data are put in the array difference which have the same dimension as dataD and dataAd. The upshot is that the function pure_fit() do not return a value but rather an array of values. By minimizing the function pure_fit(), we, in fact, minimize the array just built and thus, the differences between the calculated and experimental excess adsorption. To do that, we used the function minimize() from the package lmfit [38], which use a Levenberg-Marquardt algorithm for non-linear least-squares minimization. Explicitly, we have

⬇

params.add(’z0’, value=1, min=0, vary=True)

params.add(’eps0’, value=10000, min=0, vary=True)

result = lmfit.minimize(pure_fit, params)

print(lmfit.fit_report(result))

The function lmfit.Parameters() create a dictionary named params on the first line. Parameters z0 and eps0 are added to the dictionary on the second and third line, specifying the initial and minimal values, and setting them as adjustable parameters with vary=True. Then, the minimization takes place on the fourth line on the function pure_fit() using parameters params. The fitted values of z0 and eps0 and their uncertainties are given in the fit report in the last line. The package lmfit have the advantage of easily allowing the modification of fitting variables, simply by changing vary=True to vary=False. To give an idea of the computational time involved, the preceding fit took between 2 min (N2 with 11 data points) to nearly 7 min (CO2 with 13 data points) on an Intel i7-2600 CPU running one thread.

Finally, from the fitted values fit_z0 and fit_eps0 (given by lmfit.fit_report()), we can build an array dataFit that will contain the calculated excess adsorption to be compared with experimental dataAd.

⬇

for i in range(0, len(dataD)):

dataFit.append(N_ex_pure(T, fit_z0, fit_eps0, dataD[i], xB))

Figure 2 give an example for pure CH4, N2 and CO2 (dataset taken from Ref. [23]). The initial values of z0 and eps0, the fitted values and their numerical uncertainty are given in Table 2. The mean error between experimental data and the model are given in Table 3. Figure 3 shows the flow chart of the MPTA implementation for pure gases.

3.2.3 Pure gas complications and pitfalls.

The implementation describes until now work but can present some problems and instabilities even in the case of pure gases. For one thing, the experience had shown that finding the root of

in the function d_pure(), can strongly depend on the initial guess d0 for the density. In fact, we remarked that a poor initial guess would be responsible for the majority of instabilities in the code. It is then useful to develop a robust and systematic way of choosing this initial guess.

First, the function N_ex_pure() can be modified in order to use a discrete sum instead of the scipy.integrate.quad integration function. Let us divide the interval (0,z0) in $N$ subintervals of equal length delta. For each subinterval, starting from z0, the function d_pure() is evaluated, but this time, the result of the preceding subinterval is used as initial guess for d0 (of course, for the first subinterval, the bulk density should be a reasonable initial guess). As long as the gas do not go through a phase transition, its density is expected to be a smooth function and then, the density of the gas at z should be close to the one at z+delta. Explicitly, the function N_ex_pure() become

⬇

N = 300

delta = z0 / N

d0 = dB

integral = 0

for i in range(0, N):

  d_z = d_pure(z0 - i*delta - delta/2, T, z0, eps0, dB, xB, x_z, d0)

  integral += d_z*delta

  d0 = d_z

return integral - dB*z0

Note that we add the variable d0 to the function d_pure(), which is pass to fsolve() to be able to specify the initial guess. Also, we do the sum from z0 to 0 and not the other way around because the gas density at z0 is known and given by dB. The choice $N=300$ is arbitrary but we obtain a relatively good precision in a reasonable amount of time with this value.

Now if a phase transition occurs in the adsorbed phase, it can produce a sharp variation in the density such that the smoothness hypothesis no longer holds. In that case, the initial guess for density needs to be adjusted to account for that phase transition. To do that, the density of the gas at dew point and the density of the liquid phase are required. Those density values can be obtained from Refprop. First, we need to ensure that the gas is in the subcritical regime. This can be done by looking at the critical temperature of the gas given by function refprop.info(icomp=1)[’tcrit’]. This function returns the critical temperature of the first gas component (specified with parameter icomp=1). With the gas in subcritical regime, the vapor and liquid density at dew point are obtained with

⬇

refprop._abfl2(’TD’, T, d, xB)[’Dliq’]

In those functions, ’TD’ means that temperature and density are used as input values. Of course, the values of the liquid and vapor density does not depend on the actual fluid density, so we can just put any reasonable number like d=1 for the density. Explicitly, we can define the two parameters d_vap and d_liq as

⬇

d_vap = refprop._abfl2(’TD’, T, 1, xB)[’Dvap’]

d_liq = refprop._abfl2(’TD’, T, 1, xB)[’Dliq’]

else:

d_vap = d_max(xB)

d_liq = d_max(xB)

Note that if the gas is in the supercritical regime, no phase transition can occur. In that case, we assume that the density vary smoothly and there is no need to adjust the initial guessed density. This is why in that case, d_vap and d_liq are both set to d_max.

Now, the function d_pure() can be modified to account for those new parameters. To do that, a condition is added to check if the chemical potential at z is greater than the one for the gas with critical density d_vap. If so, and if the initial guess d0 is less than the liquid density d_liq, we will set d0 just over the liquid density (namely 1.1*d_liq) to help the solver fsolve(). Here again, there is the underlying assumption that the chemical potential is a monotonous increasing function of the density such that $\mu(\rho_{liq})\geq\mu(\rho_{vap})$ because $\rho_{liq}\geq\rho_{vap}$. The d_pure() function now become

⬇

def f(d_z):

    y = chem(T, dB, xB)[0] + DRA(z0, eps0, z) - chem(T, d_z, x_z)[0]

    return y

if chem(T, dB, xB)[0] + DRA(z0, eps0, z) >= chem(T, d_max(xB), xB)[0]:

    return d_max(xB)

if chem(T, dB, xB)[0] + DRA(z0, eps0, z) > chem(T, d_vap, xB)[0] and d0 < d_liq:

    d0 = 1.1*d_liq

return scipy.optimize.fsolve(f, d0)[0]

Those modifications to d_pure() and N_ex_pure() functions should give a more robust algorithm in case of instabilities and phases transitions. For example, Figure 4 shows the fit of the model for pure subcritical CO2 on Norit R1 activated carbon at 298 K, considering original and modified d_pure() and N_ex_pure() functions (experimental dataset taken from Ref. [40]). The modified d_pure() function exhibit a phase transition of the CO2 from gas to liquid. The fluid phase can be obtain from Refprop using

⬇

print(refprop.getphase(fluid_dict))

3.3 Gas mixtures

Numerical implementation of gas mixtures is a bit more complicated. The MPTA model assumes a common value for $z_{0}$ to limit the number of adjustable parameters, as discussed by Bjørner et al. [41]. This implies that the model must be simultaneously minimized for all pure components isotherms constituting the mixture. This can be achieved readily by building an array, as in function pure_fit(), but now containing all the gas components. For example, in a mixture of two gases A and B, the difference[i], $0\leq i<N_{A}$ will contain the first component, exactly like in the pure case, while difference[i] with $N_{A}\leq i<N_{A}+N_{B}$ will contain the second gas B. For simplicity, only binary mixtures will be considered, but generalization to more complex mixtures is straightforward.

It is useful to define the molar fraction of pure gases and mixture in separate arrays. Let us consider a bulk binary mixture of $80\%$ molar fraction for the first gas and $20\%$ molar fraction for the second. We can define the two arrays

⬇

X = [[1, 0] , [0, 1]]

Here, xB is the bulk molar fraction array while X[i] represent the molar fraction array of pure gas A when $i=0$, and pure gas B when $i=1$. For a ternary mixture, X will be given by X=[[1,0,0],[0,1,0],[0,0,1]] and so on.

In the following, we consider a binary mixture composed of $80\%$ CH4 and $20\%$ CO2. The Refprop is initialized accordingly with

With those definitions, we can redefine d_vap and d_liq as

⬇

for i in range(0, len(xB)):

if T < refprop.info(icomp=i+1)[’tcrit’]:

    d_vap.append(refprop._abfl2(’TD’, T, 1, X[i])[’Dvap’])

    d_liq.append(refprop._abfl2(’TD’, T, 1, X[i])[’Dliq’])

else:

    d_vap.append(d_max(X[i]))

    d_liq.append(d_max(X[i]))

Then, d_vap[0] and d_liq[0] will returned the values for gas A and d_vap[1] and d_liq[1] will returned the values for gas B in the binary mixture.

As explain in Section 2.2, in order to calculate the gas density in the adsorbed phase, a system of equations need to be solved in which the gas molar fraction is one of the variables. A possible problem related to this stems from the fact that when the solver tries to adjust this molar fraction, it can end up with floating point numbers that do not have an exact binary representation. This is known as the representation error and then, the program rounds the number and introduce a small error which is usually insignificant. However, Refprop functions have a consistency test that checks if the sum of the molar fractions ends up to precisely 1. Then, the ${\sim 10^{-16}}$ error introduced by the rounding can make the consistency test to fail and return errors. One way to circumvent this is to use the package decimal that provide arbitrary precision floating point arithmetic.

Moreover, the argument used to justify the automation of initial guess d0 for the density now also holds for an initial guess for the gas mixture molar fraction x0, as long as the molar fraction vary smoothly. Then, for gas mixtures, the density function now need to receive initial guesses for density and gas molar fractions. Explicitly, the function takes the form

⬇

decimal.getcontexr().prec = 100

def f(u):

    d_z = float(u[0])

    x_z_A = decimal.Decimal(u[1])

    x_z_B = decimal.Decimal(1 - x_z_A)

    x_z = [x_z_A , x_z_B]

    y1 = chem(T, dB, xB)[0] + DRA(z0, eps0[0], z) - chem(T, d_z, x_z)[0]

    y2 = chem(T, dB, xB)[1] + DRA(z0, eps0[1], z) - chem(T, d_z, x_z)[1]

    return [y1 , y2]

if chem(T, dB, xB)[0] + DRA(z0, eps0[0], z) >= chem(T, d_max(X[0]), xB)[0] or \

    chem(T, dB, xB)[1] + DRA(z0, eps0[1], z) >= chem(T, d_max(X[1]), xB)[1]:

    return [d_max(xB) , x0[0]]

if chem(T, dB, xB)[0] + DRA(z0, eps0[0], z) > chem(T, d_vap[0], xB)[0] and \

    d0 < d_liq[0]:

    d0 = 1.1 * d_liq[0]

if chem(T, dB, xB)[1] + DRA(z0, eps0[1], z) > chem(T, d_vap[1], xB)[1] and \

    d0 < d_liq[1]:

    d0 = 1.1 * d_liq[1]

return = scipy.optimize.fsolve(f, [d0 , x0[0]])

The second line set the precision for the Decimal() function. With this precision, the representation error is $\sim 10^{-51}$ which is small enough for the consistency test of the Refprop functions to succeed. Most of the properties and parameters referring to a pure gas are now described by arrays of two values, one for each gas component. Then, chem(...)[i], eps0[i], d_vap[i] and d_liq[i] refer to the component A for i=0 and component B for i=1. The variables dB, xB refer to the bulk density and molar fractions while d_z and x_z are the unknowns and refer to the density and the molar fraction of the adsorbed phase at distance z from the adsorbent surface. The function fsolve() now have to solve a system of two equations. To proceed, we define a function f(u), with u an array containing the gas density d_z and the molar fraction of the first component x_z_A to be found. This function returns an array containing the two equations to be solved. The solver now requires initial guesses in an array format to match the system of equation, which is given by [d0,x0[0]]. Remark that for a mixture of two gases, we only used the molar fraction of the first component and get the molar fraction of the second component by 1 - x0[0]. In the function f(u), the Decimal() function is used to redefine the molar fraction to x_z_A and x_z_B such that 1-x_z_A-x_z_B$\sim 10^{-51}$. Finally, the function d_mix() return an array composed of the gas density and the molar fraction of the first component.

This d_mix() function induced modifications to the N_ex_pure function that becomes

⬇

N  = 300

delta = z0 / N

d0 = dB

x0 = xB

integral_A , integral_B = 0 , 0

for i in range(0, N):

    [d_z, x_z_A] = d_mix(z0 - i*delta - delta/2, T, z0, eps0, dB, xB, d0, x0)

    integral_A += d_z*x_z_A*delta

    integral_B += d_z*(1 - x_z_A)*delta

    d0 = d_z

    x0 = [x_z_A , 1 - x_z_A]

return [integral_A - dB*xB[0]*z0 , integral_B - dB*xB[1]*z0]

In this function, d_mix() return the two variables d_z and x_z_A. Variables integral_A and integral_B represent the excess adsorption for each component in the mixture. As we can see, the terms d_zx_z_Adelta and d_z*(1-x_z_A)delta give the absolute adsorption for component A and B respectively, in the interval delta. The initial guesses d0 and x0, for density and molar fraction, are updated with newly calculated values d_z and [x_z_A ,1-x_z_A]. Finally, the (Gibbs) excess adsorption for each component are given by the absolute adsorption minus the term dBxB[0]z0 or dBxB[1]*z0.

The mix_fit() function for mixtures take the form

⬇

value = params.valuesdict()

z0 = value[’z0’]

eps0_A = value[’eps0_A’]

eps0_B = value[’eps0_B’]

difference = [ ]

for i in range(0, len(dataD_A)):

    difference.append(N_ex_pure(T, z0, eps0_A, dataD_A[i], X[0]) - dataAd_A[i])

for i in range(0, len(dataD_B)):

    difference.append(N_ex_pure(T, z0, eps0_B, dataD_B[i], X[1]) - dataAd_B[i])

return difference

Of course, N_ex_pure() function is used here and not N_ex_mix() because the purpose of mix_fit() is to optimized the fit of pure gas A and B simultaneously with an unique z0. Finally, the fit is executed with

⬇

params.add(’z0’, value=1, min=0, vary=True)

params.add(’eps0_A’, value=10000, min=0, vary=True)

params.add(’eps0_B’, value=10000, min=0, vary=True)

result = lmfit.minimize(mix_fit, params)

print(lmfit.fit_report(result))

For a CH4 / CO2 mixture on Calgon F-400 at 318.2K [23], our fitted parameters are given in Table 4. Figure 5 represents the fit for pure gases with a common z0. The mean error between experimental data pure isotherms and the fitted MPTA model are $4.62\%~{}/~{}4.51\%$ respectively.

Figure 6, presents the molar fraction of each gas component near the surface of the adsorbent. Finally, in Figure 7, we present the (Gibbs) excess adsorption for a $80\%$ CH4 / $20\%$ CO2 mixture. The figure shows the excess adsorption of each component and the total excess adsorption for the mixture. The error bars represent the expected experimental uncertainties on the measurements given in Ref. [23]. The mean error between experimental data and the model are presented in table 5. It is worth noting that the mean error of the MPTA model is of the same order as the experimental uncertainty on datasets.

To give an idea of the importance of bulk phase molar fraction values at equilibrium, we can recalculate the excess adsorption of this particular mixture using only the feeding molar composition, which differs from the measured molar fraction by at most $10\%$. The resulting mean error between experimental and calculated values increases from $4.42\%$, $4.41\%$ and $3.67\%$ to $25\%$, $71\%$ and $4.15\%$ for the CH4 component, the CO2 component and the mixture, respectively.

4 Conclusion

We presented a detailed numerical implementation of the MPTA adsorption model, covering pure gases and binaries mixtures. From our treatment of the binary mixture case, generalization to more complex mixtures is straightforward. This implementation includes some of the strategies used to obtain a robust model even in case of fast variation of fluid densities, like phases transitions. For all the mixtures in Ref. [23], the mean error between predictions of the MPTA model for mixture (Gibbs) excess adsorption and experimental data vary from $1.44\%$ ($20\%$ CH${}_{4}\,$/$\,80\%$ CO2 mixture) to $5.09\%$ ($20\%$ N${}_{2}\,$/$\,80\%$ CO2 mixture). Meanwhile, the experimental uncertainties for mixture adsorption vary from $2.4\%$ ($80\%$ CH${}_{4}\,$/$\,20\%$ N2 and $40\%$ CH${}_{4}\,$/$\,60\%$ N2 mixtures) to $4.1\%$ ($80\%$ N${}_{2}\,$/$\,20\%$ CO2 mixture).

For a complete version of the code, the interested reader can consult:

https://github.com/RaphaelGervaisLavoie/MPTA

Acknowledgements

The authors would like to thanks Ege Dundar for his help and suggestion about the implementation and also the Natural Sciences and Engineering Research Council of Canada and the Savannah River National Laboratory for financial support.