I n d . Eng. Chem. Res. 1988,27, 1474-1480
1474
Numerical Simulation of a Chromatograph Column: Linear Case Kevin J. Ward, Serge C. Kaliaguine, a n d Philippe A. Tanguy* Department of Chemical Engineering, Universit6 Laval, QuQbec,Canada G1K 7P4
Gilles J e a n CANMET Energy Research Laboratories, Energy, Mines and Resources Canada, Ottawa, Ontario, Canada K I A OGl
An algorithm for the determination of the mass-transfer parameters in the Kubin-Kucera adsorption model of a packed chromatograph column is developed in this study. While the adsorption system modeled in this investigation was linear, the methodology is general enough to be extended to nonlinear adsorption. The model mass balance equations are solved through finite differences, and the computed response curve is fitted to an experimental chromatogram by using a nonlinear least-squares routine for the optimization of the mass-transfer parameters. The mass-transfer parameters determined by the devised algorithm were found to be in excellent agreement with those given in the literature for a linear adsorption system of nitrogen on activated carbon. Early studies of gas chromatograph systems involved simple mathematical models such as the Van Deemter equation, based on a single diffusivity (Ma and Mancel, 1972). A more exact model which analyses the various mass-transfer steps in the packed column was proposed by Kubin (1965) and Kucera (1965) and was later developed by Schneider and Smith (1968). The model represents the column as a one-dimensionalpacking of spherical particles and takes into consideration the mass-transfer processes in the bulk gas phase, at the gas-particle interface, and in the pore volume of the particles. However, the Kubin-Kucera model does not take into account the bidisperse pore structure of most commercial adsorbents, which consist of small microporous particles formed into macroporous pellets or adsorbents. In such adsorbents three distinct resistances to mass transfer exist due to the mass transfer at the bulk gas-particle interface, the diffusion in the macropores, and the diffusion in the micropores. A suitable mathematical model for these systems was developed by Sarma and Haynes (1973) and is now widely used. For a particular chromatograph column model, the pertinent mass-transfer parameters can be estimated by various techniques. Classical parameter estimation techniques, in general, consist of fitting the experimental chromatogram to a mathematically predicted response curve in the Laplace, frequency, or time domain. The fitting techniques in the Laplace and frequency domains are commonly referred to as moment methods and Fourier analyses, respectively. Each of these parameter estimation methods requires derivation of the Laplace transfer function from the linear or linearized differential equations of the mass-transfer model. For nonlinear adsorption systems, in which transfer functions cannot be used, the nonlinear model must then be solved by numerical techniques to obtain the predicted response curve. Many adsorption systems, such as the adsorption of nitrogen compounds on zeolites, exhibit nonlinear behavior. The nonlinearity of a system can be established experimentally by observing the effect of varying the input signal concentration on the response signal. If the shape and retention time of the response signal are found to be sensitive to the size of the input pulse, the system is nonlinear (Ruthven, 1984). The system nonlinearity may be attributed to non-Fickian diffusion within the adsorbent
* To whom correspondence should be addressed.
pores, nonlinear intrinsic adsorption kinetics, and nonisothermal sorption. Very few studies have been reported in the literature for parameter estimation for nonlinear adsorption systems. Since nonlinear systems are frequently encountered (Ruthven, 1984), the objective of the present study is to develop a numerical methodology for physical parameter estimation in such systems. Theoretical Model For the sake of simplicity, the Kubin-Kucera model was used in the present study. The model includes mass balances for the bulk gas phase in the column and the pore volume gas phase within the adsorbent particles. Models for bidispersed adsorbents (Sarma and Haynes, 1973), on the other hand, would include another mass balance for the pore volume gas phase within the microparticles contained in the larger adsorbent particles. The inherent assumptions in the Kubin-Kucera chromatograph column model are the following: the adsorbate concentrations are in the Henry’s law region of the adsorption isotherm, the column operates isothermally, the adsorbent particles are spherical and of average radius R, the concentration profiles are always symmetrical within the adsorbent particles, the diffusion in the adsorbent particles is Fickian, the flow in the column is axially dispersed, and two mass-transfer resistances exist due to mass transfer through the film at the bulk gas particle interface and diffusion within the adsorbent pores. The mass balance in the bulk gas phase of the chromatograph column is given as
where D , is the axial dispersion coefficient (m2/s), C is the adsorbate bulk gas concentration based on the empty column volume element (mol/m3),z is the axial distance from the column inlet (m), V is the interstitial bulk gas velocity (m/s), t is the time variables (s), q-,is the packed bed void fraction, R is the average adsorbent particle radius (m),and Qf is the molar mass-transfer flux through the film at the bulk gas-particle interface (mol/(m2s)) as defined by (2) Qf = kf(C - (CJR) where kf is the fluid-to-particle mass-transfer coefficient (m/s), C, is the adsorbate pore volume concentration, based on the empty column volume element (mol/m3),and
0888-~885/88/2621-1474$01.50/0 0 1988 American Chemical Society
Ind. Eng. Chem. Res., Vol. 27, No. 8, 1988 1475
(Cp)Ris C, evaluated at the external particle surface (mol/ m3). Equation 1 is subject to the following initial and boundary conditions: C(z,O) = 0
(3)
C(0,t) = f ( t )
(4)
ac ( m , t ) az
(5)
=0
't
0
I
where f ( t ) is the input concentration disturbance (mol/(m2
where tp is the intraparticle void fraction, De is the intraparticle effective diffusivity (m2/s),pp is the adsorbent particle density (kg/m3),r is the radial distance from the center of the adsorbent particle (m), and Qad is the intrinsic adsorption rate (mol/(kg s)) given by (7)
where k,' is the adsorption rate constant (m3/(kg s)), KA' is the adsorption equilibrium constant (m3/mol), and Cad is the adsorbed equilibrium adsorbate concentration on the adsorbent surface (mol/kg). Equation 6 is subject to initial condition 8 and boundary conditions 9 and 10:
Cp(z,r,O)= 0
(8)
I
I
Xi
Xi-1
8)).
The mass balance for diffusion within an adsorbent particle is given as
0 0
Xitl
Figure 1. Typical finite difference mesh.
correspond to the computation of the values of the dependent variable (C or Cp). Since the initial and boundary conditions associated with the problem give information at the circled mesh points, it is clear from the figure that no explicit procedures could be used to solve the problem. The advantage of an implicit method over an explicit method is that implicit methods are unconditionallystable and thus permit the use of larger time steps (Vemuri and Karplus, 1981). In the developed finite difference scheme, the accumulation term for both the bulk gas and particle-gas mass balances is represented by a first-order Euler scheme for the first time step and by a second-order Gear scheme for subsequent time steps (Fortin et al., 1986). The Gear scheme enables the use of larger time steps, but it cannot be used until the second time step. This accuracy improvement results from the centering of the Gear finite difference at the current time step ( j ) while that of Euler is centered between the jth and (i-11th time steps (Vemuri and Karplus, 1981). The implicit finite difference representation of either mass balance yields a tridiagonal matrix system to be solved algebraically. The problem at hand is then to resolve an equation of the form [AIX = W
t 10) The physical parameters, contained in this model, which characterize the shape of the response curve include D,, De, kf,KAI, and k;. Parameter Estimation Methodology The parameter estimation methodology developed in the present study involves two steps: a finite difference solution of the model equations to obtain the predicted response curve, and the computation of improved parameter values by a nonlinear least-squares fitting of the predicted response curve to an experimental chromatogram. Finite Difference Scheme. Finite difference resolution of differential equations involves discretization of the mass balance equations within given domains. The resulting algebraic equations, which may be explicit or implicit, are then solved to obtain values of the independent variables at discrete nodes in each domain. An implicit finite difference scheme was developed for resolving the model equations. The standard technique used consists of specifying the time derivatives (accumulation terms) at the j t h and previous time steps and specifying the spatial derivatives also at the j t h time step. This is shown in Figure 1 where the abscissa ( x ) is the independent variable which is representative of z or r, the ordinate ( t )is the time, and the length of the spatial domain is Xo. The interior points denoted by X s then
(11)
where
W = f(Ci,,-.l)
Euler scheme
W = F(Ci,j-l, Ci,j+J
Gear scheme
(12)
(13)
Matrix A is an LU factorized form (Burden and Faires, 1985) of the matrix of dimensionless parameters derived from the finite difference equations. This factorized matrix is constant for fixed physical and resolution parameter values. Two important matrix solution techniques for such tridiagonal systems are the Thomas algorithm (Burden and Faires, 1985) and the Successive Over Relaxation method (SOR) (Vemuri and Karplus, 1981). However, the Thomas algorithm is generally more efficient, unless the number of grid points is very large, and thus was selected in this study. The computation of the response curve consists of resolving the coupled finite difference equations for the bulk gas and the particle pore volume. Curve Fitting Algorithm. A nonlinear least-squares algorithm was programmed for the fitting of the calculated response curve to the experimental chromatogram. The least-squares fitting requires the minimization of the objective function ( Y ( X ) )
Y ( X ) = ?&(XI2 i=l
X E R"
(14)
where Xis the vector of optimization parameters, n is the number of optimization parameters, m is the number of
1476 Ind. Eng. Chem. Res., Vol. 27, No. 8, 1988
curve fitting points, and 4i(X) is the residual function defined by 4i(x) = C,,,,(X,t) - Ce,*t(t) (15) where Cpred(X,t) is the response curve predicted by a mathematical model, and Cexp,(t)is the experimental response chromatogram. If the mathematical model for predicting the response curve is nonlinear in X, the least-squares problem is nonlinear. The classical solution technique known as the Gauss-Newton method uses a linear model, M,(X),about X, given as (Dennis and Schnabel, 1983) M,(x) = 4(XJ - J(X,)(X- XJ (16) where 4(X,) is the vector of residual function evaluated a t X,,X, is the vector of current parameter values, and J(X,) is the Jacobian matrix defined as J(Xc)i,j = %i(Xc)/d(Xc)j (17) In the event that m and n are equal, Mc(X)can be driven to zero, but when m is greater than n, the matrix system is overdetermined and M,(X) can only be minimized. Thus, the optimal vector, X,, is found which corresponds to the parameter values resulting in the best fit of the response curves. The classical solution of eq 16 is
30
40
50
60
70
80
90
Time ,s
Figure 2. Response curve sensitivity to the number of column nodes
(MB).
bed ( L = 0.204 m, +, = 0.38) of monodispersed, spherical m, tp = 0.59). activated carbon particles ( R = 1.0 X Experiments were conducted a t four flow rates having Reynolds numbers of 0.051, 0.11, 0.30, and 0.47, respectively.
Results Wakao and Kageui (1982) estimated the values of the mass-transfer parameters in the linear Kubin-Kucera Xt = X, - (J(Xc)TJ(Xc))-lJ(X,)T~(Xc) (18) model by curve fitting in the time domain using Fourier where X, is the vector of updated parameter values. Direct series and the Laplace inverse of the model transfer solution of eq 18, however, can cause numerical stability function. They assumed the value of the adsorption rate problems such as underflows and overflows and can square constant (k:) and calculated a value of the mass-transfer the conditioning of the problem. To avoid these problems, coefficient from the correlation the Jacobian is decomposed by a QR factorization as N S h = 2 + 1.1Nsc1/3NR,0.6 (23) discussed by Strang (1976). The iterative least-squares solution then consists of using where NSh,Ns,, and NReare the Sherwood, Schmidt, and Reynolds numbers, respectively. They determined the x, = x, + R-l(Xc)QT(X,)4(X,) (19) values of the mass-transfer parameters to be D, = 4.8 X Since the computation of the Jacobian by numerical m2/s, De = 6.3 x m2/s, and p&A' = 5.29. With differences can be very expensive, a Broyden update for these parameter values and the developed finite difference the computation of the Jacobian involving secant apscheme for the Kubin-Kucera model, numerical generation proximations was incorporated. The Jacobian is updated of the response curve was attempted. Initially, the prinaccording to (Dennis and Schnabel, 1983) cipal resolution parameters including the time-step, the number of column nodes, the number of particle nodes, (20) J(Xi) = J(Xi-1)+ Ai and the convergence tolerance were varied to establish the where A, is computed by respective effects on the predicted curve. The value of the time step which ensured numerical accuracy and stability was found to be dependent upon values of several variables, the most important of which include the number of column nodes (MB), the intraparFor the optimization of the values of the mass-transfer ticle effective diffusivity (De),the axial dispersion coeffiparameters in the Kubin-Kucera model, the calculated cient (D,), the particle-to-fluid mass-transfer coefficient response curve was fitted to the experimental chromato&), and the dimensionless equilibrium constant ( p A'). gram such that the error (12 norm) between the two reSmall time steps are required with large values of fand sponse curves defined as D,, and with small values of De and p&*'. It was also found that At is inversely proportional to MB. The response curve was observed to be quite sensitive to the number of column nodes as shown in Figure 2. As was less than a specified tolerance. the values of the number of nodes are increased, the response curve becomes sharper, but for more than 102 Experimental System particle nodes the increase in sharpness is only slight. Thus, a maximum number of nodes, defined as the Wakao and Kageui (1982) determined the governing quotient of the column length (L)and the average particle mass-transfer parameters by the chromatographic techdiameter (2R), was selected for the finite difference repnique for the adsorption of nitrogen on activated carbon. resentation. For the column in question, the maximum Thus, modeling their column provided a good test for the number of nodes was found to be 102. developed software. In coupling the finite difference resolution schemes for Their adsorption chromatographic measurements were the bulk gas and pore volume mass balances, convergence made a t 20 "C and atmospheric pressure by imposing was based upon the pore volume concentration at the nitrogen on a stream of hydrogen carrier gas in a packed
Ind. Eng. Chem. Res., Vol. 27, No. 8, 1988 1477 50
Table I. P a r a m e t e r Values for the Reference Case
0
40
c
t 0
c u
30
8 I
I
I
1
p 20
I
N
F
b
IO
0
0
IO
20
30
40
50
60
70
80
90
100
Time, s
Figure 4. Response curve sensitivity to the value of the fluid-toparticle mass-transfer coefficient.
Time ,s
Figure 3. Response curve sensitivity to the value of the intraparticle effective diffusivity coefficient.
adsorbent particle surface [ (Cp)R]. The convergence criterion is given as [ ( c p ) R l i - [(cp)Rli-l [(cp)Rli-l
< Cc
(24)
for each particle surrounded by a significant bulk gas concentration. The subscript i denotes the convergence iteration number. The convergence tolerance value (€3on this criterion was found to be dependent upon the number of finite difference nodes in the adsorbent particles (MP). With fewer numbers of particles nodes, smaller values of tCcould be used to improve the accuracy of the computed particle concentration profiles. However, this results in a large computation time. An MP value of 20 and an t, value of 0.05 were found to be the best combination to ensure accuracy with a minimal amount of computation time. With these values of the resolution parameters, the parameter values given by Wakao and Kageui, (1982) and a masstransfer coefficient value of 0.10 m/s (computed by eq 24), a time step ( A t ) ,determined by software testing, of 0.037 s was required. Based on the optimal values of the finite difference resolution parameters (MB = 102, MP = 20, 6, = 0.05, and At = 0.037 s) and the parameter values given for the chromatograph column, a predicted response curve was calculated. As will be seen later on when we discuss Figure 11, this curve labeled “before optimization” was compared with the normalized experimental chromatogram. The experimental chromatogram is normalized using CoeHpt (expressed in arbitrary chromatogram units-seconds)as the normalization factor with
Therefore, the normalized concentration Cexpt/Coexpt is expressed in inverse seconds. Although there is some deviation between the two curves, the difference can be attributed to errors in the mass-transfer parameter estimates,. The “before optimization” curve will be designated as the reference case, for which the parameter values are summarized in Table I.
0 Time, s
Figure 5. Response curve sensitivity to the value of the dimensionless equilibrium constant.
Response Curve Sensitivity to the Parameter Values. The response curve sensitivity was examined with respect to variations in the most important physical parameters including: the axial dispersion coefficient (Dm), the intraparticle effective diffusivity (De),the dimensionless equilibrium constant (p A’) and the fluid-to-particle mass-transfer coefficient (,i$ Figure 3 illustrates the response curve sensitivity to the value of the effective intraparticle diffusivity where the reference case corresponds to De = 6.3 X lo-’ m2/s. If De is increased by an order of magnitude, the response curve sharpens significantly, while decreasing De by a factor of 2 attenuates the response signal quite importantly. The observed increase in sharpness of the response curve for increasing values of De is in accordance with previous work by Haynes (1975). The response curve sensitivity to the value of the particle-to-fluid mass-transfer coefficient (k,) is shown in Figure 4 where k f is varied over an order of magnitude. As kf is reduced from 0.10 to 0.01 m/s, the response is attenuated only slightly. The effect of varying the dimensionless equilibrium constant (p&A’) on the response curve is illustrated in Figure 5 where p&A’ is varied from 2.65 to 10.58. The reference case is for a p&*‘ equal to 5.29. It was observed that as ppKA’ is increased the pulse retention time increases proportionally and the peak height is attenuated significantly. The linear relationship between the pulse retention time and P & ~ ’ is illustrated in Figure 6. Since the equilibrium constant KA’ is defined as k L / k d , it can be stated that, in
1478 Ind. Eng. Chem. Res., Vol. 27, No. 8, 1988 50 I
I
I
I
I
- +,-Exponential
'
I
I
(A= I O S
Exponential ( A =
)
0 30 5 ' ' i
T
-
'0
2
,8
6
4
IO
12
14
Time, s
PP KQ,
Figure 6. Relation between the pulse retention time and the value of the dimensionless equilibrium constant.
50
r---
Figure 8. Response curve sensitivity to the form of the input signal. 25
'
I ~
DO^ , mZ/s 4 8x
, 608
Y
~-~
~~
12 17
0
10
20
30
40
50
60
70
80
93
100
Time , S
Figure 7. Response curve sensitivity to the value of the axial dispersion coefficient.
the absence of any adsorption (ki = 0), the retention time for the chromatograph column in question would be about 15.5 s as given by the ordinate intercept in Figure 6. Consequently, a pulse retention time can provide an insight into the relative rates of adsorption and desorption in an adsorption column. This phenomena is in agreement with the findings of Galan et al. (1975). The importance of the axial dispersion coefficient (D,) on the shape of the response curve, and therefore the accuracy of the parameter estimates which would be obtained in curve fitting procedures, is manifested in Figure 7 where the reference case corresponds to D, = 4.8 X m2/s. When D,, is decreased by an order of magnitude, the response curve is sharpened quite significantly, while increasing D , by a factor of 2 broadens the response curve fairly importantly. Thus, larger D, values tend to broaden the response curve. Response Curve Sensitivity to the Form of the Input Signal. For a chromatograph column with relatively short retention times, the predicted response curve may be sensitive to the mathematical expression used to describe the input signal. While mathematically the input signal for a one-shot input is often described by a delta (impulse) function, it is more accurately represented by a fitting of the experimental chromatogram or a decaying exponential function. It has been shown by Boniface and Ruthven (1985) that a normalized input signal could be represented by an exponential function such as C,*(t) = Ae-At (26) where A is the injection time constant in inverse seconds.
Ind. Eng. Chem. Res., Vol. 27, No. 8, 1988 1479 I
50
25
Experimental curve
I N
0
After optr"ot!on Before optimization
20
'I !,
E
I5
L
8
8 .$ 0
IC
E
b
2
5
0
1 IO
20
30
40
50 60 Time, s
70
80
90
100
radial position at successive time intervals for the first particle in the column. The figure illustrates how the adsorbate concentration first increases near the external particle surface until a maximum normalized concentration of 0.0215 at 3.42 s is attained, at which point the concentration throughout the particle is uniform. After reaching this maximum concentration, the adsorbate concentration first diminishes near the particle external surface. This results in higher pore volume concentrations near the particle center until the concentration becomes zero throughout the particle after about 15.0 s. It should be noted that the concentration at the normalized particle radius of 1.0, as a function of time corresponds, almost identically to the input pulse data for this first particle. From a figure such as Figure 9, the adsorbed adsorbate concentrations ( c a d ) can be computed as well. Since the Henry's law adsorption isotherm was assumed in the Kubin-Kucera model, the (dimensionless) adsorbed concentrations can be computed from the equilibrium relation =PFA'
(cp*)r*
~
10
20- 30
40
50
60
70
80
90
100
Time, s
Figure 10. Propagation of the concentration pulse throughout the chromatograph column.
(C*ad)r*
n "0
1
(28)
where the * denotes a dimensionless quantity. The propagation of the normalized concentration pulse throughout the chromatograph column is illustrated in Figure 10 for the reference case. The concentration pulses are given a t dimensionless axial distances ( z * ) from the column inlet of 0, 0.1, 0.2, 0.4, 0.6,0.8, and 1.0. The first and latter correspond to the inlet and predicted response curves, respectively. The concentration pulse is observed to be attenuated in an exponential fashion. Also, boundary conditions 4 and 5 are satisfied in Figure 10 since the input pulse is at z* = 0 and the concentration vanishes at the end of the column. Optimization of the Parameter Values. The predicted response curve was fitted to the experimental chromatogram using the parameter values of the reference case in Table I as initial estimates. A three-dimensional search in @A', De,and D , was explored, but as was discovered by Wakao and Kageui (1982))the simultaneous determination of Deand D, was not feasible. Numerically, this nonfeasibility arises because both Deand D , influence the shape of the response curve in a similar manner as illustrated in Figures 3 and 7. An accurate value of the axial dispersion coefficient is relatively easy to obtain by a correlation for a given column. Wakao and Kageui (1982) correlated D , with the molecular diffusion coefficient for the column in question. Using their values of 4.8 X m2/s for D,, a two-dimensional search in p S A ' and Dewas attempted. Figure
Figure 11. Comparison of the experimental chromatogram with the calculated response curves before and after optimization. Table 11. Response Curve Points Selected for Fitting PurDoses time, s normalized concn 41.25 0.0143 46.25 0.0293 50.00 0.0378 53.75 0.0407 57.50 0.0377 Table 111. Comparison of the Literature and Optimized Values of p,K,' and D e Darameter lit." oDtimized 5.29 5.37 P#A' 0.63 X 10" 1.01 x 104 De, m2/s "Values obtained by Wakao and Kageui (1982) by using the method of curve fitting in the time domain for the Laplace inversion of the model transfer function.
11 shows the experimental response curve (Wakao and Kageui, 1982, p 44)) the calculated response curve based on the parameter values in Table I, and the calculated response curve fitted to the experimental curve after optimization of the values of p S A ' and De. The experimental chromatogram curve points selected for fitting purposes are given in Table 11. As illustrated in Figure 11,the calculated response curve after optimization of the pSA'and Devalues fits the entire experimental data very well. The predicted response curve was fitted to the experimental chromatogram with an 1, norm of the error, defined by eq 22, of 0.0007 after just three iterations. The computation time for this determination was 2.0 h on a Control Data CYBER 830 computer. The literature values and the optimized values of p S A ' and Deare compared in Table 111. As demonstrated by this table, the dimensionless equilibrium constant and the intraparticle effective diffusivity values were increased by 1.5% and 60%, respectively, to correct the deviation between the two response curves (Figure ll). This rapid convergence was achieved by determining the Jacobian through numerical differences until the parameter values were in close proximity to the solution. Broyden updating of the Jacobian was found to be useful when the five-point l2 norm was less than about 0.001; otherwise, the number of optimization iterations was increased. During an optimization, as the values of the governing mass-transfer parameters change, the predetermined time step may lead to instability in the finite difference resolution. It was found through software testing that the
1480 Ind. Eng. Chem. Res., Vol. 27, No. 8, 1988
smaller time steps were required in the following cases: higher D,,values, higher kf values, lower De values, and lower ppKA' values. For the finite difference resolution of the model for the chromatograph column in question, it was concluded that, if a parameter changes by a certain factor in the direction which requires a smaller time step, the time step should be reduced by about 1.33 times this factor. However, with At's as small as 6.0 x s, it was found that the finite difference resolution was stable for almost any combination of parameter values. Conclusions The development of the finite difference algorithm for solving the Kubin-Kucera adsorption equations for linear systems and fitting the calculated response curve to an experimental chromatogram through nonlinear least squares was successful. The optimized values of the physical parameters obtained for an adsorption system consisting of nitrogen on an activated carbon adsorbent were found to be in excellent agreement with those in the literature. The performance of the developed finite difference resolution scheme is dependent upon the appropriate time step (At)and the number of column nodes required (MB). The number of column nodes required for finite difference resolution is given by the nearest integer to the quotient of the column length ( L ) and the average adsorbent particle diameter (%). Thus, numerically, shorter columns with larger adsorbent particles would be preferred to reduce the amount of computation time. To estimate the values of the governing mass-transfer parameters by the conceived methodology in the linear Kubin-Kucera chromatograph column model, D, should be determined by an independent technique. The p&A' and De values can then be easily found by a two-dimensional search. The developed search algorithm, involving a nonlinear least-squares algorithm with a QR factorization of the Jacobian, performed very well in fitting the calculated response curve to the experimental chromatogram by adjusting the values of pFA'and De. The algorithm, which is a Gauss-Newton type, enables convergence in as few as three iterations when good initial parameter estimates are available. Broyden updating of the Jacobian, rather than using numerical differences, is useful in such curve fitting procedures when the parameter values are in close proximity to the solution values. Otherwise, the rapidity of the convergence is reduced. The advantages in using the developed numerical methodology for parameter estimation include the following: its general nature is amenable to the incorporation of system nonlinearities, parameter values determined by
other methods can be improved, and quantitative information concerning the particle concentration profiles and the propagation of the concentration profiles and the propagation of the concentration pulse is calculated. The findings of this study suggest that additional work should be done. Since nonlinearities could quite easily be incorporated into the finite difference scheme for the Kubin-Kucera chromatographic model, the governing mass-transfer parameters in adsorption systems consisting, for instance, of nitrogen compounds on zeolites could be estimated. Such an estimation is the ultimate goal of our present work. From the computation reported here, it can be, however, recommended that the data be obtained from chromatographic experiments conducted with short columns and large zeolite particles. This would minimize the computation time required in mathematically modeling such systems. Such a work is left for further investigation. Literature Cited Boniface, H. A.; Ruthven, D. M. "The Use of Higher Moments to Extract Transport Data from Chromatographic Adsorption Experiments". Chem. Eng. Sci. 1985, 40, 1401. Burden, R. L.; Faires, J. D. Numerical Analysis; Prindle, Weber & Schmidt: Boston, 1985. Dennis, J. E., Jr.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; Prentice-Hall: Englewood Cliffs, NJ, 1983. Fortin, A.; Fortin, M.; Tanguy, P. ''Numerical Simulation of Viscous Flows in Hydraulic Turbomachinery by the Finite Element Method. Comp. Meth. Appl. Mech. Eng. 1986,58, 337. Galan, M. A.; Suzuki, M.; Smith, J. M. "The Effect of Adsorption Characteristics on Pulse Retention Times". Znd. Eng. Chem. Fundam. 1975,14, 273. Haynes, H. W., Jr. "The Determination of Effective Diffusivity by Gas Chromatography. Time Domain Solutions". Chem. Eng. Sci. 1975, 30, 955. Kubin, M. Collect. Czech. Chem. Commun. 1965, 30, 1104, 2900. Kucera, E. J . Chromatogr. 1965, 19, 237. Ma, Y. H.; Mancel, C. "Diffusion Studies of COz, NO, NOz and SOz on Molecular Sieve Zeolites by Gas Chromatography". AZChE J . 1972, 18, 1148. Ruthven, D. M. Principles of Adsorption and Adsorption Processes; Wiley: New York, 1984. Sarma, P. N.; Haynes, H. W. "A Model for the Application of Gas Chromatography to Measurements of Diffusion in Bidisperse Structured Catalysts". Adv. Chem. Ser. 1973, 133, 205. Schneider, P.; Smith, J. M. "Adsorption Rate Constants from Chromatography". AZChE J . 1968, 14, 762. Strang, G. Linear Algebra and its Applications; Academic: New York, 1976. Vemuri, V.; Karplus, W. J. Digital Computer Treatment of Partial Differential Equations; Prentice-Hall: Englewood Cliffs, NJ, 1981. Wakao, N.; Kageui, S.Heat and Mass Transfer in Packed Beds; Gordon & Breach Science: New York, 1982.
Received for review May 22, 1987 Revised manuscript received April 4, 1988 Accepted April 15, 1988
