Inverse Reaction Chromatography. 1. Computer Simulation of

Jul 13, 2009 - Computer Simulation of Hydrogen/Deuterium Exchange on Oxide Surfaces ... Materials Research Institute, The Pennsylvania State Universit...
0 downloads 0 Views 189KB Size
Download PDF
13886

J. Phys. Chem. C 2009, 113, 13886–13893

Inverse Reaction Chromatography. 1. Computer Simulation of Hydrogen/Deuterium Exchange on Oxide Surfaces V. A. Bakaev Materials Research Institute, The PennsylVania State UniVersity, UniVersity Park, PennsylVania ReceiVed: December 12, 2008; ReVised Manuscript ReceiVed: April 17, 2009

This Article concerns the mathematical model of reaction chromatography. It considers a chromatographic column whose stationary phase, e.g., an oxide, contains OH groups on its surface. The mobile phase contains a solute, e.g., D2O, whose molecules can exchange deuterium for hydrogen atoms on the surface (H/D exchange). The system of partial differential equations describing the chromatographic process with H/D exchange in such a column and the method of characterizing their numerical solution are developed. The concentrations of HOH, HOD, and DOD at the outlet of the column vs time are obtained for frontal and pulse injections of D2O. These outlet curves show how one can determine the concentration of OH groups on the surface of an oxide and evaluate the rate constant for the H/D exchange. The method is named inverse reactive chromatography (IRC) because its purpose is to determine not the property of a mobile phase (its composition or reactivity) as in the conventional reaction chromatography but that of a stationary phase. The comparison of the computer simulations described in this paper with the experiment is presented in the companion paper (Bakaev, V. A.; Pantano, C. G. J. Phys. Chem. C 2009, the same issue). 1. Introduction Chromatography is basically a method for analysis (analytic chromatography) and separation (preparative chromatography) of gaseous and liquid mixtures. Since the early days of its application, the chemical reactions between the mixture components have plagued the method (see Introduction1). Later, chromatography found some application as a method for studying those chemical reactions (see a review of early papers2). Currently, chromatography is widely used in catalytic research as an analytic instrument coupled with a separate chemical reactor, but there are also applications where a chromatographic column is simultaneously used as a chemical reactor.3 Such an application of chromatography is called reaction (dynamic3) chromatography. Chemical reactions between the components of an analyzed mixture (solute) influence the shape of elution profiles (breakthrough curves for frontal analysis). The basic method for understanding elution profiles is the solution of the mass conservation (mass balance) partial differential equations (PDEs), which describe advection and adsorption (this Article is concerned only with adsorption chromatography) in chromatographic columns.4,5 The reaction chromatography is described by basically the same PDEs; the only difference is the source term that describes chemical reactions in the mixture. One can use the Laplace transform to easily solve those PDEs when the latter are linear.2 This situation occurs in the analytic or reaction chromatography when adsorption isotherms are linear, which occurs on a homogeneous surface at infinitesimal concentrations of the solute components. The nonlinear reaction chromatography was considered by the method of perturbation when deviations from linearity are small.1,4 Less rigorous methods based on the plate theory or the stochastic theory of adsorption have also been used in the reaction and dynamic chromatography.6 Apart from the analysis of mixtures, adsorption chromatography is also used for study of adsorbents. This application of

chromatography is usually called inverse chromatography, in particular, inverse gas chromatography (IGC) (see refs 7-9 and references therein). One of the IGC applications is measurement of adsorption isotherms and heats of adsorption.8 There is the comprehensive review of the methods for measuring adsorption isotherm in liquid chromatography (which, however, are not classified as inverse chromatography).5 The term inverse chromatography is used here to designate a method which employs the chromatographic technique to determine surface properties of an adsorbent, such as, e.g., energy distribution of adsorption sites8,9 or concentration and reactivity of active sites as described below. Inverse adsorption chromatography deals only with physical adsorption. In particular, IGC can determine not only isotherms and differential heats of physical adsorption but also the energy and entropy distribution of adsorption sites on a surface.8,9 Inverse reaction chromatography (IRC) considered in this paper relates to the reaction chromatography mentioned above as IGC relates to the conventional gas chromatography. In particular, we consider here the hydrogen-deuterium (H/D) exchange between hydroxyl groups on the surface of oxide and water isotopes (HOH, HOD, DOD) in the gas phase. The H/D exchange is the simplest case of the surface chemical reaction. Thus, the distinction of IRC from the reaction chromatography considered in refs 1, 2, 4 is in that the latter deals with chemical reactions in the mobile phase whereas the former deals with chemical reactions in the stationary phase. This distinction is reflected in the structure of the mathematical model of IRC considered in this Article. We consider here a chromatographic column packed with an oxide contained on its surface hydroxyl groups. Suppose one maintains at the inlet of the column a constant concentration of DOD (solute) and constant flow of the carrier gas through the column. It is known that at certain conditions (cf. discussion of eq 22) this concentration of solute will spread along the column with constant velocity and a sharp front (shock). Imagine an observer at the shock who moves with it along the column. He

