J. Phys. Chem. C 2007, 111, 7463-7472
7463
On Inverse Adsorption Chromatography. 1. Solution for Nonlinear, Nonideal Model of Chromatography and Adsorption Isotherm V. A. Bakaev The PennsylVania State UniVersity, Materials Research Institute, UniVersity Park, PennsylVania ReceiVed: NoVember 21, 2006; In Final Form: February 21, 2007
A mathematical model of nonlinear, nonideal adsorption chromatography of one component solute (adsorbate) which takes into account kinetics of adsorption, longitudinal diffusion in the mobile phase, and compressibility of the mobile phase (carrier gas) is formulated and solved by the method of perturbation. The solution is obtained in an analytical form and consists of two terms. The first term of the solution corresponds to nonlinear, ideal chromatography which is the foundation of the method called elution by characteristic points (ECP). At present, ECP is the basic method of the inverse adsorption chromatography which, however, neglects both the kinetics of adsorption and longitudinal diffusion. The next term of the solution gives the first order corrections to ECP. The application of the method is only briefly discussed in this paper and is explained in detail in the accompanying experimental paper (Bakaev, V. A.; Bakaeva, T. I.; Pantano, C. G. J. Phys. Chem. C 2007, 111, 7473-7486.)
1. Introduction This paper is the first part of a two part communication. It deals with the mathematical model of nonlinear, nonideal chromatography. The second part is mainly experimental. In what follows it will be referred to as II (see ref 25). The second part considers the experimental application of the solution obtained in the first part, as well as other problems of inverse adsorption chromatography and surface heterogeneity. Inverse adsorption chromatography is an experimental technique which allows one to obtain adsorption characteristics such as isotherms and heats of adsorption from chromatographic experiments. The history of determination of adsorption isotherms and heats of adsorption from chromatographic experiments is briefly reviewed by Kiselev and Yashin.1 In fact, they considered only one method which they called Glueckauf’s method.2 Now the method1 is known as the elution by characteristic points (ECP) method. A more complete analysis of chromatographic methods for measuring adsorption isotherms based on the concept of characteristic point allows one to compare the ECP method with other methods of inverse adsorption chromatography.3 Again Glueckauf was named as the first person to use chromatography for measuring isotherms of adsorption, but the priority of using specifically the ECP method in gas chromatography was ascribed to Cremer and Huber.4 Finally, the recent progress in measuring adsorption isotherms (mainly from liquids) by chromatographic means has been recently reviewed by Guiochon, Felinger, Shirazi, and Katti (see Chapter 3, Section 3.5 in ref 5). The term ECP was coined by Conder and Purnell.6 In fact, their main concern in that paper was another method which they called frontal analysis by characteristic point (FACP). Mathematical models behind these two methods are very similar,5 that for FACP being slightly simpler than for ECP. However, as mentioned by Conder and Purnell, ECP has the merit of simplicity of operation.6 This means that in the case of ECP one can usually employ a conventional chromatograph without any modification, whereas in the case of FACP, the
injection system of at least the conventional gas chromatograph (GC) should be modified. In this paper, we concider only ECP. The computational part of the original ECP method is very simple (cf. discussion of eq 31), but this simplicity is achieved at the expense of sweeping approximations. Besides, not all types of adsorption isotherms can be directly measured by the ECP method but only those without inflection points.1,3 The ECP method is based on the nonlinear ideal chromatography model. The model neglects the role of adsorption kinetics (mass transfer in the mobile phase) as well as molecular and eddy diffusion. These and other assumptions of the ECP method are discussed by Roles and Guiochon.7 We applied the ECP method for studying the energy and entropy distributions of adsorption sites for butanol on the glass fiber surfaces (see the Introduction in II). Since the specific surface areas of glass fibers are very small (about 0.2 m2/g), only the strong adsorption sites with the long residence (adsorption) times of adsorbed molecules could be probed by the method. The energy distribution of these sites is obtained from the lower concentration region of an adsorption isotherm and the ECP method is ideally suited for probing that region (see p 435 in ref 3). However, the energy distributions obtained by the original ECP method displayed unrealistic (double valued) shape in the domain of the strongest adsorption sites (cf., for example, Figure 1 in ref 8). It is known that the ECP method is accurate only if the adsorptive is not too strongly adsorbed. Otherwise adsorption and desorption may become too slow for equilibrium to be attained and the system departs from ideal chromatographic behavior (see p 435 in ref 3). The concentration (partial pressure) of the solute in the mobile phase which is in equilibrium with those strong sites is very low which also can slow down the equilibration in the column. For that reason, we initially assumed that it is the too slow kinetics of adsorption on the strongest adsorption sites that is responsible for the above-mentioned artifacts of the energy distributions.8 Thus, initially the primary goal was to modify the ECP method to take account of possibly slow adsorption kinetics at
10.1021/jp067733b CCC: $37.00 © 2007 American Chemical Society Published on Web 05/02/2007
7464 J. Phys. Chem. C, Vol. 111, No. 20, 2007 low solute concentrations in the mobile phase. This problem is partly resolved in this paper where a modified ECP method is described which takes account not only for adsorption kinetics but also for other nonideality factors neglected by the original ECP. It is shown in II that the adsorption characteristics obtained by the modified ECP method deviate considerably from those obtained by the original ECP method. However, the artifact mentioned above was not eliminated. It turned out that it has another source. The ECP method (original or modified) does not allow one to obtain an isotherm of adsorption but only its first derivative. To obtain an isotherm, one has to integrate its derivative, which gives rise to the problem of the integration constant. It is the integration constant that is the source of the artifact mentioned above (see section 2.1 in II). The mathematical models of nonlinear nonideal chromatography are comprehensively reviewed in conjunction with the problems of preparative chromatography.9 The peculiarity of analytic and preparative adsorption chromatography is that they deal with adsorption mainly on homogeneous surfaces. A typical isotherm of adsorption on homogeneous surfaces is the Langmuir isotherm. That isotherm is approximately linear at small concentration which is also the case for other isotherms on homogeneous surfaces. Accordingly, the chromatography of infinitesimal concentrations (analytical chromatography) is described by linear models, whereas that of finite concentrations (preparative chromatography) is described by nonlinear models.9 In this and the company paper II, we consider adsorption on heterogeneous surfaces. A typical isotherm of adsorption on heterogeneous surfaces is the Freundlich isotherm. That isotherm is never employed in preparative or linear chromatography because its slope tends to infinity when concentration (pressure) tends to naught. This means that “...band elution never ends and the peak tails forever” (see p 17 in ref 9). The fact that the slope of the Freundlich isotherm tends to infinity at small pressures (concentrations) is not senseless. It indicates that the domain of pressures where the slope of an isotherm becomes approximately constant (isotherm becomes approximately linear) cannot be reached under certain conditions. That was shown for adsorption of rare gases, nitrogen, and other low boiling gases on pyrex and some other heterogeneous surfaces.10 Isotherms of these gases were not linear but obeyed either the Freundlich or Dubinin-Radushkevich equations even in the ultrahigh vacuum.10 With respect to the chromatography peak tailing, it means that although it will not last strictly speaking forever one should not expect it to terminate for any reasonably long time. Thus, the chromatography model employed below has one important distinction from those used in either preparative chromatography5,9 or even in the original ECP:3 Our model is strongly nonlinear even at very small concentrations. The various kinetic models of nonlinear chromatography consider diffusion across the mobile phase stream, resistance to mass transfer across the external film around the particle, intraparticle diffusion, and kinetics of adsorption-desorption.5,9,11 The intraparticle diffusion should be excluded from this list because we consider adsorption on nonporous particles (glass fibers8 or nonporous nanoparticles such as fumed silica II). Since the role of other kinetics mechanisms mentioned above is not clear for our case of adsorption at low solute concentrations on nonporous particles in a chromatographic column, we prefer a lumped model of adsorption kinetics. There are two such models in chromatography. These are the solid film linear driving force model and the liquid film linear driving force model (see p 39 in ref 9). The former model was introduced by Glueckauf and Cotes12 and used later by many authors,9,11
Bakaev whereas the latter model was introduced somewhat earlier by Zabezhinski et al. (see ref 13 and references therein) and was not used as frequent as the former one. It was experimentally substantiated13 and employed for modeling of gas adsorbers.14 The mathematical part of the work14 which is essentially the frontal analysis (FA) model5 was later included in the monograph on the PDEs of mathematical physics (see eq 2-6.86 in ref 15). Besides, the liquid film linear driving force model is used as a basic model for adsorption kinetics in the monograph on the first-order PDEs (see eq 1.2.2 in ref 16). We select the liquid film linear driving force model as an empirical kinetics equation because the driving force of the model is the deviation of the solute concentration in the mobile phase from the equilibrium solute concentration corresponding to adsorption in the immobile phase (cf. eq 1b). The stronger an adsorption site is, the lower its equilibrium solute concentration is (at a given adsorption and temperature). Thus, when adsorption in the column occurs on the strongest adsorption sites, the equilibrium solute concentration may be very low and the concentration of solute in the mobile phase which is usually proportional to the equilibrium solute concentration is also low. This leads to the slowdown of adsorption kinetics on the stronger adsorption sites due to the decrease of the driving force for the mass transfer. In contrast to that, the driving force of the solid film linear driving force model is the deviation of the equilibrium adsorption corresponding to the solute concentration in the mobile phase from the actual adsorption in the immobile phase. Since the low solute concentration in the mobile phase does not necessarily mean that adsorption on strong adsorption sites is also low, it does not necessarily mean that the driving force of the mass transfer in the solid film linear driving force model is small. Therefore, the reason for the slowdown of the adsorption kinetics on the strong adsorption sites seems to be more explicit in the liquid film linear driving force model than in the solid film linear driving force model. Since, as mentioned above, the main motivation of this study was an attempt to take that slowdown of adsorption kinetics into account, the liquid film linear driving force model seems to be more adequate to the problem. An attractive feature of the original ECP method is that it provides a simple analytical solution both for the direct and the inverse problem of the nonlinear ideal chromatography model. The former problem is the determination of an elution profile (EP) from a known adsorption isotherm, whereas the latter problem is the determination of an unknown adsorption isotherm from EP. Inclusion of nonideality factors in the ideal chromatography model complicates considerably even the above-mentioned direct problem. Numerous numerical solutions were published.5,14,17 In particular, the effect of the mass transfer coefficient in the linear driving force model on the EP of nonlinear chromatography was discussed by numerical and analytical means.9,18 It was shown that when the kinetics of adsorption is sufficiently fast its influence on the EP of nonlinear chromatography is the same as that of longitudinal diffusion. This is a well-known fact for the linear chromatography.9,11 The influence of the longitudinal diffusion on the EP in nonlinear chromatography was extensively studied by numerical methods.5,9,17 Even more computationally involved is the inverse problem of nonideal nonlinear chromatography whose solution by the method of iterations is only discussed on the grounds of some numerical illustrations.19 A real (and yet unresolved) problem here is the convergence of iterations. Another approach is reviewed in Chapter 3, Section 3.5.6 of ref 5.
On Inverse Adsorption Chromatography. 1.
J. Phys. Chem. C, Vol. 111, No. 20, 2007 7465
In what follows, we consider a model of nonlinear adsorption chromatography which takes into account the adsorption kinetics, the longitudinal diffusion, and the compressibility of the carrier gas (for the gas chromatography). The method of its solution (section 3) is based on expansion of the concentration of the solute (adsorbate) in the mobile and stationary phases into powers of the factors which determine the kinetics of adsorption, diffusivity, and compressibility of the carrier gas. The first term of this expansion is ECP. The second term which is obtained here in an analytical form gives the first order corrections to ECP. The third and the following terms of the expansion are truncated. Then in Section 4 we discuss how those analytical formulas can be used for obtaining adsorption isotherms from EP (inverse problem of chromatography mentioned above) by iterations similar to those considered earlier.19 Finally, in section 5, we illustrate the analytical formulas obtained in section 3 and the convergence of iterations described in section 4 by some schematic numerical examples. A detailed description of the method’s application can be found in the accompanying experimental paper II. 2. Description of the Basic Model The basic PDE for nonlinear, nonequilibrium, dispersive model of chromatography is well-known5,9 / ∂2C* ∂C* Vs ∂q* ∂u*C* + / + ) D* ∂t* V ∂t* ∂z* ∂z*2
(1a)
m
It can be supplemented by the lumped kinetics equation15 (see also eq 1.2.2 of ref 16) which is now called the liquid film linear driving force model9
τ/kin
∂q* ) C* - Ceq(q*) ∂t*
(1b)
with the initial conditions
C*(z*,0) ) 0; q*(z*,0) ) 0
(2)
and the boundary condition
C*(0,t*) ) φ(t*)
(3)
Here and below the starred variables designate dimensional variables, whereas the same variables without a star as a superscript are their dimensionless counterparts determined below (cf. eq 7). In eq 1, t* > 0 is time; z* > 0 is the coordinate along the column; C* is the concentration of a one component solute in the mobile phase; q* is the concentration of the solute in the stationary phase; u* is the linear velocity of the mobile phase; V/m and V/s are volumes of the column occupied by the mobile and stationary phase, respectively; D* is the longitudinal diffusion coefficient; and τ/kin is a parameter which determines the kinetics of adsorption. Finally, Ceq(q*) is a function inverse to an equilibrium adsorption isotherm (in the units of eq 1a). The boundary conditions at the inlet of the column are determined (for the gas chromatography) by the following ordinary differential equation (ODE) which describes schematically the elution of solute from an injector by the flow of the carrier gas
V/in
dC*(0,t*) ) -F*C*(0,t*) dt* C*(0,0) )
C/max
)
N/in V/in
(4a) (4b)
and has the solution
(
F* C*(0,t*) ) φ(t*) ) C/max exp - / t* Vin
)
(5)
Here V/in is the volume of an injector at the beginning of the column and F* is the flow rate of the carrier gas (volume per unit of time). Finally, N/in is the injected quantity of solute which together with V/in determines C/max, maximal concentration of the solute in the mobile phase (see below). Now we briefly discuss the physical meaning of these equations. When D* ) 0, eq 1a expresses the conservation of the solute in the column. It differs from the corresponding mass balance eq 9.9 of ref 3 only in that u* in eq 1a is supposed to be dependent on z* due to the carrier gas compressibility. Another distinction from the model3 is that we neglect the socalled sorption effect in the gas chromatography. The effect arises when the concentration of solute in the carrier gas is comparable to the concentration of the carrier gas itself. We neglect the sorption effect in our model because C* is supposed to be small, much smaller than the concentration of the carrier gas. The latter statement calls for explanation. Usually the chromatography of small solute concentrations is the linear chromatography, whereas the finite solute concentrations chromatography which takes into account the sorption effect is the nonlinear chromatography.3 As explained in the Introduction, this is true because the analytical and preparative adsorption chromatography is carried out on homogeneous or mildly heterogeneous surfaces. In this paper, we consider a chromatography (in fact, inverse chromatography) on strongly heterogeneous surfaces where the adsorption isotherm and correspondingly the chromatography model are nonlinear even at very small concentrations (see the Introduction). This is the reason why our model is nonlinear at sufficiently small solute concentrations when the sorption effect is negligibly small. A mass balance eq 1a with D* ) 0 is still not rigorous. It neglects the molecular diffusion of solute molecules in the mobile phase and the so-called eddy diffusion.5,9 It also assumes that the areas of cross sections of the mobile phase and the stationary phase are the same all over the column. In fact in a packed column, their ratio fluctuates around its mean value V/s / V/m in eq 1a. The lump parameter D* on the right-hand side of eq 1a which is called a coefficient of longitudinal diffusion (apparent diffusion coefficient, dispersion factor) takes account of all these stochastic factors. The molecular diffusion coefficient depends on the properties of a solute molecule (molecular weight, etc.; see p 96 in ref 20), but other components do not depend on them. This allows one to estimate the value of D* in eq 1a (and that of V/in in eq 4) from the elution peak of a non adsorbing gas (see eq 6 in II). If the value of D* is known, the value of τ/kin in eq 1b can be selected so that EP of the solute do not depend on the flow rate of the carrier gas (see discussion of Table 1 in II). As explained in the Introduction, eq 1b is an empirical equation which is, at best, only semiquantitatively correct. Finally it should be mentioned that eq 4 is an idealization of the processes that occur in a real injector. The boundary condition in eq 3 is also an idealization, a more accurate boundary condition being, for example, eq 4 in ref 11. Designate now L* as the length of the column and t/m as a time which is required for a non adsorbing solute to trace the column (hold-up time, dead time) and take these variables as the units of length and time, respectively. Take also C/max eq 4b
7466 J. Phys. Chem. C, Vol. 111, No. 20, 2007
Bakaev
as the unit of the solute concentration in the mobile phase and V/mC/max as the unit of the quantity of adsorbate. Since
V/m ) jF/ot/m
V/s q* D*t/m z* C* t* ; z ) ; N ) ; D ) ; ; C ) L* t/m C/max V/mC/max (L*)2 τ/kinV/m u* / u) t ; τ ) / / (7) L* m kin t V m s
In these variables, eqs 1-5 take the form
∂2C ∂C ∂N ∂uC + + )D 2 ∂t ∂t ∂z ∂z τkin
(8)
∂N ) C - Ceq(N) ∂t
C(z,0) ) 0; N(z,0) ) 0; z > 0
( )
t C(0,t) ) φ(t) ) exp τ0
(9)
τkin x1 + u2 τkin
dC ) -ψ(C,N) ds
(14a)
dN ) ψ(C,N) dt
(14b)
where ψ(C,N) ) C - Ceq(N). This is a system of equations in full derivatives and in this respect it resembles a system of two ODEs. However, in contrast with the latter the derivatives in eq 14 are taken in different directions. The derivative with respect to the s variable eq 14a is taken along the straight lines defined by the equation z ) ut - z0, whereas the derivative with respect to t in the eq 14b is taken along the line z ) z0. These two directions are the characteristics of eq 11. When τkin is sufficiently large, this hyperbolic system of two PDEs can be easily solved numerically by the method of characteristics.21 However, when τkin in eq 11 becomes very small, the system becomes rather difficult to solve numerically. The situation resembles that known for the stiff systems of ODEs which are notoriously difficult for numerical integration. One method used for their solution is the perturbation theory22 which we will apply to our case. This method has already been used for solution of a bilinear chromatography model (see section 5 in ref 23). However, the difference between the model considered earlier23 and that considered in this paper is so considerable that it is difficult even to compare the final results. Since the parameter τkin is assumed to be small, expand the dependent variables of eq 11 in its powers
C ) C0 + τkinC1 + τkin2C2 + ...
(10)
N ) N0 + τkinN1 + τkin2N2 + ...
where t > 0 and
τ0 )
(13)
Equations 11 can be presented as two ODEs
(6)
where j is the James and Martin pressure correction factor (see pp 41 and 42 in ref 20) and F/o is the carrier gas flow rate at the outlet of a column, all of the units determined above are experimentally measurable variables. Equation 6 is again an approximation. It assumes the so-called dead volumes (cf. p 48 in ref 20) are much smaller than the gas volume in the column. This condition is observed in our loosely packed columns filled either with glass fibers or with nanoparticles diluted by much larger inert particles where j is about 0.98 (see II). We now introduce dimensionless variables
t)
C(0,t) ) φ(t)
(15)
Substitute eq 15 into eq 11a and nullify coefficients at each power of τkin to obtain
V/in jF/ot/m
(10a) jF/o
The value of F* in eq 5 is substituted for in eq 10a. This is because eq 5 and the expression for moments of a chromatographic peak which follow from it were obtained for an incompressible mobile phase (see II). For a compressible carrier gas and small values of j, a reasonable approximation for F* can be the average value of the flow rate which is jF/o.
∂Ck ∂Ck ∂Nk + +u ) 0; ∂t ∂t ∂z
k ) 0, 1, 2, ...
(16)
∑ τkkinCk - Ceq(kg0 ∑ τkkinNk)
(17)
Substitute eqs 15 in eq 11b to obtain
∑ τk+1 kin kg0
∂Nk ∂t
)
kg0
Expand the second term on the right-hand side of eq 17 in powers of τkin
3. Solution of the Direct Problem 3.1. Kinetics of Adsorption. Assume first that D ) 0 and u does not depend on z (in this approximation, u ) 1 but we formally keep it as a factor to preserve the familiar form of equations). In this case, eqs 8-10 reduce to the following system of PDEs:
Ceq(
τkkinNk) ) Ceq(N0) + τkinC′eq(N0)N1 + ... ∑ kg0
(18)
Insert eq 18 into eq 17 and nullify coefficients at powers of τkin. The coefficient at 0th and first power of τkin give
∂C ∂N ∂C + +u )0 ∂t ∂t ∂z
(11a)
0 ) C0 - Ceq(N0)
(19a)
∂N ) C - Ceq(N) τkin ∂t
(11b)
∂N0 ) C1 - C′eq(N0)N1 ∂t
(19b)
Since Ceq(N) is the inverse of an adsorption isotherm Neq(C), it follows from eq 19a that
with the initial condition (z > 0)
C(z,0) ) 0; N(z,0) ) 0 and the boundary condition (t > 0)
(12)
N0 ) Neq(C0) and C′eq(N0) )
( )
dN0 1 ) dC0 N′eq(C0)
-1
(20)
On Inverse Adsorption Chromatography. 1.
J. Phys. Chem. C, Vol. 111, No. 20, 2007 7467
Substitute the first of eq 20 into eq 16 with k ) 0 to obtain eq 21a. Equation 21b can be obtained by substituting N1 from eq 19b into eq 16 with k ) 1 and using the second equation of eq 20
f(C0) f(C0)
∂C0 ∂C0 +u )0 ∂t ∂z
(21a)
( ) (
)
∂C1 ∂C1 ∂ dN0 ∂ dN0 ∂N0 +u ) -C1 + (21b) dt ∂z ∂t dC0 ∂t dC0 ∂t
In these equations
f(C0) ) 1 +
dN0 dNeq )1+ dC0 dC0
(22)
Since C0 depends on t (and z), the derivatives on the righthand side of eq 21b are
(
)
( )
∂C0 ∂ dN0 ) f ′(C0) ∂t dC0 ∂t
(23a)
( )
∂C0 2 ∂ dN0 ∂N0 ) 2f ′(C0)[f(C0) - 1] + ∂t dC0 ∂t ∂t ∂2C0 (23b) [f(C0) - 1] ∂t2 2
Equations 21a and 21b may be supplemented by the following boundary and initial conditions that agree with eqs 9 and 10
C0(0,t) ) φ(t); C1(0,t) ) 0
(24)
C0(z,0) ) 0; C1(z,0) ) 0
(25)
where z > 0 and t > 0. Equations 21a and 21b are the PDEs of the first order. They can be solved sequentially as follows. First one can solve the quasilinear homogeneous eq 21a with the first equation in eq 24 and to find C0(z,t). With C0(z,t) being a known function, eq 21b (with eq 23) is a strictly linear nonhomogeneous PDE which together with the second equation in eq 24 can be solved for C1(z,t) by the method of characteristics. If C0(z,t) is a solution of eq 21a, then
(∂z∂t)
)-
C0
∂C0 ∂C0 f(C0) / ) ∂z ∂t u
(26)
One can integrate eq 26 at constant C0 to obtain
t ) t0 +
f(C0) z u
(27)
Here t0 is just a constant of integration. Its physical meaning becomes clear from the following: As seen from eq 27, t ) t0 when z ) 0. Thus, t0 designates the time which determines the boundary condition at z ) 0. Since the latter is given by eq 10,
( )
t0 C0 ) φ(t0) ) exp τ0
(28)
z u
dNeq )t-1 dC0
(29)
(30)
Now use eqs 7 to restore the dimensionality of variables in eq 30 and also use eq 6 to obtain
dN/eq ) jF/o(t* - t/m) dC*
(31)
This is the ECP equation for obtaining an adsorption isotherm N/eq(C*)from inverse gas chromatography (IGC).1,3 Here t* depends on C* so that t*(C*) is EP at the outlet of a column. Equation 31 explicitly solves the inverse problem of adsorption chromatography, calculation of adsorption isotherm from EP. It does not require calculation of the shock trajectory. (The complete solution of the direct problem, calculation of EP from the adsorption isotherm, would require the calculation of the shock trajectory.) This convenient property of the ECP method of IGC is retained in the modified ECP method considered below. The modified ECP method considered below requires calculation of only the first order corrections to C0(z,t), in particular C1(z,t) in eq 15. The function C1(z,t) which satisfies eqs 21b and 24 can be found by the method of characteristics as follows. If z and t in C1(z,t) depend on a parameter s, then
dC1 ∂C1 dt ∂C1 dz ) + ds ∂t ds ∂z ds
and eq 27 takes the form
t ) -τ0 ln(C0) + f(C0)
On one hand, this equation implicitly determines C0(z,t) as the solution of eq 21a with the first equation in eq 24 as the boundary condition. On the other hand, it determines the characteristic of eq 21a labeled by C0 (or by t0 eq 28). (The method of characteristics is used below for solution of eq 21b.) The first equation in eq 25 (the initial condition for C0(z,t)) is not observed by eq 29, and because of this, C0(z,t) implicitly determined by eq 29 represents only a part of the integral surface determined by eqs 21a, 24, and 25. The total integral surface contains a discontinuity (shock).9,16 The trajectory of the shock t ) tsh(z) starts at the origin of the first quadrant of the (z,t) plane (for the boundary condition of eq 10 and f ′(C0) < 0) and divides this quadrant into two parts. The solution given by eq 29 determines the integral surface only above the shock trajectory when t > tsh(z). Below the shock trajectory, the value of C0 is identically zero which satisfies both eq 21a and the initial conditions in eq 25. The boundary condition eq 24 need not be satisfied in this region because below the shock trajectory z > 0 for all values of t except only one point: the origin of coordinate system. Thus, the integral surface C0(z,t) that is a solution of eqs 21a, 24, and 25 consists of two parts. One of them (eq 29) satisfies eq 21a and the boundary condition of eq 24, whereas another (C0(z,t) ≡ 0) satisfies eq 21a and the initial condition of eq 25. These parts are separated by the shock trajectory t ) tsh(z) along which C0(z,t) has a jump. This solution is unique and called a weak solution.9,16 The weak solution is the foundation of the ECP method mentioned in the Introduction. This is the zero approximation of our model. Since usually t0 , t, one can neglect t0 in eq 27 and obtain the ECP approximation for C0(z,t).1,3 In this approximation, at z ) 1 and u ) 1 f(C0) ) t and according to eq 22
(32)
This equation coincides with eq 21b if the following system of ODEs is observed
7468 J. Phys. Chem. C, Vol. 111, No. 20, 2007
( )
dC1 dz dt ∂ dN0 + ) f(C0); ) u; ) -C1 ds ds ds ∂t dC0
Bakaev where
(
)
B(z;C0) ) 1 - z
∂ dN0 ∂N0 (33) ∂t dC0 ∂t where C0(z,t) is implicitly determined by eq 29 and N0(C0) is determined by eq 20. The plot of t vs z parametrically determined by the first two equations of this ODEs system is called a characteristic. One can exclude parameter s from the first two equations of this system to obtain
dt f(C0) ) dz u
(34)
If C0 is kept constant, eq 34 can be integrated to give eqs 27 and 29. This means that the systems of characteristics for eqs 21a and 21b coincide. Now exclude the parameter s from the second and the last of eq 33 and use eq 23 to obtain
dC1 + g(z;C0)C1 ) h(z;C0) dz
(35)
{
( ) ( )
∂C0 1 f ′(C0) u ∂t
C1(z;C0) )
1 B(z;C0)
∫0z dx h(x;C0) B(x;C0)
{
E(C0) C1(z;C0) ) F(C0) ln(B(z;C0)) + C0 B(z;C0)
}
∂2C0 ∂t2
∫0z dx h(x;C0) exp(∫zx dx′ g(x′;C0))
∫0z dx h(x;C0) exp
(
f ′(C0) u
(36)
∫zx dx′
(37)
)
∂C0 (38) ∂t
SinceC1(0;C0) ) C1(0,t0) ) 0, eq 24 is satisfied. The partial derivatives in eqs 36 and 38 (as well as in eq 46) can be obtained by differentiation of eq 29 by t and z provided that C0 ) C0(z,t)
∂2C0 ∂t ∂2C0 2
∂z
)
2
}
(42)
(
)
f(C0) - 1 C0 f ′′(C0) 1+ uτ0 f ′(C0)
where f(C0) is determined by eq 22 and B(z;C0) is determined by eq 40. We consider here only convex upward isotherms for which f ′(C0) < 0 so that the logarithm in the right-hand part of eq 42 is always real. At τ0 f 0, both B and E tend to infinity but their ratio before the braces in the right-hand side of eq 42 is finite. The second term in the braces of eq 42 is also finite when τ0 f 0, but the first term diverges logarithmically. Thus, the right-hand side of eq 42 has a logarithmic singularity as τ0 tends to zero. This peculiarity of eq 42 (and of eq 47) is discussed in section 5. 3.2. Longitudinal Diffusion. Now we assume that τkin ) 0 but D * 0. In this case, the model reduces to one quasilinear PDE
f(C)
∂C ∂2C ∂C +u )D 2 ∂t ∂z ∂z
(43)
C0 ∂C0 ) - /B(z;C0) ∂t τ0
(39a)
C0 f (C0) ∂C0 ) ∂z uτ0B(z;C0)
with the boundary condition given by eq 10 and the initial condition given by the first eq 9. In line with the previous method, we expand its solution in powers of D
(39b)
C ) C0 + DCd + ‚‚‚
C0
)
2
τ0 B(z;C0)
C0 f(C0) 2
z B(z;C0)
f(C0) - 1 [f(C0) - 1]f ′′(C0) ; F(C0) ) - 2; τ0 f ′(C )2
G(C0) )
Equation 35 together with the second equation in eq 24 determines variation of C1 along the characteristic of eq 21 labeled by C0. Its solution is
C1(z;C0) )
G(C0)
0
[f(C0) - 1]
or using the first equation in eq 36
(41)
where B(x;C0) is determined by eq 40. Finally, substitute h(x;C0) from eqs 36 into this equation and use eqs 39a and 39c to integrate it to obtain
E(C0) )
[f(C0) - 1] ∂C0 2 h(z;C0) ) 2f ′(C0) + u ∂t
C1(z;C0) )
(40)
These equations determine the corresponding partial derivatives along a characteristic determined by eq 29. Since C0 is the solution of eq 21a, it as well as its partial derivatives depend both on z and t. However, along the characteristics of that equation, C0 is constant and z and t are bound by eq 29, so that the only variable in the partial derivatives of eq 39 is z. Now substitute eqs 39a in eq 38 and integrate by x′ to obtain
where C0 is a constant parameter (for that reason it is separated from z not by a comma but by a semicolon) and
g(z;C0) )
C0 f ′(C0) uτ0
(uτ0) B(z;C0)
3
[
z
3
[
1+z
]
C02f ′′(C0) uτ0
(39c)
f(C0) + 2C0 f ′(C0) +
]
C02 [f ′′(C0) f(C0) - 2f ′(C0)2] (39d) uτ0
(44)
and obtain two PDEs by the method described above. The first equation is eq 21a for C0, and the second equation obtained by zeroing the coefficient at D is
f(C0)
∂Cd ∂Cd ∂C0 ∂2C0 +u ) -Cd f ′(C0) + 2 ∂t ∂z ∂t ∂z
(45)
The left-hand side and the first term in the right-hand side of this equation coincide with those of eq 21b and 23a. Therefore,
On Inverse Adsorption Chromatography. 1.
J. Phys. Chem. C, Vol. 111, No. 20, 2007 7469
the characteristics of eq 45 coincide with those of eqs 21 and the variation of Cd along the characteristic labeled by C0 is described by eq 41 where h(x;C0)is determined by the second term in the right-hand side of eq 45
Cd(z;C0) )
1 uB(z;C0)
∫0z dx
∂2C0 ∂x2
B(x;C0)
(46)
Substitute eq 39d in eq 46 and integrate to obtain
{
E(C0) Cd(z;C0) ) F(C0) ln[B(z;C0)] + C0 B(z;C0) G(C0) E(C0) )
z B(z;C0)
}
(47)
]
(48)
The dependence u(z) in dimensionless form is (see p 40, eq 6 in ref 20)
u(z) )
uo [P - (P - 1)z] 2
2
1/2
(49)
where P is the ratio of the inlet and outlet pressures and uo is the outlet linear velocity of the carrier gas. The average velocity u is proportional to uo
u ) juo where j )
3 P2 - 1 2 P3 - 1
(50)
(see p. 41, eqs 10 and 11 in ref 20). Equation 48 can be solved by the method of characteristics (see ref 24 and references therein). However, this full solution is difficult to combine with the methods employed above. To make the final results simpler, we consider here an approximate solution on the same level of approximation as accepted above. Assume that the pressure drop over the column is small so that δ ) P - 1 is small and expand the right-hand side of eq 49 in powers of P at P ) 1 preserving only the linear part
u(z) = uo + uo(z - 1)δ
(52)
and substitute in eq 48 to obtain two equations. The first equation again coincides with that of eq 21a and the second equation reads
∂Cδ ∂C0 ∂Cδ + uo ) -Cδ f ′(C0) - u oC 0 ∂t ∂z ∂t uo(z - 1)
with the same functions f(C0) and B(z) as in eq 42. The only difference between eq 42 and eq 47 is that f(C0) - 1 in the former is substituted for f(C0) in the latter. Besides, E(C0) in eq 47 contains an additional factor of 1/u2 in comparison to eq 42 which, however, makes no difference since in our case u ) 1. 3.3. Gradient of the Carrier Gas Velocity. Now we neglect the longitudinal diffusion and the kinetics of adsorption but assume that the linear velocity of the carrier gas u depends on z. The mass balance equation now takes the form
∂C ∂C du f(C) +u ) -C ∂t ∂z dz
C ) C0 + δCδ
f(C0)
f(C0) f(C0) f ′′(C0) f(C0) ; F(C0) ) - 2; G(C0) ) 2 2 uτ0 u τ0 f ′(C0) C0 f ′′(C0) 1+ f ′(C0)
[
Now, in line with the above approaches, expand C in powers of δ
(51)
Thus, according to this equation which should be valid at sufficiently small δ, the velocity of the mobile phase increases linearly from the inlet to the outlet of the column. One can see from Figure 2 in ref 20 that this approximation gives a semiquantitatively correct result already at δ e 0.5.
∂C0 (53) ∂z
Characteristics of this equation again coincide with those of eq 21 and the variation of Cδ along the characteristic labeled by C0 (eq 29) again obeys eq 35 with the same value of g(z;C0) as for C1 andCd. Thus, the value of Cδ can be obtained from eq 41 which now takes the form
uoC0 Cδ(z;C0) ) B(z;C0)
[
∫0z dx 1 + x C-0 1
]
∂C0 B(x;C0) (54) ∂x
The right-hand side of this equation can be explicitly integrated using eqs 39b, 40, and 50 to obtain
[
Cδ(z;C0) u f(C0) z z ) (f(C0) C0 2uτ0 B(z;C0) j uτ0
]
C0 f ′(C0)) - 1 (55) 4. Solution of the Inverse Problem: Isotherm of Adsorption In general, one should consider C as being a function of all three parameterssτkin, D, and δssimultaneously. Expanding it in these parameters, one obtains
C ) C0 + τkin
∂C ∂C ∂C +D +δ + ‚‚‚ ∂τkin ∂D ∂δ
(56)
Here the partial derivatives are taken assuming that only one parameter is varied, the other two being zero. Thus in eq 56, ∂C/∂τkin is C1 eq 42, ∂C/∂D is Cd eq 47, and ∂C/∂δ is Cδ eq 55 and eq 56 takes the form
C ) C0 + τkinC1 + DCd + δCδ + ‚‚‚
(57)
Designate the experimental EP as t ) texp(C). The characteristic labeled by C0 intersects the line z ) 1 at the time determined by eq 29. Thus, the EP can be determined by the equation (eq 29 with z ) 1, u ) 1, t ) texp(C))
f(C0) ) texp(C0 + τkinC1 + DCd + δCδ) + τ0 ln(C0) (58) where C1, Cd, and Cδ are determined by eqs 42, 47, and 55, respectively, and depend on C0, f(C0), and the derivatives of f(C0) up to the second order. Given texp(C), eq 58 is, in fact, a nonlinear ODE of the second order for the unknown function f(C0). To select one of its solutions, one has to determine the values of f(C0) and f ′(C0) for some specific value of C0. We do not follow this route here. Instead we use the method of iterations suggested by the form of eq 58: First nullify τkin, D, and δ to make the right-hand side of eq 58 a known function of C0 (experimental EP) and obtain the first approximation for f(C0) which is basically the
7470 J. Phys. Chem. C, Vol. 111, No. 20, 2007
Bakaev
Figure 1. Corrections to ECP. τ0 ) 0.1: diamonds - C1/C0 eq 42; squares - Cd/C0 eq 47; triangles - Cδ/C0 eq 55. τ0 ) 0.01: solid line - C1/C0 eq 42.
ECP approximation. Then for each C0, one can use eqs 42, 47, and 55 to determine C1, Cd, and Cδ and to find the next approximation of f(C0) from eq 58 again using the experimental EP now shifted along the concentration axis by C1, Cd, and Cδ multiplied by τkin, D, and δ, respectively. This procedure can be repeated with the refined values of f(C0), etc. The experimental variable C0 is usually not a continuous variable with respect to time but a n-dimensional vector C B0 ) {C10,C20,...,Cn0} because some chromatographic detectors measure solute concentrations at discrete moments of time. Therefore f(C0) is also a vector whose i component f i is f(Ci0). The algorithm of the iterative solution of eq 58 described above can be symbolically written down as
Bf k+1 ) Aˆ Bf k
(59)
Here nonlinear operator Aˆ expresses the fact that the right-hand side of eq 58 depends not only on the value of C0 but also on the function f(C0). Equation 59 is a particular case of the Richardson iterations with the acceleration parameter λk ) 1
Bf k+1 ) Bf k + λk(Aˆ Bf k - Bf k)
(60)
The problem of convergence for this equation is discussed in the next section for a particular isotherm of adsorption. There is still another way to solve the inverse problem.5 If an adsorption isotherm can be expressed by an equation with adjustable parameters, such as for instance eq 61, one can try to obtain these parameters by minimizing difference between the left and right-hand sides of eq 58. 5. Discussion The experimental application of the above results is considered in detail in the accompanying paper II (see ref 25). Here we only discuss some peculiarities of the method. Consider an isotherm of adsorption described by the following equation:
N0 ) exp[a0 + a1 ln C0 + a2 (ln C0)2]
(61)
This is eq 8 from part II (see ref 25) presented in dimensionless variables of the present paper where the coefficients ak with k
> 2 are truncated. The values of the parameters a0, a1, and a2 correspond to Figure 3 of part II (pyridine on the silica surface at 100 °C). The values of f(C0) and its derivatives that enter eqs 42, 47, and 55 were calculated using eqs 22 and 61. The left-hand sides of eqs 42, 47, and 55 are presented in Figure 1. First, it is seen from Figure 1 that the corrections to ECP due to the adsorption kinetics (C1/C0) and that due to the longitudinal diffusion (Cd/C0) at τ0 ) 0.1 are close to one another. This means that one can approximately substitute Cd/ C0 for C1/C0 and use an effective dispersion factor (τkin + D) instead of D in eq 58. Thus, the correction due to adsorption kinetics can be approximately expressed through that of longitudinal diffusion and vice versa. Such an additivity of spreading factors is well-known in linear chromatography.9,11,20 As far as eqs 42 and 47 are close to one another, such an additivity is approximately observed in our model too. Second, the values of C1/C0 and Cd/C0 depend on τ0 and even become logarithmically divergent as τ0 tends to naught (see discussion of eq 42). This is illustrated in Figure 1 where τ0 ) 0.1 corresponds to the value obtained in our experiments in part II, whereas the smaller values of τ0 in Figure 1 are chosen just for comparison. One can see that when the value of τ0 decreases 10 times, the value of C1/C0 increases by about 50%. Besides, the domain of concentrations for the plot of C1/C0 in Figure 1 shrinks considerably when τ0 decreases. This is because τ0 equals an amount of solute (in the units of V/mC/max) injected in the column (cf. eqs 4b, 6, and 10a or integrate eq 10). The smaller the amount that is injected, the smaller maximum concentration of EP at the outlet of the column one obtains. The initial concentration of solute at the inlet of the column always equals unity in our units because its dimensional counterpart C/max is chosen as the unit of concentrations. However, the value of C/max is not limited in the model and can tend to infinity if the injected amount is finite but the time of injection determined by τ0 tends to naught. This, however, cannot happen in a real gas chromatograph because as soon as C/max exceeds the concentration of the saturated vapor at the column temperature, the gaseous solute condenses and the boundary condition at the outlet of the column drastically changes. Thus, the above-mentioned logarithmic singularity in
On Inverse Adsorption Chromatography. 1.
J. Phys. Chem. C, Vol. 111, No. 20, 2007 7471
Figure 2. Convergence of iterations in eq 59. Solid line f(C0); squares, first approximation for f(C0); diamonds, second approximation for f(C0).
eqs 42 and 47 is not a real artifact of the model because the model itself is valid only at sufficiently large values of τ0. Now we discuss the convergence of iterations described by eqs 59 and 60. First we use eqs 22 and 61 with parameters mentioned above to calculate f(C0). The plot of this function is shown in Figure 2 as a line. Then we use eqs 58, 42, 47, and 55 to calculate texp(C0) using the following values of parameters: τkin ) 0.01; D ) 0.02; δ ) 0.02. These parameters approximately correspond to those used in the experimental part of this work; see part II. We take this texp(C0) as a simulated EP and try to calculate from it f(C0) by the methods described in sections 3 and 4. Consider eq 59 (eq 60 with λk ) 1). Take the value of texp(C0) + τ0 ln(C0) as the first approximation for f(C0). (The approximation is very close to the ECP which is f(C0) = texp(C0).) This approximation is designated by squares in Figure 2 (the first approximation for f(C0) is designated as Bf1 in eq 59). The second approximation for f(C0) is designated as Bf2 in eq 59. To obtain Bf2 from Bf1, one has to numerically differentiate texp(C0) + τ0ln(C0) and to use eqs 58, 42, 47, and 55. After that, Bf3 can be obtained from Bf2 in the same manner, etc. The average deviation of Bfk+1 from Bfk is
|Bf k+1 - Bf k| )
1
n
∑
fk(Ci0)
- 1| | n i)1 f (Ci ) k+1 0
(62)
where fk(Ci0) is the value of the kth approximation for f(C0) at the ith value of C0. For k ) 1, this deviation is about 19% (deviation of squares in Figure 2 from diamonds), but for k ) 2, this deviation is only 0.3%. The third approximation for f(C0) almost coincides with its second one. The second approximation for f(C0) is shown in Figure 2 as diamonds and almost coincides with the true value of f(C0). In this example, one knows in advance the real value of f(C0), and one can see that when two successive iterations of that function are close to each other they are close to the real value of the function (at least in this instance). The ECP approximation of EP (the right-hand side of eq 30) depends only on the derivative of an adsorption isotherm. Moreover, it depends only on adsorption in the domain of concentrations which appear at the outlet of the column. The
maximum concentration at the outlet of the column is much smaller than unity: the maximum concentration at the inlet of the column. The former concentration depends on the length of the column. For example, in the example in Figure 2, the maximum concentration of EP is about 5 × 10-3 which is 200 times less than the maximum concentration at the inlet of the column. However, if one could measure EP inside the column (at z < 1) or use shorter columns, one would prolong the isotherm in the domain of higher concentrations. Thus, the chromatographic processes inside the column depend on the isotherm of adsorption in the domain of higher concentrations than those which appear at the outlet of the column. For this reason, an attempt to refine the ECP method by computer simulation of longitudinal diffusion required extrapolation of the isotherm obtained from ECP approximation in the domain of higher concentrations.19 In the modified ECP method considered above, such an extrapolation is implicitly included through the higher derivatives of the adsorption isotherm which are required to calculate corrections to C0 in eqs 42 and 47. For if one knows the value of adsorption together with its first, second, and third derivatives at some concentration (as is the case in eqs 42 and 47), one can extrapolate the corresponding isotherm to higher concentrations. Thus, although the original ECP approximation considers only the situation at the outlet of a chromatographic column, the modified ECP considered in this paper takes account of processes that occur inside the column up to the injector. Finally, it should be mentioned that extrapolation of an adsorption isotherm to higher concentrations (pressures) is reliable only when one knows in advance the equation of the isotherm. This is the case with adsorption on heterogeneous surfaces. In this case, there are only a small number of empirical equations such as eq 61 or its modifications which describe experimental data.10 6. Conclusion The purpose of this paper is to refine the ideal chromatography model which is the foundation of the ECP method of inverse chromatography. (The ideal adsorption chromatography model neglects the kinetics of adsorption and other nonideality factors which our new model takes into account.) The ECP
7472 J. Phys. Chem. C, Vol. 111, No. 20, 2007 method is applied further to the study of the energy distributions of strong adsorption sites on heterogeneous surfaces of silicate glasses (glass fibers) and such nanoparticles as, for example, fumed silica (see the accompanying paper II, ref 25). This determines the peculiarity of the chromatography model accepted in this communication. First, isotherms of adsorption on strong adsorption sites of a heterogeneous surface are convex upward (have a negative second derivative) which is a necessary precondition for the applicability of modified ECP method developed in this paper. Second, the peculiarity of those isotherms allows one to consider the nonlinear chromatography and still assume that the concentration of a one component solute is small. Third, the influence of such nonideality factors as longitudinal diffusion and kinetics of adsorption on elution profiles is also assumed to be small in comparison with that of isotherm nonlinearity. The corrections for the ideal chromatography model are obtained for two sources of nonideality: the kinetics of adsorption and the longitudinal diffusion. These corrections are obtained in an analytical form (see eqs 42 and 47) as the linear terms of the expansion of the model solution into powers of the kinetics of adsorption and longitudinal diffusion coefficients. For a numerical illustration considered in this paper (Figure 1), these corrections are almost additive (as in the linear chromatography). (The correction for the carrier gas compressibility is also obtained by the same method, see eq 55, although its more accurate versions are published in the literature.) The nonideality corrections mentioned above depend on the first, second, and third derivatives of an adsorption isotherm and on the constant that determines the length of an injection pulse. An isotherm of adsorption derivative (which is the first step for calculation of energy distribution) can be obtained from an elution profile by this modified ECP method by two ways: One can either solve the second order ordinary differential equation or use the iteration method. Only the latter way is tested here, and the convergence of iterations is shown to be excellent for a numerical example considered in Figure 2 (however, in some other cases iterations do not converge; see part II, ref 25).
Bakaev References and Notes (1) Kiselev, A. V.; Yashin, Ya. I. Gas-Adsorption Chromatography; Plenum Press: New York, 1969; Chapter 4. (2) Glueckauf, E. J. Chem. Soc. (London) 1947, 1302. (3) Conder, J. R.; Young, C. L. Physicochemical Measurement by Gas Chromatography; Wiley: Chichester, U.K., 1979; Chapter 9. (4) Cremer, E.; Huber, H. Angew. Chem. 1961, 73, 461. (5) Guiochon, G.; Felinger, A.; Shirazi, D. G.; Katti, A. M. Fundamentals of PreparatiVe and Nonlinear Chromatography, 2nd ed.; Elsevier: Amsterdam, 2006. (6) Conder, J. R.; Purnell, J. H. Trans. Faraday Soc. 1969, 65 (3), 824. (7) Roles, J.; Guiochon, G. J. Chromatogr. 1992, 591, 245. (8) Bakaev, V. A.; Bakaeva, T. I.; Pantano, C. G. J. Phys. Chem. B 2002, 106, 12231. (9) Guiochon, G.; Lin, B. Modeling for PreparatiVe Chromatography; Academic Press: New York, 2003. (10) Rudzinski, W.; Everett, D. H. Adsorption of Gases on Heterogeneous Surfaces; Academic Press: London, 1992; Chapter 2. (11) Golshan-Shirazi, S.; Guiochon, G. J. Chromatogr. 1992, 603, 1. (12) Glueckauf, E.; Coates, J. I. J. Chem. Soc. (London) 1947, 1315. (13) Zabezhinski, Ya. L. Zhurnal fizicheskoi khimii 1943, 17, 32 (in Russian). (14) Tikhonov, A. N.; Zhukhovitskii, A. A.; Zabezhinski, Ya. L. Zhurnal fizicheskoi khimii 1946, 20, 1113 (in Russian). (15) Tychonov, A. N.; Samarski, A. A. Partial Differential Equations of Mathematical Physics; Holden-Day Inc.: San Francisco, 1964; Vol. 1, p 140 (translation from Russian). (16) Rhee, H.-K.; Aris, R.; Amundson, N. R. First-Order Partial Differential Equations. Theory and Application of Single Equations; Prentice Hall: Englewood Cliffs, NJ, 1986; Vol. 1, Chapters 5 and 6. (17) Shkilev, V. P. Russ. J. Phys. Chem. 1999, 73, 633. (18) Lin, B.; Golshan-Shirazi, S.; Guiochon, G. J. Phys. Chem. 1989, 93, 3363. (19) Shkilev, V. P.; Bogillo, V. I. Russ. J. Phys. Chem. 1998, 72, 1502. (20) Guiochon, G.; Guillemin, C. L. QuantitatiVe gas chromatography (for laboratory analyses and on-line process control); Elesevier: Amsterdam, 1988. (21) Ames, W. F. Numerical Methods for Partial Differential Equations, 3rd ed.; Academic Press: Boston, 1992; p 252. (22) L. F. Shampine, Numerical Solution of Ordinary Differential Equations; Chapman & Hall: New York, 1994; p 392. (23) Goldstein, S.; Murray, J. D. Proc. R. Soc. London, A 1959, 252, 360. (24) Shkilev, V. P. Russ. J. Phys. Chem. 1999, 73, 1281. (25) Bakaev, V. A.; Bakaeva, T. I.; Pantano, C. G. J. Phys. Chem. C 2007, 111, 7473-7486.
