Numerical implementation of the multicomponent potential theory of adsorption in Python using the NIST Refprop database
Raphaël Gervais Lavoie, Mathieu Ouellet, Jean Hamelin, Pierre Bénard
aa r X i v : . [ c ond - m a t . o t h e r] F e b Numerical implementation of the multicomponentpotential theory of adsorption in Python using the NISTRefprop database
Rapha¨el Gervais Lavoie ∗ , Mathieu Ouellet, Jean Hamelin, and Pierre B´enardInstitut de recherche sur l’hydrog`ene, Universit´e du Qu´ebec `a Trois-Rivi`eres,D´epartement de chimie, biochimie et physique, Universit´e du Qu´ebec `aTrois-Rivi`eres. Abstract
In this paper, we present a detailed numerical implementation of the multicomponentpotential theory of adsorption which is among the most accurate gas mixtures adsorptionmodels. The implementation uses the NIST Refprop database to describe fluid propertiesand applies to pure gases and mixtures in both subcritical and supercritical regimes. Thelimitations of the model and the issues encountered with its implementation are discussed.The adsorption isotherms of CH / CO mixture are modeled and parameterized asimplementation examples. Keywords: adsorption; mixture adsorption; multicomponent adsorption; potential theoryof adsorption; MPTA; density functional theory
Gas adsorption in porous material is relevant to a wide variety of industrial and scientificapplications, ranging from gas purification and separation to gas storage in stationary ormobile applications. While adsorption experiments of pure gases can be performed readily,experiments on gas mixtures can be challenging, expensive and time-consuming [1]. Thissituation motivates the utilization of theoretical models to predict the behavior of gas mixturesin the presence of adsorbent.Several approaches have been proposed to model the adsorption isotherms of both puregases and mixtures [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. In this paper, we willconcentrate on the potential theory of adsorption (PTA), a two-parameter thermodynamicmodel developed by Shapiro and Stenby [16], based on the pore filling approach of Polanyi’stheory of adsorption [17]. For gas mixtures adsorption, the PTA model generalized into themulticomponent potential theory of adsorption (MPTA), which can deliver great performance,especially when considering disordered adsorbents such as activated carbon (an overview ofmixture adsorption models is presented in Ref. [18]). For a N components gas mixture, thesimplest MPTA model required the adjustment of only N + 1 parameters, obtained from purecomponents adsorption isotherms.The objective of this paper is to present a detailed numerical implementation of the MPTAmodel. In this spirit, we will write down explicit lines of code, written in Python3, an ∗ Email: [email protected] µ B , chem() function Chemical potential of the bulk phase (J/mol) µ Ad , chem() function Chemical potential of the adsorbed phase (J/mol) ε , DRA() function Adsorbent surface potential (J/mol) ρ B , db float Fluid density in the bulk phase (mol/L) ρ Ad , d_z float Fluid density in the adsorbed phase (mol/L) ε , eps0 float, array Characteristic energy of adsorption (J/mol) z , z float Microporous volume (cm /g) z , z0 float Limiting micropore volume (cm /g) β 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) 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 wherethe 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 significantlyaffected by the adsorbent material. The thermodynamic properties of the bulk phase aredescribed by the equation of state (EoS) of the fluid.
The main assumption of the MPTA model is to suppose that the fluid–surface interaction iscompletely described by a local potential field ε , generated by the surface. The potential ε is attractive near the surface of the adsorbent and tends to vanish far from the surface. Theconstitutive equation of the MPTA model is given by Refs. [15, 16] µ B ( T, ρ B ) = µ Ad ( T, ρ Ad ) − ε, (1)2here µ B and ρ B are the chemical potential and the fluid density in the bulk phase, while µ Ad and ρ Ad are the locals chemical potential and fluid density in the adsorbed phase. The bulkphase properties are assumed to be constant while the adsorbed phase properties vary withposition [16]. Hence, the density ρ Ad is usually very high inside pores, close to the adsorbent,and decrease to reach ρ B far from the adsorbent. Given Eq. (1), the local thermodynamicproperties in the adsorbed phase are uniquely determined from properties of the bulk phaseand the interaction potential ε . Since we consider only equilibrium situations, the temperatureis constant.Along with Eq. (1), the additional information needed is an explicit form for the potential ε . We use the Dubinin–Radushkevich–Astakhov (DRA) equation [19], a two-parameter ( ε and z ), semi-empirical potential originates from the theory of volume filling micropores [20, 21]and which can represent a wide range of energy distributions ε ( z ) = ε (cid:16) ln z z (cid:17) β . (2)Here, ε is the characteristic energy of adsorption of the adsorbate on the adsorbent, aparameter representing the magnitude of the fluid-solid interaction. The variable z has unitsof a volume and z is the limiting micropore volume of the adsorbent which is the volume ofthe pore that significantly contributes to the adsorption. The ratio z/z represents the fractionof the microporous volume associated to an energy ε ( z ). β is a measure of the heterogeneity ofthe adsorbent and is usually set to 2 for activated carbon, but can eventually become a fittingparameter (in this paper, we will use β = 2). Under certain conditions, the variable z can beseen as the distance from the surface of the adsorbent. This is the case when the microporousstructure 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 fromthe surface of the adsorbent as in Ref. [16].The presence of the divergence in Eq. (2) as z → ρ Ad ( z ) from the chemicalpotentials. Then, the (Gibbs) excess adsorption N ex (which is what’s experimentallymeasured), is calculated by N ex ( ρ B ) = Z z ( ρ Ad ( z ) − ρ B ) dz. (3)We should emphasize a crucial point here: even if ε and z have well-defined physicalmeaning, these two parameters are better used as fitting parameters in the model. Values forthe characteristic energy of adsorption ( ε ) for physisorption on microporous adsorbents rangefrom 5–15 kJ/mol while typical limiting micropore volume ( z ) range from 0.1 to 2 cm /g.To fit the model to experimental datasets, initial guess of ε and z are made. Thedifferences between the calculated excess adsorption ( N ex ( ρ B )) and the experimental data atvarious pressures are then simultaneously minimized. The result of the minimization betweencalculated and experimental values of excess adsorption translates into optimal values for ε and z .In the literature, many authors used the fugacity f instead of the chemical potential. Inthat case, Eq. (1) takes the form f Ad ( T, ρ Ad ( z )) = f B ( T, ρ B ) exp (cid:18) ε ( z ) RT (cid:19) . (4)Fugacity being essentially the exponential of the chemical potential, we choose to work with thelater. Also, it is worth noting that when working with the chemical potential, the logarithmicsingularity of Eq. (2) is integrable, which is not the case when working with fugacity.3 .2 Mixtures For gas mixtures, the situation is more challenging, particularly for the description of the localchemical potential in the adsorbed phase. The molar fraction of the component i will be noted x i such that P Mi x i = 1 where M is the number of pure components in the mixture. Whilethe molar fraction of the bulk phase can be well known from experimental measurements (gaschromatography for instance), the molar fraction in the adsorbed phase is unknown and mayvary substantially from the bulk phase. In gas mixtures, each gas component interacts withthe surface through is own potential ε i ( z ) = ε i (cid:16) ln z z (cid:17) β , (5)where i = 1 . . . M . Due to the selectivity of adsorbent material, the adsorption of onecomponent will be favored, affecting the local molar fraction of the mixture in the adsorbedphase.The introduction of the molar fractions adds unknown variables to the problem. Instead ofsimply inverting Eq. (1) to find ρ Ad ( z ), we need to solve simultaneously the following systemof equations [22] µ iB ( ρ B , x iB ) + ε i ( z ) − µ iAd ( ρ Ad ( z ) , x iAd ( z )) = 0 , (6)for x iAd ( z ) and ρ Ad ( z ), where the superscript i refer to the gas component of the mixture. Fromnow on, the chemical potential also depend on molar fraction x i of gas components ( x iB is themolar fraction of the component i of the bulk phase at equilibrium, i.e. when the adsorptionprocess is complete). In Eq. (6), ρ B , x iB and ε i ( z ) are known and ρ Ad ( z ) and x iAd ( z ) are theunknowns. The dimension of the system of equation to solve is M . The set of x iAd ( z ) represent M − P Mi x iAd = 1 decrease by 1 the number of independentvariables) and there is one more independent variable ρ Ad ( z ) which, of course, is the same forall gas components.Once the x iAd ( z ) and ρ Ad ( z ) are determined, the next step is to calculate the excessadsorption for each component: N iex ( ρ B ) = Z z (cid:0) ρ Ad ( z ) x iAd ( z ) − ρ B x iB (cid:1) dz. (7)Finally, the total amount of excess gas adsorbed is simply given by N ex ( ρ B ) = M X i =1 N iex ( ρ B ) . (8)As for the pure gases, to adjust the model to experience, initial values of ε i and a common z must be postulated. Because z is considered to be an adsorbent property (the limitingmicropore volume), we use a single value for all gas components. The excess adsorption iscalculated for each gas component from the postulated initial values. Then, the differencesbetween the calculated and experimental pure excess adsorption is minimized simultaneouslyfor all experimental bulk phase densities and all gas components. This finally gives the optimalvalues of parameters ε i and a common z .A key point here is that even for mixtures, the adjustment of the model to experimentaldata is made using pure gas adsorption isotherms , not mixture data [22]. Once the parametersare adjusted from pure gas isotherms, then gas mixture adsorption can be simulated. Itis important, however, to understand that x B should be the bulk phase molar fraction atequilibrium which may differ greatly from the initial or feeding molar fraction . For example, in4ef. [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 themodel predictions. Consequently, one should be careful when using the MPTA model on datawith unknown bulk phase equilibrium molar fractions. Up to now, there is one thing that we left aside: the way we will evaluate the chemicalpotentials µ iB and µ iAd . In order to accurately describe the thermodynamic properties of thefluid, one should use an accurate equation of state (EoS). A well-known EoS commonly usedfor this kind of purpose is the Soave–Redlich–Kwong (SRK) EoS [24, 25]. However, in thispaper, we choose to use the NIST–REFPROP (Refprop) database [26] which give presumablythe most accurate results. For the purpose of this paper, the underlying references used bythe Refprop are Refs. [27, 28, 29]. Credit goes to Refs. [30, 31] for the Python integration ofRefprop. To quantify the precision of the model, we will use the absolute average deviations (AAD)between experimental and fitted data for excess adsorption:
ADD ( N ) = 100 k k X j =1 (cid:12)(cid:12)(cid:12)(cid:12) N expex ( ρ jB ) − N calex ( ρ jB ) N expex ( ρ jB ) (cid:12)(cid:12)(cid:12)(cid:12) % . (9)In this equation, ρ jB is one of the experimental bulk phase density in the set of k values { ρ B , ρ B , . . . , ρ kB } . N expex ( ρ jB ) and N calex ( ρ jB ) are the experimental and calculated gas adsorbeddensity at bulk phase density ρ jB , respectively. The proposed numerical integration being not very computationally expensive, Python (version3) [32, 33] might be the best choice of language to obtain a simple, concise and open-sourceimplementation. Indeed, with all the embedded functions and algorithms (principally fromSciPy [34] and NumPy [35] packages), most of the numerical operations are easier to doin Python than in other languages. Moreover, the Jupyter (IPython) notebook [36] andmatplotlib [37] package provide great and convenient tools for graphical visualization andanalysis.For readers unfamiliar with Python, the language is simple enough to serve as pseudocodefor other languages as well. Also, Section 3.1 give a short introduction to Python. Notethat some programming convention will be overlooked for space consideration and that theproposed implementation is oriented on simplicity and concision rather than performance.
We present here some of Python concepts needed for the rest of this paper. Firstly, the requiredpackage 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 numpy, scipy, lmfit, decimal import matplotlib.pyplot% matplotlib inline
The last two lines are required for graphical visualization directly inside the notebook.5ython support object oriented programing, which means that properties or sub-functionscan be call using dotted notation object.property() . For instance, all the Refprop functionswill 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 valuesadded when required: y = [ ]y.append(5)
Values are extracted with · [i] , where i is the position in the array, starting from 0. In ourexample, x[1]=2 and y[0]=5 . The length of the array x is given by len(x) , and the array canby printed out explicitly with print (x) .A function f ( x, y, z ) = x + y − e z can be define by def f(x, y, z):function = x**2 + y**2 - numpy.exp(z) return function Note that the indentation here is crucial: Python use indentation to structure the code insteadof 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 x, y, z = 5, ’yes’, [1, 2, 3]
On this single line, the value 5 is assign to the variable x , the string ’yes’ is assign to y andthe array [1,2,3] is assign to z . Let us write x = numpy.arange(0, 50)y = x**2 The first line creates an array x containing the values from 0 to 49 while the second line createsan array y where each value is the square of the ones in the x array. In Python, there is adistinction between a list and an array, but we will overlook it since it is basically a matterof performance. The term array will then be used for both lists and arrays. Finally, a simplepoint plot of y as a function of x can be made with fig = matplotlib.pyplot.plot(x, y)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 whereelements are extracted using a key instead of an index. Then, to extract the numerical valueof a variable in a dictionary data structure, we need to add the key associated to that variablein brackets after the function. For example, if the function f() return the dictionary f = {’P’: 5, ’T’: 298, ’SubCrit’: True}, the command f[’P’] will return the numerical value 5.
In this section, we implement the simplest case of a supercritical pure gas far from its criticalpoint. The possibles complications will be added successively in the following sections.
First, we need to import the Refprop package, define the installation path and initialize thefluid that we want to work with (methane in our example):6 mport refproprefprop.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 dataP = numpy.loadtxt(’data.txt’, dtype=’float’, usecols=[0])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 canbe obtained using refprop.press(T,d,x) , with T the temperature (in K), d the density (inmol/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 phaseand 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 ρ and eps0 for ε for convenience. To transformthe pressure array dataP into a density array dataD , we first define a function pressure thatdepend on density as def pressure(d, T, x, p): return refprop.press(T, float(d), x)[’p’] - p Then, for each value of pressure in dataP[i] ( i = 0 . . . N −
1, where N is the number ofdata), we solve the equation pressure(d,T,xB,p)=0 to find the associated density d . We usethe function scipy.optimize.fsolve(f,x0,args=(...)) for that purpose (which is based onMINPACK’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 . Ifthe function f is multidimensional (as in our case), the parameter args=(...) enable thespecification 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 densitywill work at this stage, so we take d0=1 . Explicitly, this give dataD = [ ] 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 knownvalues of T , xB and p , given by args=(T,xB,dataP[i]) , starting from an initial guess d0=1 , andwith i ∈ (0 , . . . , N − dataD in the last line. Notethat fsolve() return an array containing numpy.float64 data type here. To get the numericalvalue only (meaning a standard float data type), we have to use float(d[0]) , which take thefirst element of the array d and ensures a float data type. Python float data type correspondto 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 thefunction refprop.limitx() and extract the maximum density allowed:7 ef d_max(x): 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 gasmixture with molar fraction array given by x . For simplicity, let us define a function chem as def chem(T, d, x): 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 returnthe numerical value of the chemical potential. This notation will be quite useful for mixturesbelow.We can also implement the DRA potential Eq. (2) (with β = 2) as def DRA(z0, eps0, z): 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 thesurface of the adsorbent can be written as chem(T, d_z, x_z)[0] = chem(T, dB, xB)[0] + DRA(z0, eps0, z) 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 wewill 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). Toextract the density d_z , we need to provide a value for z in the range 0 ≤ z ≤ z
0, and solvethe 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 d_pure(z, T, z0, eps0, dB, xB, x_z): def f(d_z):y = chem(T, dB, xB)[0] + DRA(z0, eps0, z) - chem(T, d_z, x_z)[0] return y i f 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 ≤ z ≤ z
0. On the fifth line,we added a test condition to ensure that the algorithm will not exceed the maximum densityallowed by the Refprop. By putting a cut-off on the chemical potential, we in fact assume thatthe chemical potential is a monotonous increasing function of the density, which means that µ max correspond to µ ( ρ max ). It is this cut-off that takes care of the discontinuity in the DRApotential in Eq. (2) as z →
0. The error introduce by this cut-off is negligible over the wholerange 0 ≤ z ≤ z . In fact, instead of using the maximal density value d_max , the results of theimplementation remain unchanged if we use a drastic cut-off value ρ = 0.Figure 1 is an example of the density profile given by d_pure() as a function of the ratio z/z0 for pure CH , N and CO on activated carbon Calgon Filtrasorb 400 at 318.2 K with abulk gas density of 0.45 mol/L (corresponding to bulk pressure arround 1.2 MPa). Values of ε and z used are the one given in table 2. As said before, the experimental data are takenfrom Ref. [23, run 1]. 8 .0 0.2 0.4 0.6 0.8 1.0Ratio z/z0051015202530 A d s o r b e d p h a s e d e n s i t y ( m o l / L ) ρ B = 0 . mol/L, ǫ ,CH = 8568 J/mol, z ,CH = 0 . cm /g, ǫ ,N = 7020 J/mol, z ,N = 0 . cm /g, ǫ ,CO = 7541 J/mol, z ,CO = 0 . cm /g. CH4N2CO2
Figure 1: Density profile of pure CH , N and CO near the surface of Calgon F-400 at 318.2K.Finally, we can integrate the function d_pure on z from to z0 to obtain the adsorbed gasquantity. We use the function scipy.integrate.quad() for that purpose: def N_ex_pure(T, z0, eps0, dB, xB):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 thetemperature, the parameters z0 and eps0 , the bulk gas density dB (in mol/L) and the molarfraction xB=[1] . Note that in the last line, we subtract the gas quantity that would be presentin the absence of adsorbent, namely dB*z0 , in order to obtain the excess adsorption and notthe absolute adsorption. 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 theMPTA model to experimental dataset for pure gases.First of all, minimization algorithms are built to minimize a function. Then, we have tobuild a function, based on our calculated excess adsorption, which can be minimized withrespect to parameters z0 and eps0 . To do that, we can define the function pure_fit wherearguments will be pass by a dictionary data structure this time (as required by the chosenminimization algorithm): 9 ef pure_fit(params):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 structurecontaining z0 and eps0 . For each experimental value of density contain in dataD , the functiontakes the difference between the excess adsorption given by N_ex_pure() , evaluated with thevalues z0 and eps0 , and the experimental adsorption contain in dataAd . Those differencesbetween calculated and experimental adsorption data are put in the array difference whichhave 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 calculatedand experimental excess adsorption. To do that, we used the function minimize() from thepackage lmfit [38], which use a Levenberg-Marquardt algorithm for non-linear least-squaresminimization. Explicitly, we have params = lmfit.Parameters()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. Param-eters z0 and eps0 are added to the dictionary on the second and third line, specifying theinitial 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 thelast line. The package lmfit have the advantage of easily allowing the modification of fittingvariables, simply by changing vary=True to vary=False . To give an idea of the computationaltime involved, the preceding fit took between 2 min (N with 11 data points) to nearly 7 min(CO 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 canbuild an array dataFit that will contain the calculated excess adsorption to be compared withexperimental dataAd . dataFit = [ ] for i in range(0, len(dataD)):dataFit.append(N_ex_pure(T, fit_z0, fit_eps0, dataD[i], xB)) Table 2: Initial and fitted values of the parameters for pure gases.Parameter Initial value fitted value fit error % z ,CH ( cm /g ) 1 0.269 0 . z ,N ( cm /g ) 1 0.256 0 . z ,CO ( cm /g ) 1 0.306 2 . ε ,CH ( J/mol ) 10 000 8568 0 . ε ,N ( J/mol ) 10 000 7020 0 . ε ,CO ( J/mol ) 10 000 7541 3 . , N and CO (dataset taken from Ref. [23]). Theinitial values of z0 and eps0 , the fitted values and their numerical uncertainty are given inTable 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.Table 3: Precision of the fits for pure gases on Calgon F-400Gas Mean error (ADDn) Experimental uncertainty [23] CH N CO E x e ss a d s o r p t i o n ( m o l / k g ) Error: CH4 0.36 % N2 0.59 % CO2 4.71 % Mean 1.89 %
CH4N2CO2
Figure 2: Result of the fit for pure CH , N and CO on Calgon F-400at 318.2K for individual values of z and ε . The implementation describes until now work but can present some problems and instabilitieseven in the case of pure gases. For one thing, the experience had shown that finding the rootof 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 instabilitiesin the code. It is then useful to develop a robust and systematic way of choosing this initialguess. 11ure gas adsorption data:dataP, dataAdEquation of state:Refprop µρ Max z Initial guesses: z , ε Surface potential:DRAConvert P to ρ : dB Calculate ρ ( z ): d_pure() Calculate excess adsorption:
N_ex_pure(dB)
Compare
N_ex_pure(dB) with experiment: pure_fit()
Minimize the difference:
N_ex_pure - dataAd
New z , ε Final resultFigure 3: Flow chart of the MPTA implementation for pure gases.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 subintervalsof 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 thegas do not go through a phase transition, its density is expected to be a smooth function andthen, the density of the gas at z should be close to the one at z+delta . Explicitly, the function N_ex_pure() become def
N_ex_pure(T, z0, eps0, dB, xB):N = 300delta = z0 / Nd0 = dBintegral = 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*deltad0 = 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 ableto specify the initial guess. Also, we do the sum from z0 to 0 and not the other way aroundbecause the gas density at z0 is known and given by dB . The choice N = 300 is arbitrary butwe 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 variationin the density such that the smoothness hypothesis no longer holds. In that case, the initialguess for density needs to be adjusted to account for that phase transition. To do that, thedensity of the gas at dew point and the density of the liquid phase are required. Those density12alues can be obtained from Refprop. First, we need to ensure that the gas is in the subcriticalregime. 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 firstgas component (specified with parameter icomp=1 ). With the gas in subcritical regime, thevapor and liquid density at dew point are obtained with refprop._abfl2(’TD’, T, d, xB)[’Dvap’]refprop._abfl2(’TD’, T, d, xB)[’Dliq’] In those functions, ’TD’ means that temperature and density are used as input values. Ofcourse, 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 definethe two parameters d_vap and d_liq as i f T < refprop.info(icomp=1)[’tcrit’]: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 guesseddensity. 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 dothat, a condition is added to check if the chemical potential at z is greater than the one forthe gas with critical density d_vap . If so, and if the initial guess d0 is less than the liquiddensity d_liq , we will set d0 just over the liquid density (namely ) to help the solver fsolve() . Here again, there is the underlying assumption that the chemical potential is amonotonous increasing function of the density such that µ ( ρ liq ) ≥ µ ( ρ vap ) because ρ liq ≥ ρ vap .The d_pure() function now become def d_pure(z, T, z0, eps0, dB, xB, x_z, d0): def f(d_z):y = chem(T, dB, xB)[0] + DRA(z0, eps0, z) - chem(T, d_z, x_z)[0] return y i f chem(T, dB, xB)[0] + DRA(z0, eps0, z) >= chem(T, d_max(xB), xB)[0]: return d_max(xB) i f 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 robustalgorithm in case of instabilities and phases transitions. For example, Figure 4 shows the fit ofthe model for pure subcritical CO on Norit R1 activated carbon at 298 K, considering originaland 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 CO from gas to liquid. Thefluid phase can be obtain from Refprop using fluid_dict = refprop._abfl2(’TD’, T, d, x) print (refprop.getphase(fluid_dict)) E x e ss a d s o r p t i o n ( m o l / k g ) Fit for CO2 with unmodified d_pure functionFit for CO2 with modified d_pure function
Figure 4: Result of the fit on experimental dataset for pure CO on Norit-R1 at 298Kwith and without modifications to d pure() and N ex pure() functions. Numerical implementation of gas mixtures is a bit more complicated. The MPTA modelassumes a common value for z to limit the number of adjustable parameters, as discussed byBjørner et al. [41]. This implies that the model must be simultaneously minimized for all purecomponents isotherms constituting the mixture. This can be achieved readily by building anarray, as in function pure_fit() , but now containing all the gas components. For example, in amixture of two gases A and B , the difference[i] , 0 ≤ i < N A will contain the first component,exactly like in the pure case, while difference[i] with N A ≤ i < N A + N B will contain thesecond gas B . For simplicity, only binary mixtures will be considered, but generalization tomore complex mixtures is straightforward.It is useful to define the molar fraction of pure gases and mixture in separate arrays. Let usconsider a bulk binary mixture of 80% molar fraction for the first gas and 20% molar fractionfor the second. We can define the two arrays xB = [0.80 , 0.20]X = [[1, 0] , [0, 1]] Here, xB is the bulk molar fraction array while X[i] represent the molar fraction array of puregas 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% CH and 20% CO . TheRefprop is initialized accordingly with refprop.setup(u’def’, u’METHANE’, u’CO2’) With those definitions, we can redefine d_vap and d_liq as14 _vap , d_liq = [ ] , [ ] for i in range(0, len(xB)): i f 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] willreturned 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, asystem 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 thismolar fraction, it can end up with floating point numbers that do not have an exact binaryrepresentation. This is known as the representation error and then, the program rounds thenumber and introduce a small error which is usually insignificant. However, Refprop functionshave a consistency test that checks if the sum of the molar fractions ends up to precisely 1.Then, the ∼ − error introduced by the rounding can make the consistency test to fail andreturn errors. One way to circumvent this is to use the package decimal that provide arbitraryprecision floating point arithmetic.Moreover, the argument used to justify the automation of initial guess d0 for the densitynow also holds for an initial guess for the gas mixture molar fraction x0 , as long as the molarfraction vary smoothly. Then, for gas mixtures, the density function now need to receive initialguesses for density and gas molar fractions. Explicitly, the function takes the form def d_mix(z, T, z0, eps0, dB, xB, d0, x0):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] i f 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]] i f 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] i f 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, therepresentation error is ∼ − which is small enough for the consistency test of the Refpropfunctions to succeed. Most of the properties and parameters referring to a pure gas are nowdescribed 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 andrefer to the density and the molar fraction of the adsorbed phase at distance z from theadsorbent surface. The function fsolve() now have to solve a system of two equations. To15roceed, we define a function f(u) , with u an array containing the gas density d_z and the molarfraction of the first component x_z_A to be found. This function returns an array containingthe two equations to be solved. The solver now requires initial guesses in an array format tomatch the system of equation, which is given by [d0,x0[0]] . Remark that for a mixture oftwo gases, we only used the molar fraction of the first component and get the molar fractionof the second component by . In the function f(u) , the Decimal() function is used toredefine the molar fraction to x_z_A and x_z_B such that ∼ − . Finally, thefunction d_mix() return an array composed of the gas density and the molar fraction of thefirst component.This d_mix() function induced modifications to the N_ex_pure function that becomes def
N_ex_mix(T, z0, eps0, dB, xB):N = 300delta = z0 / Nd0 = dBx0 = xBintegral_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*deltaintegral_B += d_z*(1 - x_z_A)*deltad0 = d_zx0 = [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 cansee, the terms d_z*x_z_A*delta and d_z*(1-x_z_A)*delta give the absolute adsorption forcomponent A and B respectively, in the interval delta . The initial guesses d0 and x0 , for densityand 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 adsorptionminus the term dB*xB[0]*z0 or dB*xB[1]*z0 .The mix_fit() function for mixtures take the form def mix_fit(params):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 = lmfit.Parameters()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)) / CO mixture on Calgon F-400 at 318.2K [23], our fitted parameters are givenin Table 4. Figure 5 represents the fit for pure gases with a common z0 . The mean errorbetween experimental data pure isotherms and the fitted MPTA model are 4 . / . / CO mixtureon Calgon fitted from pure gas isotherms.Parameter Initial value Fitted value Fit error % z ,common ( cm /g ) 1 0.300 1 . ε ,CH ( J/mol ) 10 000 7475 3 . ε ,CO ( J/mol ) 10 000 7767 3 . / 20% CO mixture. The figure shows the excess adsorption of each component andthe total excess adsorption for the mixture. The error bars represent the expected experimentaluncertainties on the measurements given in Ref. [23]. The mean error between experimentaldata and the model are presented in table 5. It is worth noting that the mean error of theMPTA model is of the same order as the experimental uncertainty on datasets. E x e ss a d s o r p t i o n ( m o l / k g ) Error: CH4 4.62 % CO2 4.51 % ǫ ,CH = 7475 J/mol, ǫ ,CO = 7767 J/mol, z ,common = 0 . cm /g. CH4CO2
Figure 5: Result of the fit for pure CH and CO on Calgon F-400 at 318.2K.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 feedingmolar composition, which differs from the measured molar fraction by at most 10%. Theresulting mean error between experimental and calculated values increases from 4 . . .
67% to 25%, 71% and 4 .
15% for the CH component, the CO component and themixture, respectively. 17 .0 0.2 0.4 0.6 0.8 1.0Ratio z/z00.00.20.40.60.81.0 A d s o r b e d p h a s e m o l a r f r a c t i o n ρ B = 0 . mol/L, ǫ ,CH = 7475 J/mol, ǫ ,CO = 7767 J/mol, z ,common = 0 . cm /g. CH4CO2
Figure 6: Adsorbed phase molar fraction of each component for a80% CH / 20% CO mixture on Calgon F-400 at 318.2K.Table 5: Mean error between experimental data and the MPTA model fora 80% CH / 20% CO mixture on Calgon.Model mean error (ADDn) Experimental uncertainty [23]CH E x e ss a d s o r p t i o n ( m o l / k g ) Error: CH4 4.42 % CO2 4.41 % Mixture 3.67 %
CH4 componentCO2 component80% CH4 / 20% CO2 Mixture
Figure 7: Excess adsorption for a 80% CH / 20% CO mixture on Calgon F-400 at 318.2K.18 Conclusion
We presented a detailed numerical implementation of the MPTA adsorption model, coveringpure gases and binaries mixtures. From our treatment of the binary mixture case,generalization to more complex mixtures is straightforward. This implementation includessome of the strategies used to obtain a robust model even in case of fast variation of fluiddensities, like phases transitions. For all the mixtures in Ref. [23], the mean error betweenpredictions of the MPTA model for mixture (Gibbs) excess adsorption and experimental datavary from 1 .
44% (20% CH / 80% CO mixture) to 5 .
09% (20% N / 80% CO mixture).Meanwhile, the experimental uncertainties for mixture adsorption vary from 2 .
4% (80%CH / 20% N and 40% CH / 60% N mixtures) to 4 .
1% (80% N / 20% CO 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 theimplementation and also the Natural Sciences and Engineering Research Council of Canadaand the Savannah River National Laboratory for financial support.
References [1] O. Talu, “Needs, status, techniques and problems with binary gas adsorptionexperiments,”
Advances in Colloid and Interface Science , vol. 76, pp. 227–269, 1998.[2] W. Henry, “Experiments on the quantity of gases absorbed by water, at differenttemperatures, and under different pressures,”
Philosophical Transactions of the RoyalSociety of London , vol. 93, pp. 29–276, 1803.[3] H. Freundlich, “Over the adsorption in solution,”
Journal of Physical Chemistry , vol. 57,no. 385471, pp. 1100–1107, 1906.[4] I. Langmuir, “The adsorption of gases on plane surfaces of glass, mica and platinum,”
Journal of the American Chemical society , vol. 40, no. 9, pp. 1361–1403, 1918.[5] S. Brunauer, P. H. Emmett, and E. Teller, “Adsorption of gases in multimolecular layers,”
Journal of the American chemical society , vol. 60, no. 2, pp. 309–319, 1938.[6] M. Temkin and V. Pyzhev, “Recent modifications to langmuir isotherms,”
Acta Physico-Chimica Sinica , 1940.[7] R. Sips, “On the structure of a catalyst surface,”
The Journal of Chemical Physics , vol. 16,no. 5, pp. 490–495, 1948.[8] R. Sips, “On the structure of a catalyst surface. II,”
The Journal of Chemical Physics ,vol. 18, no. 8, pp. 1024–1026, 1950.[9] P. Kisliuk, “The sticking probabilities of gases chemisorbed on the surfaces of solids,”
Journal of Physics and Chemistry of Solids , vol. 3, no. 1-2, pp. 95–101, 1957.[10] J. Toth, “State equations of the solid-gas interface layers,”
Acta chimica AcademiaeScientiarum Hungaricae , vol. 69, no. 3, pp. 311–328, 1971.1911] M. Dubinin and V. Astakhov, “Development of the concepts of volume filling of microporesin the adsorption of gases and vapors by microporous adsorbents,”
Russian ChemicalBulletin , vol. 20, no. 1, pp. 3–7, 1971.[12] D. M. Ruthven, “Sorption of oxygen, nitrogen, carbon monoxide, methane, and binarymixtures of these gases in 5a molecular sieve,”
AIChE Journal , vol. 22, no. 4, pp. 753–759,1976.[13] S. Suwanayuen and R. P. Danner, “A gas adsorption isotherm equation based on vacancysolution theory,”
AIChE Journal , vol. 26, no. 1, pp. 68–76, 1980.[14] S. Doong and R. Yang, “A simple potential-theory model for predicting mixed-gasadsorption,”
Industrial & engineering chemistry research , vol. 27, no. 4, pp. 630–635,1988.[15] A. Myers, “Thermodynamics of adsorption in porous materials,”
AIChE Journal , vol. 48,no. 1, pp. 145–160, 2002.[16] A. A. Shapiro and E. H. Stenby, “Potential Theory of Multicomponent Adsorption,”
Journal of Colloid and Interface Science , vol. 201, no. 2, pp. 146–157, 1998.[17] M. Polanyi, “The Potential Theory of Adsorption,”
Science , vol. 141, no. 3585, pp. 1010–1013, 1963.[18] E. Dundar, R. Zacharia, R. Chahine, and P. B´enard, “Modified potential theory formodeling supercritical gas adsorption,”
International Journal of Hydrogen Energy , vol. 37,no. 11, pp. 9137–9147, 2012.[19] M. M. Dubinin, “Fundamentals of the theory of adsorption in micropores of carbonadsorbents: characteristics of their adsorption properties and microporous structures,”
Pure and Applied Chemistry , vol. 61, no. 11, pp. 1841–1843, 1989.[20] N. D. Hutson and R. T. Yang, “Theoretical basis for the Dubinin-Radushkevitch (DR)adsorption isotherm equation,”
Adsorption , vol. 3, no. 3, pp. 189–195, 1997.[21] M. Dubinin, “Generalization of the theory of volume filling of micropores tononhomogeneous microporous structures,”
Carbon , vol. 23, no. 4, pp. 373–380, 1985.[22] M. A. Monsalvo and A. A. Shapiro, “Study of high-pressure adsorption from supercriticalfluids by the potential theory,”
Fluid Phase Equilibria , vol. 283, no. 1, pp. 56–64, 2009.[23] M. Sudibandriyo, Z. Pan, J. E. Fitzgerald, R. L. Robinson, and K. A. Gasem, “Adsorptionof methane, nitrogen, carbon dioxide, and their binary mixtures on dry activated carbonat 318.2 k and pressures up to 13.6 mpa,”
Langmuir , vol. 19, no. 13, pp. 5323–5331, 2003.[24] O. Redlich and J. N. S. Kwong, “On the thermodynamics of solutions; an equation ofstate; fugacities of gaseous solutions.,”
Chemical reviews , vol. 44, no. 1, pp. 233–244, 1949.[25] G. Soave, “Equilibrium constants from a modified Redlich-Kwong equation of state,”
Chemical Engineering Science , vol. 27, no. 6, pp. 1197–1203, 1972.[26] E. Lemmon, M. Huber, and M. McLinden, “NIST Standard Reference Database 23:Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 9.1,National Institute of Standards and Technology, Standard Reference Data Program,Gaithersburg,” 2013. 2027] U. Setzmann and W. Wagner, “A new equation of state and tables of thermodynamicproperties for methane covering the range from the melting line to 625 K at pressures upto 100 MPa,”
Journal of Physical and Chemical reference data , vol. 20, no. 6, pp. 1061–1155, 1991.[28] R. Span, E. W. Lemmon, R. T. Jacobsen, W. Wagner, and A. Yokozeki, “A referenceequation of state for the thermodynamic properties of nitrogen for temperatures from63.151 to 1000 K and pressures to 2200 MPa,”
Journal of Physical and Chemical ReferenceData , vol. 29, no. 6, pp. 1361–1433, 2000.[29] R. Span and W. Wagner, “A new equation of state for carbon dioxide covering the fluidregion from the triple-point temperature to 1100 k at pressures up to 800 mpa,”
Journalof physical and chemical reference data , vol. 25, no. 6, pp. 1509–1596, 1996.[30] B. J. Thelen, “Python (3.x & 2.x) API for NIST Standard Reference Database 23 (a.k.a.Refprop) on Linux and Windows.” Accessed: 2016-11-11.[31] B. Wernick, “Linking REFPROP with Python Applications.” Accessed: 2016-11-11.[32] T. E. Oliphant, “Python for scientific computing,”
Computing in Science & Engineering ,vol. 9, no. 3, pp. 10–20, 2007.[33] K. J. Millman and M. Aivazis, “Python for scientists and engineers,”
Computing inScience & Engineering , vol. 13, no. 2, pp. 9–12, 2011.[34] E. Jones, T. Oliphant, P. Peterson, et al. , “Open source scientific tools for python,” 2001.[35] S. v. d. Walt, S. C. Colbert, and G. Varoquaux, “The numpy array: A structure forefficient numerical computation,”
Computing in Science & Engineering , vol. 13, no. 2,pp. 22–30, 2011.[36] F. P´erez and B. E. Granger, “Ipython: A system for interactive scientific computing,”
Computing in Science & Engineering , vol. 9, no. 3, pp. 21–29, 2007.[37] J. D. Hunter, “Matplotlib: A 2d graphics environment,”
Computing in Science &Engineering , vol. 9, no. 3, pp. 90–95, 2007.[38] M. Newville, T. Stensitzki, D. B. Allen, and A. Ingargiola, “LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python,” Sept. 2014.[39] IEEE Computer Society, “IEEE 754-2008 Standard for floating-point arithmetic,” 2008.[40] F. Dreisbach, R. Staudt, and J. Keller, “High pressure adsorption data of methane,nitrogen, carbon dioxide and their binary and ternary mixtures on activated carbon,”
Adsorption , vol. 5, no. 3, pp. 215–227, 1999.[41] M. G. Bjørner, A. A. Shapiro, and G. M. Kontogeorgis, “Potential theory of adsorption forassociating mixtures: possibilities and limitations,”