10.1021/jp810982w CCC: $40.75  2009 American Chemical Society Published on Web 07/13/2009

Inverse Reaction Chromatography

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13887

sees before him fresh portions of hydroxyl groups on the surface of oxide, which take part in the H/D exchange with solute until all the DOD molecules of solute at the moving shock will convert into HOH molecules. These HOH molecules will move along the column intact, but another portion of solute behind the shock will continue H/D exchange with hydroxyl groups on the surface. The model described below is, in fact, just a translation of this qualitative picture into a mathematical language. The central assumption of the model is eq 17, which actually excludes chromatographic separation of the solute components (HOH, HOD, and DOD). Thus, although the model uses the conventional mathematical means of chromatography, mass balance equations eqs 11, adsorption isotherms (identical to all components of solute) eqs 19 and 34, and so forth, it is not a conventional chromatographic model: It can also be called a dynamic H/D exchange. Equation 17 allows one to convert the conventional equations of chemical kinetics eqs 3-6 and the mass balance eqs 11 into eqs 24-27. These are a combination of PDEs and ordinary differential equations (ODEs) written down in the laboratory coordinate system. Then, we transform them into the coordinate system moving with the shock velocity along the column as in the qualitative picture described above. This transformation converts PDEs eqs 24 into ODEs eqs 29. These ODEs are then numerically solved together with eq 25 by the predictor-corrector algorithm which is frequently employed for the numerical solution of the ODEs systems. This is the physical meaning of the algorithms described in section 2.2. The mathematics of a first-order PDE and a system of such PDEs is considered in refs 10, 11 with illustrations taken mainly from the chromatography and chemical engineering. These illustrations are limited mainly to the ideal chromatography model,5 which neglects adsorption kinetics and axial dispersion because inclusion of the latter would increase the order of PDEs. More complete models of chromatography, their physicochemical substantiations, and computer simulations are reviewed in ref 5. The model presented below is based on the ideal chromatography model described in refs 5, 10. It is essentially nonlinear as is the water isotherm on oxides and borrows from refs 5, 10 only the theory of shock formation. All the other features of the model described in this paper are different from what is described in refs 5, 10. We assume that only adsorbed molecules react with active sites on the surface and there is no reaction between molecules themselves in either mobile or stationary phase. Under these assumptions, we determine elution profiles (or rather the breakthrough curves, since our border condition correspond to frontal analysis of chromatography) presented in Figure 3. It should be emphasized here that the physical meaning of the plots presented in Figure 3 is different from the elution profiles (breakthrough curves) of genuine chromatography because, as mentioned above, chromatographic separation is excluded from our model. Their meaning is briefly discussed in section 4.

tSiOH + DOD ) tSiOD + HOD

(1)

tSiOH + HOD ) tSiOD + HOH

(2)

2. Description of Model 2.1. Proton-Deuterium Exchange on Oxide Surfaces. Consider a chromatographic column packed with an oxide, e.g., silica, through which an inert carrier gas, e.g., helium, containing deuterium oxide (DOD) as a solute flows. The surface of silica contains silanol groups (tSiOH), which can exchange hydrogen for deuterium creating tSiOD groups on the surface and DOH and HOH (water) molecules in the carrier gas. This exchange is described by the reversible reactions:

We neglect the exchange between silanol groups and gasphase molecules because usually the concentration of the latter is considerably smaller than that of adsorbed molecules as well as direct exchange between adsorbed molecules (see the end of section 3). In what follows, we designate concentrations of HOH, HOD, and DOD in the solute by the subscripts 1, 2, and 3, respectively. The concentrations of H- and D-silanol groups (tSiOH and tSiOD) are designated by the subscripts H and D, respectively. Designated by Qi (i ) 1, 2, 3), QH, and QD are the concentrations in the stationary phase of the corresponding molecules, silanol groups, and deuterium exchanged silanol groups, respectively. (If F is the oxide bulk density, Σ its specific surface area, and Γi adsorption of the ith molecules per the unit of surface, then Qi ) FΣΓi.) The kinetic equations corresponding to the reactions of eqs 1 and 2 can be written as

∂QH ) k[(QD(Q1 + Q2) - QH(Q3 + Q2)] ∂t

(3)

δQ1 ) k[QHQ2 - QDQ1] δt

(4)

δQ2 ) k[QDQ1 + QHQ3 - QHQ2 - QDQ2] δt

(5)

δQ3 ) k[QDQ2 - QHQ3] δt

(6)

