Anal. Chem. 1996, 68, 2147-2154
Numerical Solution of the Complete Mass Balance Equation in Chromatography Gitti L. Frey and Eli Grushka*
Department of Inorganic and Analytical Chemistry, The Hebrew University, Jerusalem, Israel
A new approach to simulate the movement of bands through a chromatographic column is presented. Similar to the Craig distribution model, the mass balance equation is divided into two equations describing two successive processes. The first equation includes two effects: solute diffusion in the mobile phase and migration of the solute band with the mobile phase. The second equation deals with the distribution of the solute between phases, i.e., the adsorption isotherm. The partial differential equations are integrated numerically over time and space using two methods. The first approach is a finite difference method. In the second approach, the propagation operator is expanded in a Chebyshev series, where large time steps can be used. The rate of adsorption and desorption is determined by the size of the time increment. By varying the size of the time step, it is possible to study kinetic effects. The influences of sample size, injection width, rate of mass transfer, and mobile phase velocity on the elution profile were studied. Simulations using the modified Craig approach with either of the two numerical procedures showed that the solutes behaved in the chromatographically expected manner. Moreover, with linear adsorption isotherms, direct relationships between HETP, as well as retention times, and experimental parameters could be established. The importance of computer simulation of the chromatographic process is that it deepens our theoretical understanding of that process and assists us in selecting optimal experimental conditions for a given separation. Frequently, the computer simulation is based on the numerical solution of the chromatographic mass balance equation. The derivation of that mass balance equation and the assumptions made to obtain it have been reviewed extensively (e.g., ref 1 and references therein) and need not be repeated here. The mass balance equation is written as follows:
∂Cs ∂Cm ∂2Cm ∂Cm +F +u )D 2 ∂t ∂t ∂z ∂z
(1)
where Cm is the solute concentration in the mobile phase, Cs is the solute concentration in the stationary phase, F is the phase ratio, z is the distance, t is the time, u is the mobile phase velocity, and D is the axial diffusion coefficient of the solute in the mobile phase. (1) Guiochon, G.; Golshan-Shirazi, S.; Katti, A. M. Fundamentals of Preparative and Nonlinear Chromatography; Academic Press: Boston, MA, 1994. S0003-2700(96)00220-X CCC: $12.00
© 1996 American Chemical Society
In the case of linear chromatography, and under the assumption of fast mass transfer, an analytical solution to the equation is obtained. Under these conditions, the band profile has a Gaussian shape. However, if the mass transfer is not fast or if the isotherm is more complex, an analytical solution is, most frequently, not possible, and numerical analysis is needed. The mass balance equation is a second-order differential equation, and several assumptions are usually made to facilitate the numerical solution. The most common numerical method used is finite difference. By assuming that the column efficiency is infinite (i.e., neglecting longitudinal diffusion), and using a finite difference procedure, the ideal model is obtained.2,3 Since actual columns have a finite efficiency, and since it is not possible to calculate numerically exact solutions of the ideal model,2,4 a semiideal model has been proposed by Guiochon and coworkers.5,6 The semiideal model reintroduces the concept of column efficiency in the calculation by replacing the resistance to mass transfer and the axial diffusion with an apparent diffusion. However, this apparent diffusion process is not included directly in the actual numerical solution, the assumption being made that the numerical error involved in the solution compensates for its omission.5 Using the finite difference procedure, the semiideal model works only when the space increment is equal to the column HETP. However, this limitation is of mathematical origin and does not have any physical justification.7 Early theories of chromatography assumed immediate equilibration of solute between the mobile and stationary phases. This assumption requires an extremely slow flow rate through the column. When the flow is slow, diffusion cannot be neglected. When the flow is fast, diffusion is not significant; however, some sort of a rate expression must be used to take into account the finite time of solute adsorption and desorption. Lapidus and Amundson emphasized the importance of longitudinal diffusion on peak shape.8 Similarly, we focused our efforts on solving numerically the complete mass balance equation. In the model suggested here, longitudinal solute diffusion is not neglected. Moreover, no rate expression is required to describe the chromatographic process. The finite exchange time between the two phases is accomplished by combining the mass balance equation with the Craig distribution model.9 (2) Rouchon, P.; Schonauer, M.; Valentin, P.; Guiochon, G. Sep. Sci. Technol. 1987, 22, 1793. (3) Czok, M.; Guiochon, G. Anal. Chem. 1990, 62, 189. (4) Rouchon, P.; Schoenauer, M.; Valentin, P.; Vidal-Madjar, C.; Guiochon, G. J. Phys. Chem. 1985, 89, 2078. (5) Guiochon, G.; Golshan-Shirazi, S.; Jaulmes, A. Anal. Chem. 1988, 60, 1856. (6) Katti, A. M.; Guiochon, G. A. In Advances in Chromatography; Giddings, J. C., Grushka, E., Brown, P. R., Eds.; Dekker: New York, 1992; Vol. 31, p 1. (7) Lin, B.; Guiochon, G. Sep. Sci. Technol. 1989, 24, 31. (8) Lapidus, L.; Amundson, N. R. J. Phys. Chem. 1952, 56, 1952.
Analytical Chemistry, Vol. 68, No. 13, July 1, 1996 2147
Two different numerical algorithms were used to compute this model. One numerical method is based on backward (in the time space) and central (in the distance space) finite differences. The second algorithm used Chebyshev expansion (discussed in detail in the next section) in the time space.10 The accuracy and stability of the two algorithms were studied. Also examined were the influences of sample size, width of injection, and mobile phase velocity. Qualitative results obtained, i.e., peak shapes, confirm that the suggested model is appropriate. In addition, quantitative results relating peak shape and width to experimental parameters are empirically obtained and discussed. The aim of the present work is to examine the feasibility of using the Chebyshev approach. To conserve computation time in this preliminary study, we have examined cases with gaseous diffusion coefficients and short columns. Future work will investigate cases more typical of liquid chromatography and capillary electrophoresis. THEORETICAL SECTION In the present treatment, a modified Craig distribution model is used. In this model, the continuous column is replaced with a segmented column. Initially, each segment contains only the stationary phase. Mobile phase with solute is introduced into the first segment. After equilibrium has been achieved, the mobile phase with the solute in it is trensfered, with velocity u, from the first segment into the next, and the first segment is filled with a fresh portion of mobile phase without the solute. A new equilibrium is achieved, now in two segments. After the equilibrium, the mobile phase from one segment is moved ahead to the next one. The whole process is repeated until the solute is eluted from the column. In this model, the chromatographic process is described as the sum of two steps. The first step is the moving of the mobile phase from one segment to the next, equivalent to the flow of mobile phase through the column. In the second step, the solute is distributed between the two phases and a new equilibrium is reached, equivalent to adsorption and desorption of the solute in the column. In the original Craig approach, diffusion was not taken into account. In the present treatment, longitudinal diffusion of the solute in the mobile phase is included in the segment-to-segment transfer stage. Also, in the present model, the movement from segment to segment is done with the velocity of the mobile phase. In other words, in the modified Craig model described here, the column segments are not physically separated from one another as in the original model. A somewhat similar model of dividing the column into segments and separating the solute movements into two parts, one to take care of diffusion and forward flow and the other to account for adsorption, is discussed by Jaulmes and Vidal-Madjar.11 The modified Craig model can be described by the following set of equations:
∂2Cm ∂Cm ∂Cm )D 2 -u ∂t ∂z ∂z
(2)
Cs ) f(Cm)
(3)
rium between the phases. Since in this model equilibrium is reached at each time increment, ∆t, of the numerical integration, the size of that increment defines the rate of adsorption. A small time increment mimics fast equilibrium, and a large time increment illustrates slow kinetics. Alternatively, because solute distribution is assigned to an independent equation, it is possible to simulate the distribution between the phases once every few time steps, stimulating a chromatographic system with very slow adsorption-desorption kinetics. In the numerical analysis of the modified Craig model described here, the column is divided into length segments, ∆z, and the separation time into ∆t segments. The result is a concentration grid, C(z,t), detailing the solute concentration in the mobile and stationary phases in the zth slice of the column after the tth time interval. C(z,t) is calculated in two stages, in accordance with the Craig steps. (a) Changes in analyte concentration in mobile phase due to flow and diffusion are calculated for every ∆z slice. This is accomplished by computing eq 2 at grid points C(1...n,t). (b) Analyte concentration in mobile phase Cm(z,t) and in the stationary phase Cs(z,t) are summed and placed into the isotherm function (eq 3), after which a new equilibrium is reached in all grid points C(1...n,t). The time scale is then advanced to t + 1, and the two steps are repeated to calculate C(1...n,t+1). Equation 3 is simple to compute and does not require much computer time. Solving step 1, i.e., eq 2, while also not difficult, consumes most of the computation time. As will be discussed shortly, a finite difference method can be employed for the integration. However, two major difficulties arise with this approach. First, the space increment, ∆z, must equal the quantity 2D/u (i.e., HETP); otherwise the results are not trustworthy. It should be stressed that this limitation in the size of ∆z is of mathematical origin and not a limitation of the physical model itself. Second, there is a ∆t-related computational error which affects the accuracy of the peak shape and retention time. For best results, ∆t must be as small as possible, a constraint which inceases the computation time. As will be shown later, the size of ∆t dictates the kinetics of mass transfer between the phases. Therefore, the size of ∆t can influence calculated H values due to numerical errors and due to adsorption-desorption kinetics effects. To avoid these constraints, a different algorithm is suggested. This numerical procedure was used previously to simulate the dynamics of two-dimensional barrier crossing12 and to calculate eigenfunctions and eigenvalues of the Schrodinger equation.13 It is based on propagating the time scale via a Chebyshev series. In this case, there is no numerical error which depends on ∆t, and relatively large time increments can be applied. Presently, we will compare the Chebyshev algorithm with the finite difference approach. I. Finite Difference. Since the forward finite difference approximation is conditionally stable, we use instead a combination of central and backward difference approximations to simulate eq 2. All finite difference approximations lead to a truncation error, defined as the difference between the true derivative and the finite difference approximation of it. This truncation error can be estimated by a Taylor series expansion.14 The central difference
The second equation is the isotherm which dictates the equilib(9) Craig, L. C. J. Biol. Chem. 1943, 150, 33. (10) Pines, E.; Huppert, D.; Agmon, N. J. Chem. Phys. 1988, 88, 5620. (11) Jaulmes, A.; Vidal-Madjar, C. Anal. Chem. 1991, 63, 1165.
2148
Analytical Chemistry, Vol. 68, No. 13, July 1, 1996
(12) Agmon, N.; Kosloff, R. J. Phys. Chem. 1987, 91, 1988. (13) Kosloff, R.; Tal-Ezer, H. Chem. Phys. Lett. 1986, 127, 223. (14) Celia, M. A.; Gray, G. G. Numerical Methods for Differential Equations; Prentice Hall: Englewood Cliffs, NJ, 1992.
approximation has the smallest error. However, this approximation needs information from more nodes (grid points) than do the forward or backward approximations. Thus, the central difference approximation requires longer computation time. By using the backward difference operator for the LHS of eq 2, and using the central difference operator for the RHS, we can rewrite eq 2 as
[
Cm(z,t-1) ) -
]
[
] ]
D∆t u∆t 2D∆t Cm(z+1,t) + 1 + Cm(z,t) + 2 2∆t (∆z) (∆z)2 D∆t u∆t Cm(z-1,t) (4) (∆z)2 2∆z
[
Equation 4 is used to estimate the solute concentration in the mobile phase at every z position along the column at time t. Next, the isotherm expression (eq 3) is employed to calculate the new equilibrium concentrations of the solute in the mobile and stationary phases at the same time interval t. At this point, the time increment is advanced to t + 1, and the whole process is repeated. II. Chebyshev Expansion. To solve eq 2 using the Chebyshev expansion algorithm, it is necessary to rewrite that equation in an operator form:
∂Cm(z,t) ∂t
then eq 9 and 10 can be solved using a forward difference approximation to give
A1/2 )
e-ν(z+1/2) ν(z+1) [e Cm(z+1,t) - eν(z)Cm(z,t)] ∆z
(13)
e-ν(z-1/2) ν(z) [e Cm(z,t) - eν(z-1)Cm(z-1,t)] ∆z
(14)
and
A-1/2 )
Substituting eqs 7, 8, 13, and 14 into eq 5 gives10
∂Cm(z,t) ∂t
) Ω(z/z - 1)Cm(z-1,t) + Ω(z/z + 1)Cm(z+1,t) [Ω(z - 1/z) + Ω(z + 1/z)]Cm(z,t) (15)
where Ω(z/z - 1) is the transition probability per unit time from z - 1 to z. The symbol “/” in the argument of Ω does not mean division but rather the separation between two neighboring distance intervals. That transition probability is equal to
Ω(z/z - 1) ) ) OCm(z,t)
D eu∆z/2D (∆z)2
(16)
(5) and Ω(z/z + 1) is the transition probability per unit time from z + 1 to z and equals
where
∂ ∂ O ) D euz/D e-uz/D ∂z ∂z It is convenient to define a local reduced velocity, v(z), as
ν(z) ) -uz/D
Ω(z/z + 1) )
(6)
(7)
To obtain the finite difference portion (for the spatial dimension) of the Chebyshev method, the operator in eq 6 is written as
D (A - A-1/2) ∆z 1/2
(8)
(17)
According, the transition probabilities from z to the two nearest neighbors z - 1 and z + 1 are
Ω(z/z - 1) + Ω(z/z + 1) )
O)
D e-u∆z/2D (∆z)2
D D eu∆z/2D + e-u∆z/2D (∆z)2 (∆z)2 (18)
The propagation in time is based on the formal solution of eq 5:
Cm(z,t+1) ) e∆tOCm(z,t)
(19)
where
∂ A1/2 ) e-ν(z) eν(z)Cm(z+1/2,t) ∂z
(9)
and
Following Kosloff and Tal-Ezer,13 we expand the exponent in eq 19 in a Chebyshev series. This series converges uniformly in the interval [-1,1], which allows propagation in relatively large time steps. In the Chebyshev expansion, nb
-ν(z) ∂
A-1/2 ) e
∂z
ν(z)
e
Cm(z-1/2,t)
e
(10)
∆tO
)e
-R∆t/2
∑a T
n)0
If we define
(
2O
n n
R
)
+I
(20)
where 1
1
ν(z + /2) = /2[ν(z) + ν(z + 1)]
(11)
and
a0 ) I0(R∆t/2)
(21)
an ) In(R∆t/2)
(22)
Analytical Chemistry, Vol. 68, No. 13, July 1, 1996
2149
and for n > 0,
ν(z - 1/2) = 1/2[ν(z) + ν(z - 1)]
(12)
In is the nth modified Bessel function of the first kind. R is the upper bound [dD(π/∆z)2] for the range of eigenvalues of the operator O, and Tn is the nth Chebyshev polynomial,
Tn(cos θ) ) cos(nθ)
(23)
calculated by its recursion relation,
Tn(X) ) 2XTn-1(X) - Tn-2(X)
(24)
were T0(X) ) I is the identity operator, T1(X) ) X, and X ) 2O/R + 1 (for further details regarding this algorithm, see refs 10, 12, 13, 15, and 16). III. Isotherm Used. Systems including linear and Langmuir isotherms have been investigated. A linear isotherm is given by
Cs(z,t) ) KCm(z,t)
C(0,t) ) 0
(30)
(25)
where K is the partition coefficient. The numerical solution to the mass balance equation for linear chromatography can be pursued in two ways. In the first way, the isotherm expression is introduced directly into the mass balance equation. In this case, eq 1 can be written as (when F ) 1)
∂2Cm ∂Cm ∂Cm ) D′ 2 - u′ ∂t ∂z ∂z
Since the mass balance equation requires boundary and initial conditions for a physically meaningful solution, we need to define these constraints irrespective of the method of solution used. The initial and boundary conditions which we chose were identical for both methods of solution. Initial Conditions. The solute concentrations at all positions in the column at time t ) 0 are the initial conditions. Different profiles can be used: Gaussian, Dirac pulse function, step function, plug, etc. We used a plug function as the initial condition and studied the influence of plug height (sample size) and width (injection width) on the elution profile. Boundary Conditions. It is necessary to set boundary conditions for the two ends of the column. At the column inlet, i.e., z ) 0, we need to prevent the solute from leaving the column. For that purpose, a reflection boundary condition is set at z ) 0:
(26)
where
D′ )
D 1+K
(27)
u′ )
u 1+k
(28)
and
Equation 26 is identical to eq 2, but D and u are replaced by D′ and u′. Solving eq 26 eliminates the need to use two equations to simulate linear chromatography. Both numerical approaches, the finite difference and the Chebyshev, can be used in the numerical solution of eq 26. The second way to a numerical solution to the mass balances equation is the modified Craig procedure. Here, the isotherm is not included in the propagation operator. Rather, the isotherm is employed to calculate equilibrium concentrations at every slice of the column subsequent to the solution of eq 2. The second isotherm studied was the Langmuir isotherm:
At the column end, i.e., z ) n, the boundary should be an absorbing one, that is, a solute sink. To ensure that the sink does not interfere with the detection, the actual column length was made much larger than the detection length. In the subsequent discussion, the term “length” will mean detection length. COMPUTATIONAL DETAILS The computations have been performed on a 486 PC running at 33 MHz. All programs were written in Visual Basic. The programs were not optimized. Therefore, computation time for the finite difference algorithm varied from a few seconds to 1 h, while the Chebyshev expansion algorithm, being more complex, could take up to several hours for the solution. We have studied the effects of sample size, injection width, and mobile phase velocity on the retention time and plate height, both calculated from the moments of the simulated peaks.17 Since all the numerical calculations were done on a PC, we limited the simulations to chemical systems with large diffusion coefficients, i.e., gas chromatography and very short columns. For liquid systems, the Chebyshev algorithm required long computation times. Since the aim of the present work was to check and compare two algorithms as well as two approaches to the simulation, we felt that no loss in generality is suffered from the use of a solute with large D value.
In this case, the complete mass balance equation cannot be rearranged to yield an equation similar to eq 26, and we need to resort to the modified Craig approach.
RESULTS AND DISCUSSION Size of Time and Space Increments. The sizes of the time and space increments determine the accuracy and stability of the numerical method used. Since two different algorithms were employed, it is necessary to study the effect of discretization of the distance and time coordinates on the results obtained from each algorithm. The influence of increment size was studied by plotting HETP and retention time versus ∆t and ∆z in four cases: (a) unretained solute, i.e., no distribution between two phases (only eq 2 is solved); (b) retained solute with a linear adsorption isotherm, where the isotherm is inserted in the mass balance equation (eq 26 is solved); (c) retained solute with a linear adsorption isotherm using the modified Craig method (eqs 2 and 25 are solved); and (d) retained solute with a Langmuir adsorption isotherm using the modified Craig method (eqs 2 and 29 are solved).
(15) Nadler, W.; Schulet, K. J. Chem. Phys. 1986, 84, 4015. (16) Agmon, N. J. Chem. Phys. 1984, 80, 5049.
(17) Grushka, E.; Myers, M. N.; Schettler, P. D.; Giddings, J. C. Anal. Chem. 1969, 41, 899.
Cs(z,t) )
2150
aCm(z,t) 1 + bCm(z,t)
Analytical Chemistry, Vol. 68, No. 13, July 1, 1996
(29)
Figure 1. HETP as a function of ∆t for finite difference algorithm. The relevant parameters are column length, 10 cm; D, 0.1 cm2/s; u, 1 cm/s; ∆z, 0.2 cm. Line A is for no adsorption; line B is for linear isotherm with K ) 2 using eq 26; line C is the same as line B but using the modified Craig approach; and line D is for the Langmuir isotherm with the constants a ) 2 and b ) 1. The injection width was equal to ∆z; solute amount was arbitrarily set equal to 1.
As mentioned before, a limitation of the finite difference method is that, for an accurate and stable solution of eq 2, ∆z must equal 2D/u. This requirement, which can be proven rigorously, was also observed by Guiochon and Lin.7 The accuracy of the simulation depends also on the size of ∆t. The smaller the ∆t, the better the approximation, the best solution being in the limit of ∆t ) 0. To calculate the exact value of HETP of our simulations, we solved the relevant equations at varying ∆t and extrapolated to zero time increment. Figure 1 shows plots of HETP as a function of ∆t when the finite difference algorithm is used for the four cases discussed above. The figure shows that HETP is a linear function of ∆t in all cases. As expected from eqs 2 and 26, the HETP for the cases of no adsorption and of ideal linear adsorption isotherm (i.e., infinitely fast equilibration of the solute between the mobile and stationary phases) should be identical. Indeed, Figure 1 shows that, in the limit of ∆t ) 0, the lines corresponding to no adsorption and to ideal linear adsorption (lines A, B, and C in Figure 1) converge to the theoretical value of 0.2 cm. It is immaterial whether the adsorption isotherm is included directly in the mass balance equation or whether it is coupled to it. However, the slopes of these three lines are not equivalent, which means that the three modes of calculations have different errors associated with them. The slope of line A, for the case of an unretained solute (K ) 0), represents the numerical error caused by the finite difference approximation. The slope of line B, for the case of linear isotherm included in the mass balance equation, i.e., eq 26, is also due to numerical error. However, since the effective diffusion coefficient and effective velocity in eq 26 are smaller than the true values in eq 2, the associated errors are also smaller, resulting in a smaller slope. In fact, the slope ratio of line A to B should be 1 + K, which for the data in Figure 1 is 3. The slope of line A is 0.96, and that for line B is 0.32, and the ratio between them is, indeed, 3. Similar results were obtained for other simulations with different K values. The slope of line C, belonging to retained solute with a linear adsorption isotherm solved using the modified Craig approach, is larger than the slopes in the two previous cases. The reason for the higher slope of line C is as follows:
Figure 2. Effect of ∆t on the peak shape in linear chromatography. Finite difference algorthm using the modified Craig approach. Same parameters as in Figure 1.
Early models discussing linear chromatography predict a Gaussian shape elution band, providing the mass transfer coefficient is large. However, if the mass transfer coefficient is small, yielding a column efficiency of 0), HETP shows the expected minimum. As can be seen in Figure 7, in the case of slow kinetics, the plate height increases with the partition coefficient. Also, the plate height for Langmuir isotherm with a given a value is higher than the equivalent linear isotherm with K ) a. CONCLUSIONS The modified Craig approach described here allows us to solve numerically the mass balance equation in chromatography without neglecting the molecular diffusion term. It also allows us to simulate kinetic effects by varying the size of the time discretization. Whether the linear or Langmuir adsorption isotherm was
2154 Analytical Chemistry, Vol. 68, No. 13, July 1, 1996
used, the results of the calculations behaved as expected for chromatographic systems. The numerical solution was carried out using two different approaches. One approach is based on a combination of backward and central finite difference. The second approach uses a Chebyshev expansion to propagate in the time dimension. The first procedure suffers from the limitation that the distance discretization is limited to one value only, i.e., 2D/u. Using any other value of ∆z results in the loss of material during the simulation. In the Chebyshev procedure, this restriction is relaxed, and the ∆z is limited only by the desired accuracy. Another advantage of the Chebyshev approach is that the size of ∆t does not limit the accuracy of the results. Thus, large time increments can be used to save computing time. The two numerical approaches to the solution of the chromatographic mass balance equation can also be summarized as follows. The Chebyshev approach is more complex and more difficult to program than the finite difference method. Moreover, the computation time is longer with the Chebyshev approach. However, the great advantage of that approach is that the results of the simulation are independent of the size of the time interval ∆t. The present study was limited to chromatographic systems with a solute having a large diffusion coefficient and a short column. We are now examining cases more closely resembling liquid chromatography. In addition, we are also attempting to solve the complete mass balance equation, using the Chebyshev procedure, which will include directly in it the distribution between the two phases irrespective of the isotherm used and of the mass transfer coefficient. ACKNOWLEDGMENT We thank Prof. R. Kosloff and Prof. N. Agmon for fruitful discussion and help with the Chebyshev algorithm. This research was supported by Grant No. 88-00021 from the United StatesIsrael Binational Science Foundation (BSF), Jerusalem, Israel.
Received for review March 6, 1996. Accepted March 26, 1996.X AC960220O X
Abstract published in Advance ACS Abstracts, May 1, 1996.