These are the conventional second-order chemical kinetics equations for bimolecular reactions of H/D exchange (t is time). The rate constant k refers to the H/D exchange between a silanol group and an adsorbed molecule. These rate constants are assumed to be the same for the direct and inverse reactions. Since the ratio of the rate constants for the direct and inverse reactions equals the equilibrium constant, we actually assume that the equilibrium constant for each of the reactions of such as eqs 1 and 2 is unity. The latter means that we neglect the difference in vibration energies and entropies between the reactants and the products of reactions eqs 1 and 2 which arise due to the mass difference between H and D. (In other words, we neglect both kinetics and equilibrium isotope effects.) The values of Qi, QH, and QD are proportional to the surface concentrations of the corresponding components. All these values depend on time and the coordinate z along the column. The negative terms in the brackets on the right-hand side of eq 3 correspond to the reaction on the left-hand sides of eqs 1 and 2 which decrease the concentration of H-silanol groups, and the positive term corresponds to the inverse reactions which increase this concentration. The derivatives on the left-hand sides of eqs 4-6 have nonstandard notations to emphasize that they refer not to the full variations of Qi with time (see below) but only to those their parts which correspond to the reactions of eqs 1 and 2. Equations 3-6 are not independent:

13888

J. Phys. Chem. C, Vol. 113, No. 31, 2009

Bakaev

δQ1 δQ2 δQ3 + + )0 δt δt δt

(7)

∂QH δQ3 δQ1 ) ∂t δt δt

(8)

(9)

(10)

The values of Qi(z, t) and Ci(z, t) obey the system of mass conservation (balance) PDEs:

ε

∂Ci ∂Qi ∂uCi δQi + (1 - ε) ) -ε + (1 - ε) ∂t ∂t ∂z δt

(11) Here, ε is the porosity of column packing (fraction of the column volume occupied by gas), (1 - ε)is the fraction of the column volume occupied by the stationary phase (oxide), and Ci is the concentration of the ith molecule (ith component of the solute) in the mobile phase. The first term on the righthand side of eq 11 describes the accumulation of the ith component of the solute delivered by the carrier gas into an infinitesimally thin slice of the column, and the second term describes the addition to it due to the chemical reactions described by eqs 4-6. The terms on the left-hand side of eq 11 describe the total rate of change of that component in the slice. Equation 11 is a conventional equation of reaction chromatography;4,10 however, in our case it is not complete. One can see that eqs 4-6, which constitute the right-hand side of eq 11, depend on QH and QD, which obey the separate eq 3. This is the peculiarity of the mathematical model of inverse reaction chromatography that distinguishes it from conventional and reactive chromatography. The additional differential equation eq 3 plays the central role in this model describing the change of the surface in the process of H/D exchange. In contrast to that, the surface does not change in reaction chromatography, although it can catalyze reactions occurring in the mobile phase.4,10 Equations 11 and 3 should be supplemented by the initial and boundary conditions, which are chosen here as

Qi(z, 0) ) 0;

Ci(z, 0) ) 0;

C1(0, t) ) C2(0, t) ) 0;

QD(z, 0) ) 0

C3(0, t) ) φ(t)

C(z, 0) ) 0

C(0, t) ) φ(t)

Differentiate this equation by t and use eq 7 to obtain eq 8, which means that the latter is the proton conservation equation. In addition, the summary number of H- and D-silanol groups equals the number of active sites. Thus, if Qtot is the concentration of the active sites in the stationary phase, then

QH + QD ) Qtot

∂C ∂Q ∂uC + (1 - ε) +ε )0 ∂t ∂t ∂z Q(z, 0) ) 0;

The above reactions do not change the total concentration of protons on the surface which can be written as

2Q1 + Q2 + QH ) const

ε

(12) (13)

These mean that initially one has a clean (empty of adsorbate) column with only H-silanol groups on the oxide surface and the concentration of pure deuterium oxide at the inlet of the column is maintained at the level of φ(t). Add eqs 11-13 together and use eq 7 to obtain

(14) (15) (16)

In these equations, C and Q designate the summary concentration of all components in the mobile and stationary phases, respectively. These equations look like conventional equations of chromatography for a one-component solute,4,5,10 but, in fact, they have another meaning, since C and Q are not concentrations of individual components. Now, assume that in eq 11

Ci ) Cxi ;

Qi ) Qxi ;

i ) 1, 2, 3

x1 + x2 + x3 ) 1

(17) (18)

This is the central assumption of the model. It means that the mole fractions of the ith components xi are the same in the mobile and stationary phases. In other words, individual isotherms of adsorption for all components are identical. This assumption is in accord with the neglect of the isotope effects in eqs 3-6 and a reasonable approximation for adsorption of isotopes as one can see, e.g., from Figure 1,12 which shows that adsorption isotherms of water (HOH) and deuterium oxide (DOD) on the silica surface are almost identical. Equations 17 and 18 exclude chromatographic separation of components (HOH, HOD, and DOD) from our model. They transform the meaning of eqs 14-16, converting them into real equations of one-component chromatography which are independent of eqs 3-6 describing H/D exchange. Equations 14-16 can be easily solved if one connects C and Qby an equilibrium adsorption isotherm equation

Q ) Q(C)

(19)

reducing eqs 14-16 and 19 to the ideal chromatography model.5 Neglect also the compressibility of the carrier gas: du/dz ) 0. The latter approximation is valid for liquid chromatography and is often also made also in gas chromatography,10 with the necessary correction being made at the final stage of calculations (see, e.g., ref 7). With these assumptions, eq 14 reduces to the basic equation of the ideal chromatography model4,5,10

∂C ∂C +u )0 (1 + 1 -ε ε dQ dC ) ∂t ∂z

(20)

It is known that eq 20 with the initial and boundary conditions of eqs 15 and 16 and a convex upward isotherm can have a discontinuous solution (shock).4,5,10 When the shock appears (which depends on the exact shape of the border condition eq 16), its propagation speed is determined by the following solute conservation equation

εuCA dt ) [εC + (1 - ε)Q(C)]A dz

(21)

The left-hand side of the equation is the flux of solute over the border of the shock, and the right-hand side is the quantity

Inverse Reaction Chromatography

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13889 2.2. Numerical Solution by the Method of Characteristics. Equations 25-27 are conventional equations of chemical kinetics similar to eqs 3-6, but eq 24 has the appearance of a mass balance equation of reaction chromatography. However, the first coefficient in the left-hand side of eq 24 coincides with that in the right-hand side of eq 22. This means that characteristics of eq 24 are determined by the same eq 22 as the shock pass, with the difference only being the initial condition: For the shock pass, the initial condition is t(0) ) 0, while for the characteristics of eqs 24, the initial condition is t(0) ) t0 where t0 is a parameter determining (labeling) the characteristic. In what follows, we call it a t0-characteristic. Along a t0-characteristic, the left-hand side of eq 24 is proportional to a full derivative

dxi ∂xi ∂xi ) + ush dt1 ∂t ∂z

(28)

where ush ) dz/dt as determined by eq 22. Using eq 28, one can rewrite eqs 24, 26, and 27 as Figure 1. Calculation along characteristics.

of solute which is required to fill the volume A dz of the bed (A is the cross section of the bed). This gives the so-called jump condition for the particular case of propagation of a shock in an empty of adsorbate column

dt 1 - ε Q(C) 1 ) 1+ dz u ε C(z, t)

[

]

(23)

which means that the total concentration of adsorbate does not change in the surface chemical reactions. Substitute eq 17 into eq 11 and use eqs 14 and 23 to obtain

δxi

t) 1 - ε Q(z, t) +u ) (1 + 1 -ε ε Q(z, ) C(z, t) ∂t ∂z ε C(z, t) δt ∂xi

∂xi

(24) Equations 3, 4, and 6 reduce to

∂QH ) kQ[QD(x1 + x2) - QH(x2 + x3)] ∂t

(29a)

dx3 ) λk(QDx2 - QHx3) dt1

(29b)

(22)

(For other derivations of the jump condition in a more general form, see in ref 10, pp 213, 333, 344, 363, 395.) Equation 22 describes the trajectory of a shock-shock pass.10 Insert eq 17 into eq 7 and use eq 18 to obtain

δQ/δt ) 0

dx1 ) λk(QHx2 - QDx1) dt1

(25)

δx1 ) k[QHx2 - QDx1] δt

(26)

δx3 ) k[QDx2 - QHx3] δt

(27)

Thus, assumptions of eqs 17 and 18 split our model into two systems of equations: Equation of ideal chromatography model (eq 20) and equations of isotope exchange 25-27.

λ)

(1 - ε)Q εC + (1 - ε)Q

(29c)

Here, t1 ) t(z; t0)is the equation of a t0-characteristic that is the solution of eq 22 with the initial condition t1(0) ) t0. These equations should be solved together with eq 25 where the derivative in the right-hand side is taken along the line parallel to the t-axis (z ) z0). We call this line a z0-characteristic. Now, eq 29, which together with eq 25 and eq 18 determines the isotope exchange in our model, looks like a conventional chemical kinetics equation. As mentioned in the Introduction, this is achieved by transfer from the laboratory coordinate system to the coordinate system moving with the shock velocity. Their physical meaning becomes especially clear when the shock velocity is constant. This occurs when the isotherm of adsorption is convex upward and the concentration of the third component (DOD in our case) is maintained constant at the inlet of the column, which means that eq 16 is

φ(t) ) 0 if t < 0

and

φ(t) ) C03 if t g 0

(30)

Equations 14-16 with this shape of φ(t) are called the Riemann problem, the saturation of clean column problem (see ref 10, p 207), or the frontal analysis in chromatography (the elution profile is then called the breakthrough curve).4 The solution of this problem is quite simple (see ref 10, p 221): The rectangular plug of solute concentration moves along the column with constant velocity

/(

ush ) u 1 +

(1 - ε)Q(C03) εC03

)

(31)

13890

J. Phys. Chem. C, Vol. 113, No. 31, 2009

Bakaev

Thus, coefficients Q in eq 25 and λ in eq 29c are constants, and eq 29 differs from eq 24 in that the latter is written down in the laboratory coordinate system whereas the former is written down in the coordinate system moving with the constant velocity ush. Equation 29 shows that if one considers a portion of solute in a coordinate system moving along the column with the constant velocity ush, then the composition of that portion changes only due to its interaction with active sites on the surface, with the exchange of components between the neighboring portions of solute being excluded. In particular, along the shock pass z ) usht which is also the 0-characteristic (t0 ) 0), the moving front sees the clean column where QD ) 0 and QH ) Qtot eq 10. Thus, eq 29 with the initial conditions of x1(0) ) 0 and x3(0) ) 1 has a simple solution along the 0-characteristic (shock pass).

x1(t1) ) 1 - (1 + t1λkQtot) exp(-λkQtott1); x3(t1) ) exp(-λkQtott1)

(32)

The equation of the t0-characteristic is z ) ush(t1 - t0) and that of z0-characteristics is z ) z0. All of the t0-characteristics lay above the 0-characteristic, and along all of them, C ) C30, Q ) Q(C30), and λ are constants. To numerically calculate xi(z,t) and QH(z,t) along these characteristics, it is necessary first to discretize the z- and t-axes: z ) k∆z, t ) n∆t, as shown in Figure 1 where t0-characteristics are now labeled by n and z0-characteristics are labeled by k. The point (z, t) becomes in these notations (k, n) and t0- and z0-characteristics become n- and k-characteristics, respectively. Assume that the values of xi(k, n) and QH(k, n) have already been calculated along the n-characteristic as well as along the initial part of the (n + 1)-characteristic up to the point (k, n + 1). Thus, we know those values at the points (k, n + 1) and (k + 1, n) and want to calculate them at the point (k + 1, n + 1) (these points are designated by circles in Figure 1). This can be easily done assuming that the value of QH in eq 29 and those of xi in eq 25 are constants. Then, eqs 25 and 29 are linear ODEs with constant coefficients which have wellknown exact solutions. Those solutions together with xi(k, n + 1) as the initial values for eqs 29 and QH(k + 1, n) as that for eq 25 can be used to obtain xi(k + 1, n + 1) and QH(k + 1, n + 1). The constant QH in eq 29 can be evaluated as an average value between QH(k + 1, n) and QH(k + 1, n + 1) by a prediction-correction method (a version of the midpoint Runge-Kutta13): First, predict QH(k + 1, n + 1) as, e.g., QH(k + 1, n + 1) ) QH(k, n + 1) + QH(k + 1, n) - QH(k, n) and calculate x1(k + 1, n + 1) and x3(k + 1, n + 1) by eqs 29a and 29b with the midpoint value of QH. Then, use the latter values of x1 and x3 to evaluate their midpoint values between the points (k + 1, n) and (k + 1, n + 1) and use eq 25 to correct QH(k + 1, n + 1). Continue these iterations until the maximal relative deviation between succeeding values of QH(k + 1, n + 1), x1(k + 1, n + 1), and x3(k + 1, n + 1) becomes less than a prescribed value. In the computer simulations (see below), that level of convergence was chosen to be 0.0001% and the number of iterations never exceeded 10 (usually that level of convergence was reached in 4-5 iterations). Despite a relatively fast convergence of the above-mentioned iterations, the midpoint method is not very accurate, as shown in the proton conservation test described below (see discussion of Figure 3). A more accurate method assumes a linear variation of QH and xi between the neighboring points. One can find the

slopes of those linear dependencies by the same predictioncorrection method that was described above, but the solution of eqs 29 and 25 with the same initial conditions as above now should be obtained numerically. The adaptive Runge-Kutta subroutine odeint13 was employed to this end. The algorithm described above requires the initial values of QH(0, n) and xi(0, n) on every n-characteristic. These are given by the following equations:

x1(0, n) ) 0;

x3(0, n) ) 1; QH(0, n) ) Qtot exp(-kQn∆t) (33)

It also requires the values of x1(k, 0), x3(k, 0), and QH(k, 0) along the 0-characteristic. The former are given by eq 32. Besides, QH(k, 0) ) Qtot along the 0-characteristic. Finally, the algorithm described above can be applied to more complex border conditions than that of eq 30. In what follows, we consider also the border condition having the shape of a pulse with a tail. In this case, the values of Q in eq 25 and λ in eq 29 are not constants but vary along characteristics, and the characteristic themselves are not straight lines as was the case for the border conditions of eq 30. Nevertheless, the algorithm described above can be applied to this more complex case with only small modifications (see Appendix). 3. Computer Simulations To start computer simulations, one has to determine first of all the adsorption isotherm equation eq 19. We chose the empirical Freundlich isotherm equation (see ref 5, pp 96 and 372) in the following form

Q ) RQm(C/Cs)β

(34)

Here, Cs is the concentration of saturated vapor at the temperature of isotherm; Qm is the adsorbate concentration in the stationary phase at the monolayer coverage of the surface; R and β are adjustable parameters. This form of the Freundlich equation is chosen to diminish the dependence of its parameters on the temperature and to establish a reasonable scale for the coefficient R. The Freundlich equation is sometimes used in the theory of chromatography (see ref 10, p 190, or ref 5, p 372), but the Langmuir equation is used far more often. The reason for that is both analytic and preparative adsorption chromatography methods require the surface of column packing to be as energetically homogeneous as possible. Since the surfaces of the majority of high surface area solids are energetically heterogeneous, they are usually modified to make them more homogeneous and useful for chromatography. Thus, the analytic and preparative chromatography5 actually deals with modified surfaces, whereas inverse chromatography deals with original surfaces of high surface area solids. The Langmuir equation describes adsorption on energetically homogeneous surfaces, whereas the Freundlich equation does this for energetically heterogeneous surfaces. We use the latter because we deal with inverse chromatography in this paper. Besides, it should also be emphasized here that isotherm of adsorptions determines only the values of some constants such as constant λ in eq 29a or constant Q in eq 33 and its shape is unimportant for our model. Another peculiarity of this computer simulation is in the border condition. From the computer simulation point of view, the simplest border condition is the frontal one corresponding

Inverse Reaction Chromatography

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13891

Figure 2. Border condition of eq 35.

Figure 3. Elution profiles for the frontal (symbols) and pulse (lines) boundary conditions.

to eq 30. However, from the experimental point of view, the realization of this border condition requires considerable modification of the inlet of the commercial gas chromatograph (GC) we employed in the experimental part of this work.14 Thus, in this computer simulation we consider a border conditionsthe shape of φ(t) in eq 16sthat can be approximately reproduced with commercial GC without modification of its inlet:

at the outlet of the column packed with an adsorbent whose adsorption isotherm with respect to all isotopes of water (HOH, HOD, and DOD) is described by eq 34 with parameters RQm ) 30.3 mol/m3 and β ) 0.8. (These parameters correspond to the temperature of 100 °C, and the isotherm normalization to the unit of volume of bulk adsorbent.14) For the pulse boundary condition, the injection profile (the shape of φ(t) in eq 35) are shown in Figure 2. For the frontal boundary condition, the partial pressure of DOD is the same as that on the plateau of Figure 2 but tp ) ∞. The partial pressures of HOH and HOD in Figure 2 are zero in both cases. The length of the column L is 22 cm. The other parameters of eq 31 are ε ) 0.507; u ) 2.35 cm/s. These parameters provide the position of the HOH shock in Figure 3 at 25.9 s in agreement with the experiment.14 The accuracy of the computer simulations presented in Figure 3 can be checked with the help of the conservation of protons test. Conservation of protons in the H/D exchange is expressed by eqs 8 and 9. With respect to Figure 3, it means that the total quantity of protons eluted from the column should be equal to Qtotstotal quantity of OH-groups in the column. The EPs can be expressed in units of the flow densities at the outlet of the column as εuC(t) xi(t). In these units, the double area under the HOH plot plus the area under the HOD plot in Figure 3 give the total amount of protons H eluted per unit cross section area of the column. For the complete elution, which is obviously the case for the line plots in Figure 3, those areas equal the total amount of OH-groups per unit cross section area of the column, which is (1 - ε)QtotL (cf. eq 10). For example, the sum of above-mentioned areas in Figure 3 corresponds to 2.88 mol/m2, whereas the quantity of OH groups per unit cross section area of the column is 2.91 mol/m2, the deviation being 1%. However, if the less accurate midpoint Runge-Kutta method of calculation along characteristics mentioned in section 2.2 is used, the deviation increases to 10-30%. Besides, it was also observed in the process of program code development that the programming errors lead to even larger discrepancies between the two above-mentioned quantities. A peculiar feature of Figure 3 is that the HOH, HOD, and DOD plots intersect at one point and the HOD plot has a maximum at that point. The latter fact follows from the first one, which can be explained as follows. One can see from eq 18 that dx2/dt1 ) -dx1/dt1 - dx3/dt1. Add up eqs 29a and 29b and put x1 ) x2 ) x3 to obtain dx2/dt1 ) 0. Since the ordinate of HOD plot in Figure 3 is proportional to x2 and its abscissa linearly depends on t1, this explains why the HOD plot has an extremum at the point of intersection of three plots in Figure 3.

φ(t) ) 0t < 0;

φ(t) ) C030 < t < tp ; φ(t) ) C03 exp[-(t - tp)/τ]t > tp

(35)

This is a broad pulse of D2O with the plateau concentration of C30 and duration of tp with a sharp front border and exponentially decreasing tail as shown in Figure 2 (partial pressure in this figure corresponds to the concentration C03). Such a pulse arises when one maintains the inlet temperature below that of the column.14 This is another peculiarity of the inverse chromatography, since such a regime of a GC inlet is unusual for conventional gas chromatography. The parameters in Figure 2 are tp ) 220 s, τ ) 17 s, and the plateau of partial pressure corresponds to the saturation of water vapor at the temperature 22.7 °C.14 3.2. Elution Profiles (EPs). The EPs (breakthrough curves) at the outlet of the column can be calculated as C(t1)xi(t1) at the time t1 ) t(L; t0). Here, C(t1) and xi(t1) are the summary concentration and mole fraction of the ith component calculated at the end (z ) L) of the t0-characteristic, t1 ) t(z; t0) being the equation of the t0-characteristic. To calculate these elution profiles, one has to organize a loop where t0 varies with a step of ∆t0 and to calculate the t0-characteristic as well as variation of C and xi (i ) 1, 2, 3) along that characteristic starting from their initial values determined either by eq 30 (frontal boundary condition) or by eq 35 (pulse boundary condition). The methods of these calculations are described in section 2.2 and the Appendix. The results of these calculations (computer simulations) are presented in Figure 3. The original EPs of the component concentrations at the outlet of the column are expressed in units of concentration. To facilitate the comparison with experiment, they are renormalized in Figure 3 to the units of the mass selective detector used in the experimental part of this work14 by multiplying by the experimental flow rate of the carrier gas and dividing by the calibration constant (see discussion of Figure 5 in the companion paper14). Figure 3 presents the EPs

13892

J. Phys. Chem. C, Vol. 113, No. 31, 2009

Figure 4. The dependence of HOH elution profile in Figure 3 on the rate constant k in (dm3 mol-1 s-1).

The rate constant k in eqs 29 and 25 was 60 M-1 s-1 (dm3/ (mol s)). The plots in Figure 3 correspond to this rate constant. The variation of the HOH elution profile versus the rate constant k is shown in Figure 4. It is seen that the increase of the rate constant above 60 M-1 s-1 does not significantly change the plot shape. This limits the upper border of the k values, which can be accurately determined from this experiment. The border depends, of course, on the other parameters determining the column characteristics, such as the carrier gas flow rate, and so forth. The model described above neglects the possibility of direct H/D exchange among HOH, HOD, and DOD molecules. The reason for such an assumption is that concentration of these molecules on the surface and a fortiori in the gas phase is smaller than the concentration of hydroxyl groups on the oxide surface when the concentration of DOD at the inlet of the column is sufficiently small. Inclusion of this direct exchange in the model would require introduction at least one additional rate constant. Another assumption is that all hydroxyl groups have the same rate constant of H/D exchange. All these assumptions can be easily eliminated from the model, but such a refinement of the model can be justified only by the comparison with experiment. The experimental part of the work14 qualitatively confirmed Figure 3 but revealed some experimental problems which should be overcome prior to the refinement of the model. 4. Conclusion The mathematical model considered in this paper was called inverse reaction chromatography. It relates to reaction chromatography (RC) by using partial differential equations (PDEs) of mass balance with a source term describing chemical reactions. However, the chemical reactions considered in our model are not those between components of mobile phase, as is the case in conventional RC, but between the components of the mobile phase and the stationary phase. This peculiarity of the model justifies the adjective “inverse” in its title and changes the structure of the RC mathematical model by adding to the system of PDEs describing RC an ODE describing kinetics of chemical reactions between mobile and stationary phases (eq 25). Besides, we considered the simplest chemical reactions of isotope (H/D) exchange that allowed us to make the central assumption of the model (eq 17), which actually removed the

Bakaev chromatographic separation. This assumption allowed us to reduce the system of PDEs to the system of ODEs (eqs 25 and 29), which is much easier to solve numerically. The latter is, however, not a conventional system of ODEs, because the righthand sides of eq 25 and eqs 29 are full derivatives taken in different directions, which we called t0- and z0-characteristics. We employed for numerical solution of this system of ODEs the predictor-corrector algorithm which is often used for the numerical solution of ordinary system of ODEs.13 The final solution is presented in Figure 3. Its correctness is guaranteed by the proton conservation test which allows one to choose a better version of the predictor-corrector algorithm or detect errors in the computer code (see discussion of Figure 3). However, the fulfillment of this test does not guarantee that all assumptions of the model are adequate. This can be checked only by comparison of Figure 3 with experiment. Such a comparison is described in the companion paper.14 Although the experimental part of the work qualitatively confirms Figure 3, it also revealed some experimental problems in the study of H/D exchange with silica which should be overcome in the course of future work. The model considered in this paper can also be considered a dynamic method for studying H/D exchange between water vapor and an oxide surface. In a static version of the method, one brings DOD vapor in contact with an oxide surface and observes equilibration by diffusion. If the kinetics of H/D exchange is fast, slow diffusion would conceal it. The dynamic method considered in this paper facilitates molecular exchange between mobile and stationary phases in the column enabling one to study faster H/D exchange kinetics than with the static method. Acknowledgment. This work was supported with funding from the National Science Foundation, through grant DMR0809657. Appendix Modification of the Method of Characteristic for the Pulse Border Condition. The method described in the section 2.2 should be slightly modified for the border condition of eq 35. The characteristics of eq 20 are straight lines obeying the following equation

t ) t0(C) +

1 1 - ε dQ 1+ z u ε dC

(

)

(1A)

The value of C is constant along a characteristic and determines its slope on the z, t plane, with z being an abscissa. The slopes of all the characteristics are identical for the frontal boundary condition of eq 30, that is, they are parallel straight lines. The same is true for the pulse boundary condition of eq 35 for the t0-characteristics which proceed from the stretch of the t-axis in Figure 2 below the point tp (t0 < tp). However, for the characteristics labeled by the value of t0 larger than tp the value of C (proportional to the partial pressure in Figure 2) exponentially decrease and the values of dQ/dC in eq 1A correspondingly increase (for the parameters of eq 34 used in the paper). One can determine the solution of eq 20 C(z, t) from eqs 35 and 1A

t0(C) ) tp + τ ln(C03 /C)

(2A)

Inverse Reaction Chromatography

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13893

when t0 > tp. In this case, eqs 35 and 1A with fixed values of z and t can be easily solved for C by the Newton-Raphson method taking an initial value of C obtained from eqs 34 and 1A at t0(C) ) tp. The characteristics labeled by these values of t0 are situated in the region of the z, t plane above the straight line determined by the equation

(

)

0 1 1 - ε dQ(C3) z t ) tp + 1 + u ε dC03

(3A)

Below that straight line, the value of C is either C30 or zero. The region where C ≡ 0 is below the shock path which is determined by eq 22. It starts at the coordinate origin and proceeds as a straight line in the region of the parallel characteristics (where C ) C30) until it intersects the upper of them determined by eq 3A. After that, the shock path proceeds from that intersection point as a solution of eq 22. The solution was obtained numerically with the help of the adaptive Runge-Kutta subroutine odeint.13 To calculate the right-hand side of eq 25, one has to calculate by the above method C(z, t) for the values of z and t prescribed by odeint and use eq 34. As mentioned in the section 2.2, the shock pass itself is the 0-characteristic of eq 24. Other characteristics indicated in section 2.2 as t0-characteristics can be calculated by the same method as the 0-characteristic, the only difference being that they proceed not from the (0, 0) point but from the (0, t0) points. The variation of x1, x3, and QH along a t0-characteristic for the pulse border condition of eq 35 can be calculated by the method described in section 2.2 for the frontal border condition of eq

30. The only difference is that in the former case the value of C is not constant along a t0-characteristic as in the latter one and has to be calculated by the method described above. Correspondingly, the values of λ in eq 29a (for k1 ) 0; see below) and Q in eq 25 should be calculated using eq 34. References and Notes (1) Lin, B.; Song, F.; Guiochon, G. J. Chromatogr., A 2003, 1003, 91. (2) Langer, S. H.; Yurchak, J. Y.; Patton, J. E. Ind. Eng. Chem. 1969, 61, 10. (3) Trapp, O. J. Chromatogr., A 2008, 1184, 160. (4) Guiochon, G.; Lin, B. Modeling for PreparatiVe Chromatography; Academic Press (An imprint of Elsevier Science): Amsterdam, 2003. (5) Guiochon, G.; Felinger, A.; Shirazi, D. G.; Katti, A. M. Fundamentals of PreparatiVe and Nonlinear Chromatography, 2nd ed.; Elsevier, Academic Press: Amsterdam, 2006. (6) Tapp, O. Anal. Chem. 2006, 78, 189. (7) Bakaev, V. A. J. Phys. Chem. C 2007, 111, 7463. (8) Bakaev, V. A.; Bakaeva, T. I.; Pantano, C. G. J. Phys. Chem. C 2007, 111, 7473. (9) Brendle, E.; Papirer, E. In Powders and Fibers. Interfacial Science and Applications; Nardin, M., Papirer, E., Eds.; CRC Press: Boca RatonLondon-New York, 2007. (10) Rhee, H.-K.; Aris, R.; Amundson, N. R. First-Order Partial Differential Equations: Vol. 1. Theory and Applications of Single Equations; Prentice Hall: Englewood Cliffs, NJ, 1986. (11) Rhee, H.-K.; Aris, R.; Amundson, N. R. First-Order Partial Differential Equations: Vol. 2. Theory and Application of Hyperbolic Systems of Quasilinear Equations; Prentice Hall: Englewood Cliffs, NJ; 1989. (12) Zhuravlev, L. T. Colloids Surf., A 1993, 74, 71. (13) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recepies in FORTRAN. The Art of Scientific Computing; 2nd ed.; Camridge University Press: Cambridge, 1992. (14) Bakaev, V. A.; Pantano, C. G. J. Phys Chem C 2009, DOI: 10.1021/jp810984r.

JP810982W