Binary retention time method for rapid determination of competitive Langmuir isotherm parameters

Separation and Purification Technology 67 (2009) 344–354

Contents lists available at ScienceDirect

Separation and Purification Technology journal homepage: www.elsevier.com/locate/seppur

Binary retention time method for rapid determination of competitive Langmuir isotherm parameters Arvind Rajendran ∗ , Wenda Chen Nanyang Technological University, School of Chemical and Biomedical Engineering, 62 Nanyang Drive, Singapore 637459, Singapore

a r t i c l e

i n f o

Article history: Received 10 December 2008 Received in revised form 25 March 2009 Accepted 26 March 2009 Keywords: Chromatography Equilibrium theory Competitive Langmuir adsorption isotherm Parameter estimation

a b s t r a c t A method, based on the equilibrium theory of chromatography, to estimate competitive Langmuir isotherm parameters by measuring the retention time of the two shock fronts from pulse injections of binary mixtures has been developed. This method eliminates the need for detector calibration and pure component injections; and computational time required is significantly reduced allowing for rapid parameter estimations. The method has been tested on published results and it has been found that parameters obtained from this “binary retention time method” are comparable to those estimated from more rigorous techniques such as the inverse method. The effectiveness and limitations of this method in rapid determination of isotherm parameters are discussed. © 2009 Elsevier B.V. All rights reserved.

1. Introduction The determination of adsorption isotherms and their parameters represents one of the fundamental challenges that concern chromatographers. Isotherm parameters are the means of understanding the adsorption thermodynamics of solute molecules on the surface of the stationary phase and are the basis for design, scale-up and optimization of adsorption based separation processes [1–7]. Several experimental methods have been developed to determine isotherm parameters each with varying levels of complexity [1,3,8]. The choice of a method is typically a trade-off between the efforts required to measure these parameters and their accuracy. Though static measurements, e.g., gravimetry, volumetry, etc. are increasingly used in the gas adsorption, the techniques used in chromatography have been predominantly dynamic. In short, dynamic methods involve the flow of a fluid (either containing the solute of interest, or the solvent) through a column packed with the stationary phase (i.e., adsorbent) of interest. A known disturbance is introduced at the column inlet, and from the response of the column measured, the adsorption isotherm parameters are calculated. Dynamic methods are preferred as the efforts needed to design and operate them are rather modest compared to the accuracies they yield. The more frequently used dynamic methods are inverse methods [9], elution by characteristic points [10], frontal analysis [11] and perturbation methods [12]. All the methods, barring the perturbation method, typically require performing detector calibrations,

∗ Corresponding author. Tel.: +65 63168813; fax: +65 67947553. E-mail address: [email protected] (A. Rajendran). 1383-5866/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.seppur.2009.03.046

i.e., to convert the detector signal to concentration units and the assumption of an isotherm model. Though methods have been developed for performing detector calibrations, they still represent an extra-step towards the determination of isotherm parameters [13,14]. Other methods, such as the frontal velocity [15,16] and hroot [17], avoid detector calibration and use the information of the velocity of the shock fronts. These methods are suited for step inputs and estimate the isotherm parameters by comparing the calculated velocity of the shock fronts with the experimentally measured one. It is worth mentioning that in order to calculate the velocity of the shock front, it is essential to assume an isotherm model. However, a shortcoming of these methods is that owing to the nature of the input, i.e., a step, large quantities of solute are required and an isotherm model has to be assumed a priori. Further, extending the method to the case of a pulse input is not straightforward as its solution is more complex compared to the case of a step input. Golshan-Shirazi and Guiochon developed the retention time method which overcame one of the limitations of the h-root and frontal velocity method [18]. This method uses the retention time of the shock front and the Henry constant to estimate Langmuir isotherm parameters of pure solutes from elution profiles of pulse injections. This method eliminates the need for detector calibration thereby simplifying the experimental efforts. Further, since a pulse injection is used, the amount of sample required for the isotherm measurement step is minimized. Although it provides rapid determination of isotherm parameters, it is applicable only to single component injections. An extension to the case of binary injections is not simple as the elution of a mixture of two solutes that compete for the adsorption sites on the stationary phase is quite different from the case of the elution of single solutes [19]. In several practical situations, pure substances are often not avail-

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

Nomenclature c C h H J k K L n N N O P Q R S t

v x z

fluid phase concentration characteristic on the physical plane corresponding to simple wave transition parameter defined as h = N1 /N2 Henry constant objective function parameter defined as k = ( − 1)/(K1 ) equilibrium constant length of the chromatographic column adsorbed phase concentration saturation capacity matrix defined as N = [nij ] state corresponding to (c1 = 0, c2 = 0) state corresponding to (c1 = c10 , c2 = c20 ) intermediate state in the chromatographic cycle intermediate state in the chromatographic cycle characteristic on the physical plane corresponding to a shock transition time interstitial velocity dimensionless axial coordinate axial coordinate

Greek symbols  selectivity, defined as  = H2 /H1  characteristic on the hodograph plane b void fraction of the column  parameter defined as dc1 /dc2 1 , 2 solutions of the fundamental differential equation given by Eq. (A.4)  Eigen value of matrix N ˜ parameter in Eq. (A.20)  column phase ratio,  = (1 − b )/b  parameter defined as  = H2  slope of a characteristic on the physical plane

image of a shock transition on the hodograph plane  dimensionless time parameter defined as = c10 /c20 ω characteristic parameters defined by Eqs. (A.11) and (A.12) Subscripts and superscripts 0 conditions corresponding to the injected concentration calc calculated exp experimental i component index k experiment index S corresponding to a shock transition

able (e.g., in enantioseparations, the sample is often in the form of a racemic (50:50) mixture, and the production of the pure substances for isotherm determination represents a significant experimental effort. The equilibrium theory of chromatography is an elegant tool for the analysis of the dynamics of both fixed-bed and moving bed chromatography [4,19,20,21]. The theory originally developed for the classical Langmuir adsorption isotherm has been recently extended to a generalized class of Langmuir isotherms [22]. The common application of the equilibrium theory is for “normal problems” that concern the calculation of band profiles of disturbances introduced at the column inlet, e.g., pulse, step, etc., for case where

345

the adsorption isotherm parameters are known. Estimating the isotherm parameters from the measured elution profile represents an “inverse problem”. Methods such as the competitive frontal analysis and elution by characteristic points solve the inverse problem. The inverse method which is increasingly used in the recent times is also a way of solving such problems [9]. In short, using a set of assumed isotherm and mass transfer parameters the elution profile is numerically calculated using a suitable chromatographic model and compared with the experimentally measured one. The parameters are then adjusted using an optimization routine so as to minimize the difference between the experimentally measured and calculated band profiles. This procedure allows for the incorporation of deviations from equilibrium theory, such as incorporation of mass transfer effects and axial dispersion. However, the determination of the parameters using the inverse method is time consuming as it requires the numerical solution of the partial differential equations for every set of parameters over the search space of the decision variables. Further, good initial guesses and range of values that the decision variables are key aspects that determine the speed of convergence of the inverse method. This paper develops the “binary retention time method (BRTM)” to determine competitive Langmuir isotherm parameters, a case very often encountered in practical situations. The method can be considered as an advancement of the retention time method, developed by Golshan-Shirazi and Guiochon [18], to the case of binary mixtures. It uses only the retention times of the two shock fronts resulting from the pulse injection of binary mixtures, and the Henry constants which can be easily measured through injection of dilute mixtures. Since the method is developed for pulse injections of binary mixtures, it eliminates the need to prepare pure components and uses minimal amount of the solute. Further, since only the shock retention times are measured, no detector calibration is needed. These features are particularly attractive in real world situations where using very small quantities of binary solutes, it is essential to evaluate the performance of chromatography as a potential route for scaling up the production of pharmaceutical compounds [23]. The paper is organized in the following manner. The local equilibrium theory of binary chromatography for a competitive Langmuir isotherm, particularly in perspective of the chromatographic cycle, is introduced. Then, the key aspects of the experimental protocol is developed. Finally, the BRTM is applied to examples from the literature and the estimated parameters are compared to those evaluated by other techniques. 2. Chromatographic cycle The chromatographic cycle consists of describing the movement of solute bands along a column for a situation when a binary mixture of known concentration is injected at the inlet of a clean column for a finite amount of time 0 . After 0 , the injection of the solute is stopped and the solvent continues to flow, thereby eluting the pulse. Mathematically, the initial and boundary conditions can be written in the following fashion at  = 0,

0≤x≤1:

at x = 0,

0 <  < 0 :

at x = 0,

 ≥ 0 :

c1 = 0, c2 = 0

(state O)

c1 = c10 , c2 = c20

c1 = 0, c2 = 0

(state P)

(state O)

(1) (2) (3)

The solution to this problem can be obtained by considering the problem as a combination of a saturation of an initially clean column, followed by elution. This problem has been addressed, within the frame of the equilibrium theory, by several authors in the past [19,24,25]. Appendix A reviews the basics of the equilibrium theory for the case of competitive Langmuir isotherms that is necessary for the discussion of the present work. In this section, the solu-

346

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

tion scheme developed by Rhee et al. [19] is reviewed and some important observations are made. Fig. 1 shows the development of a chromatographic cycle for an arbitrary set of parameters . Fig. 1(a) and (b) shows the solution of the chromatographic cycle on the (ω1 , ω2 ) and on the hodograph planes, respectively. In the figures, P and O denote the feed and initial states, respectively, while R and Q represent the two intermediate states. The chromatographic cycle consists of two simple wave transitions (1 , 2 ) and two shock transitions ( 1 , 2 ). The solution of the problem is shown in the hodograph plane as thick lines where the image of the path is given by 2

1



2 1 →O O→R→P→Q

(4)

It is worth noting that once the initial and final states are specified, the solution on the hodograph plane and the (ω1 , ω2 ) plane depends only on the isotherm parameters and not on the operating conditions. Once the solution is constructed on the hodograph plane, it can be transferred to the physical plane through the use of Eqs. (A.17) and (A.20). The solution of the problem on the physical plane shown in Fig. 1(c). As discussed before, the two shock fronts denoted as S1 and S2 originate at (x = 0,  = 0) corresponding to the 1 and

2 shocks, respectively. Also, the two simple wave transitions, 1 and 2 , represented by the C1 and C2 characteristics, respectively, originate at x = 0,  = 0 . For short injection times, which is very often the case in many practical situations, the shock fronts and waves interact within the chromatographic column and results in the chromatogram that is measured by a detector at the column outlet. It is worth noting that this case provides the solution for the general case. By allowing 0 to be large enough, the solutions of the saturation and desorption problems can be determined. Using the trajectories of the transitions on the physical plane, the elution profile can be computed. One such example evaluated at x = x is shown in Fig. 1(d). In this study, we use only property of shocks and hence the details of other interactions, that can anyhow be found elsewhere [19], are not discussed. Let us first discuss the trajectory of the 2 shock S2 depicted in Fig. 1(c). From the origin till the point B, the 2 shock travels with a constant velocity, a situation similar to the saturation of the column. Between the points B and D, the shock interacts with the 1 characteristics. Since the shock and wave belong to different families, the interaction decelerates the shock and the waves are refracted. It is worth noting that concentration of c1 jumps across the 2 shock. From the point D, the shock moves at a constant velocity until it meets a 2 characteristic at G. Since the shock and

Fig. 1. Development of a chromatographic cycle for an arbitrary system (a) Image of the solution on the (ω1 , ω2 ) plane. (b) Image of the solution on the hodograph plane. (c) Plot of the solution on the physical plane. (d) Calculated elution profile at x = x .

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354 Table 1 Coordinates and equations of key points and trajectories of the two shocks on the physical plane. These equations are summarized from the literature [19]. Point B

xB =

Point D

xD =

Point G

xG =



0 ω0 ω0 (1−ω0 ) 1 2 1

;

0  2 (1−ω0 ) 1

;

(−1)2 ω0 ω0 1 2

0 (1−ω0 ) 1

2

xE =

S2 shock OB

 = (1 +

ω10 ω20 )x

S2 shock BD

=

ω20 x −

S1 shock OE

S1 beyond E

= =

1+



E =

;

 

  

S2 shock beyond G

0 1−ω0

ω

2



 x 

1+

 0

−

x(ω1 ) = xE +



1+

 2 (1−ω0 )

+



1 (−1)2 ω0 ω0 1 2

[ω0 ]

2



2

xG



0 (ω0 −ω0 )(1+ω0 /) 2 1 1 ω0 ω0 (1−ω0 )(1−ω0 ) 1 2 1 1

 2 1 ω0

0

−1

+ 0 + x

1

0 (1−ω0 )(1−ω0 ) 1

2 + 0 + x

2

(−1)ω0 ω0 1 2

 x



1 (−1)2 ω0 1

(x − xD ) + D



ω0 1

1

1−ω0

2

Point E





1 ω0 ω0 1 2

G = 0 +

;

0 (ω0 −ω0 ) 2 1 ω0 ω0 (1−ω0 )(1−ω0 ) 1 2 1 1

=



1+

D = 0 1 +

(−1)ω0 ω0 (1−ω0 ) 1

S2 shock DG

B =

⎡

 (1−ω1 )

⎢ (ω ¯ ) 1 ⎣ ω11 + 1−ω 1

⎤

ω1 (ω ¯ 1) [ω1 ]2

⎥

dω1 ⎦

ω0 1

2

(ω1 ) = E + (ω ¯ 1 ) + (1 + [ω1 ] )(x(ω1 ) − xE ) 2

(ω ¯ 1 ) = [1 + [ω1 ] ]xE +

0 (1−ω0 )(ω0 −1)(ω1 )2 1

1

ω10

2

(ω0 −1)ω0 ω0 1

− E + 0

2

≤ ω1 ≤ 1/

the wave belong to the same family, they cancel each other and the shock is further decelerated. Now, let us consider the propagation of the 1 shock, S1 , which originates at O. The shock moves at a constant velocity until it meets the refracted 1 wave at E. Since the wave and the shock belong to the same family, they cancel each other and the 1 shock starts to decelerate. The equations of various transitions and those of key points on the (x, ) plane is given in Table 1. As can be seen from the table, the trajectories of the two shocks on the physical plane are piecewise continuous and can be computed once the isotherm parameters and operating conditions are specified. It is worth noting that the location of the points B, D, G and E and the entire trajectory of the 2 shock are all related to 0 , ,  and to the parameters ω10 and ω20 , which correspond to the feed concentrations alone! While the trajectory of the 2 shock is obtained in a closed form, that of the 1 shock beyond the point E is obtained only in a parametric form requiring numerical integration. Further, similar to 2 , the trajectory of the 1 shock is related to 0 , ,  and to the parameters ω10 and ω20 which correspond to the feed concentration and the operating conditions. For a normal problem when the operating conditions and isotherm parameters are available, the chromatographic cycle can be first solved on the hodograph plane and then on the physical plane as described above and elution profiles can be conveniently obtained. 3. Isotherm parameter evaluation All the above discussions concern the “normal problem”, i.e., those which allow for the description of the chromatographic cycle for a set of isotherm parameters. However, the determination of isotherm parameters from a set of experimentally measured chromatograms is an “inverse problem”. The solution of “inverse problems” is particularly complicated for the following reasons. Firstly, the measured elution profile is recording of the detector

347

output as a function of time. Hence, it is not possible to transform the information about the concentrations to the hodograph plane. The conversion to concentration units can be performed either by calibrating a detector, or by collecting samples for offline analysis. In both cases, the calibration does constitute an extra experimental step. Secondly, in many cases, the detector used for the experiments may not be able to distinguish between the two solutes, e.g., a UV detector used for chiral chromatography. Hence, it might not be possible to obtain the concentration profile of the individual solutes unless a suitable detector or offline analyses are involved. Finally, obtaining isotherm parameters from the measured retention time of the shocks, which is the objective of this work, is not straightforward for the reasons explained below. It can be seen from Table 1 that the equations to describe the shock fronts are described by variables that relate to the experimental system and operating conditions such as 0 ,  and  and the parameters related to the isotherm, ω10 and ω20 . The experiexp exp mental retention time of the shocks, S1 and S2 are measured at the exit of the column, i.e., at x = 1. Hence, if the isotherm parameters are not known, the unknown variables in the equations of the shock trajectories are ω10 and ω20 . It might appear that exp exp these can be evaluated from the measurement of S1 and S2 by using appropriate equations given in Table 1. However, this is not straightforward for the following two reasons. Firstly, it is not possible to directly invert the equation describing the 1 shock trajectory beyond point E. Secondly it is important to note that the equations of both the shocks are piecewise continuous. Hence, based on the measured chromatogram it will not be possible to unambiguously choose the correct equation for the inversion. In other words, since the location of the points B, D, G and E are not known a priori, it will not be possible to determine which segment of the trajectory exits the column, i.e., crosses x = 1. In order to overcome these limitations, the inverse problem is recast as an optimization problem to evaluate the isotherm parameters. This forms the basis of the two methodologies developed below. The first methodology, termed “Generalized binary retention time method (G-BRTM)” is applicable for both mass (same injection volume with different injected concentrations) or volume (same injected concentrations with different injected volumes) overload conditions. The second methodology, termed “Volume overload binary retention time method, (V-BRTM)”, is applicable for volume overload injections only. 4. Generalized binary retention time method (G-BRTM) The G-BRTM consists of two experimental steps. The first experimental step (step 1) involves the measurement of the Henry constants of the two solutes which can be achieved by injecting a dilute binary mixture. At dilute conditions, the competition for the adsorption sites becomes negligible and the two solutes travel with a velocity independent of each other. Under these conditions, the Henry constant can be calculated from the retention time of the solute peak using the following relationship: Hi =

R,i − 1 

(5)

where R,i is the peak retention time. Step 2 involves performing overloaded pulse injections, i.e., injecting a pulse of a concentrated mixture. From the measured chromatograms, the retention times of the two shock fronts are to be measured. In practical situations, mass transfer resistances and axial dispersion lead to diffused shock fronts, and hence the retention time of the shock fronts has to be estimated. Good approximations can be obtained by identifying the inflection points of the elution profiles by computing the maxima of the gradient (prefer-

348

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

objective function, J, is evaluated: J=

n 

exp

2

exp

2

calc calc [(S1,k − S1,k ) + (S2,k − S2,k ) ]

(6)

k=1

Fig. 2. Methods of estimating shock retention times for an experimental chromatogram.

exp

able for smooth experimental data) or by finding the time, S1 , that equalizes the hatched areas A1 and A2 for the shock front of exp the lighter component and S2 , the time that equalizes the areas A3 and A4 for the shock front of the heavier component as described in Fig. 2 [11]. The competitive Langmuir isotherm has four independent parameters, H1 , H2 , N1 and N2 . While the Henry constants, H1 and H2 can be measured in step 1, the saturation capacities, N1 and N2 are to be evaluated from step 2. As discussed before, since the inverse problem is not amenable for a straightforward solution, the problem is cast as an optimization problem. In this method, the optimizer chooses values for the decision variables, N1 and N2 , within a range specified by the user. Using the measured values of H1 and H2 , the assumed values of N1 and N2 , along with other operating parameters,  and , the values of ω10 and ω20 can be calculated from their appropriate definitions. The location of the points B, D, G and E on the physical plane can hence be calculated using the expressions listed in Table 1. Once these points are identified, the suitable equation of the shock trajectory that exits the column, i.e., x = 1, can be identified using the criteria listed in Table 2. The evaluation of the retention time of the 2 shock is rather straightforward. Upon identification of the suitable segment calc can be obtained by evaluatof the trajectory, the retention time S2 ing the expression at x = 1. For the 1 shock, should xE ≥ 1, then the trajectory is described by the equation of segment OE from which calc can be computed at x = 1. However, if x ≤ 1, then the trajecS1 E calc tory of the shock is obtained in a parametric form and hence S1 has to be obtained numerically. This can be performed in two ways. The entire trajectory beyond point E can be obtained by calculating x(ω1 ) and (ω1 ) for discrete values of ω1 with ω10 ≤ ω1 ≤ 1/. Upon computation of the trajectory, the retention time corresponding to x = 1 can be evaluated by interpolation. Alternatively, the equation describing x(ω1 ) can be solved numerically for ω1 (x = 1) and subsequently (ω1 ) can be evaluated by using ω1 (x = 1) to obtain calc . These calculated retention times,  calc and  calc , are then comS1 S1 S2 exp exp pared with the experimentally measured ones S1 and S2 and the

Table 2 Criteria to determine the choice of piecewise continuous equation to evaluate retention time of the shock. The equations of different segments are given in Table 1. Shock S2

S1

Condition

Segment of shock trajectory that crosses x = 1

xB ≥ 1

OB

xB < 1 ≤ x D

BD

xD < 1 ≤ xG

DG

xG < 1

Beyond G

xE ≥ 1

OE

xE < 1

Beyond E

The above equation denotes the sum of squared errors between the calculated and measured shock retention times of n experiments. The optimizer continues to search for values of decision variables until the minimum value of J is obtained. A common issue in optimization methods is the convergence at local minima. In order to reduce the risk of convergence, non-gradient based schemes such as genetic algorithms can be used [26]. However, it is worth mentioning that numerical results, not reported here, indicate that the surface of J as a function of N1 and N2 , indicate the presence of a unique minima. 5. Volume overload binary retention time method (V-BRTM) The G-BRTM can be used without limitations on the type of injections. However, from a practical point of view, especially in industrial applications, where the time to perform isotherm measurements are rather limited, volume overloaded injections seem more convenient. Firstly, there is no necessity to prepare additional samples with varying concentrations. Secondly, the injection of samples at the maximum concentration (as dictated by solubility) using different injection loops can be readily carried out. In the following a two-step methodology, henceforth called as “Volume-overload binary retention time method (V-BRTM)” is proposed which enables rapid determination of competitive Langmuir isotherms for volume overloaded injections. This first step in the V-BRTM, similar to the G-BRTM (step 1) described earlier, requires the injection of dilute mixtures to measure Henry constants. The second step (step 2) involves multiple injections of samples with identical concentrations but with different injection volumes. All other operating parameters are to be kept constant. As will be shown shortly, this experimental procedure reduces the complexity in parameter estimation and makes way for rapid determination of parameters. 5.1. V-BRTM for the general case of h = N1 /N2 = / 1 The following analysis is developed for a general case, where the saturation capacities of the two solutes are allowed to be nonidentical. As described, earlier the solution of the chromatographic cycle is conveniently represented in the hodograph plane and on the (ω1 , ω2 ) plane. Since the injected concentration is fixed for a set of volume overload experiments, the solution of the problem on the hodograph plane and on the (ω1 , ω2 ) plane will be identical for all the injections, i.e., the image of the solution on these planes are identical. It can be further seen from Table 1 that all the characteristic points and equation of the shock trajectories is described by the parameters ω10 and ω20 , which correspond to the injected concentration of the solutes. Hence, one can foresee that for such experiments, instead of using N1 and N2 as decision variables, it is more suitable to use ω10 and ω20 as decision variables. It is particularly attractive since these variables, unlike the saturation capacities, are bound within a finite range, as defined by Eq. (A.13). It is thus possible to generalize the solution to any arbitrary system and no information about the range of saturation capacities etc. need to be provided by the user. The optimization problem is hence to minimize the objective function defined by Eq. (6) with ω10 and ω20 being the decision variables. Once ω10 and ω20 are determined by the optimizer as described earlier, the Langmuir isotherm parameters can be obtained by solving Eqs. (A.3), (A.14) and (A.15).

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

5.2. V-BRTM for the special case of h = N1 /N2 = 1

should lie in the range:

The analysis developed in the previous sub-section is applicable for a general class of systems where there is no constraint on the equality of N1 and N2 . However, to be thermodynamically consistent, it is often important to evaluate the isotherm parameters with the requirement of N1 = N2 . In this case, when the Henry constants of the two solutes are fixed, of the two isotherm parameters, N1 and N2 , we have only one independent parameter. In terms of the transformed variables, this indicates that ω1 and ω2 are no longer independent. The relationship between ω1 and ω2 can be obtained by solving Eqs. (A.6) and (A.7) with the condition h = N1 /N2 = 1: ω1 =

ω2 ( + ) − − 1 ω2 ( +  2 ) − − 

349

(7)

where = c10 /c20 . It is worth recollecting that the ranges of values that ω10 and ω20 can take is give by Eq. (A.13). In order to satisfy the range of values that ω1 can take, i.e., 0 < ω1 ≤ 1/, the values of ω20

1+ ≤ ω20 ≤ 1  +

(8)

The problem is now considerably reduced to finding the value of ω20 which minimizes the objective function given by Eq. (6). Once ω20 is determined, the corresponding value of ω10 can be computed from Eq. (7) and the isotherm parameters can be evaluated by solving Eqs. (A.3), (A.14) and (A.15). 6. Optimization algorithm The optimization problem is to minimize the objective function given in Eq. (6). The choice of decision variables depends on the choice of methods adopted for the solution. In the case of the GBRTM, N1 and N2 are used as decision variables. In the case of the V-BRTM, ω10 and ω20 are used as decision variables for the general case, while ω20 is used for the special case. A suitable optimizer that can solve such problems can be used. The algorithm developed for the case studies considered in this study is discussed next.

Fig. 3. Experimental pulse response (symbols) and calculated profiles from equilibrium theory (lines) of a racemic mixture of 1-phenyl-1-propanol on Chiralcel OD under SFC conditions. Dotted lines indicate the location of the shock fronts. Total injection concentration of the racemic mixture: (a) 50 g/L, (b) 100 g/L, (c) 250 g/L, and (d) 500 g/L. Experimental results are from Ottiger et al. [14]. Other details are given in Table 3.

350

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354 Table 3 Operating conditions of different systems considered for the case study. Case

Parameter

Case 1

Column length, L (cm) Porosity, b Column hold-up time, L/v (min) Number of plates Injected volume, V (␮L) Henry constant, Hi

Case 2

R enantiomer

S enantiomer 25.00 0.722 2.722

2380 20.50 4.20

Column length, L (cm) Porosity, b Column hold-up time, L/v (min) Number of plates Injected concentration, ci0 (g/L) Henry constant, Hi

2090 20.50 5.28 25.00 0.77 3.20

8000 2.015 7.24

8200 2.068 6.04

As discussed earlier, the first step in the protocol is the estimation of the Henry constants, which for the case studies considered were obtained from the reported values. The optimization was performed in MATLAB® . For both the methods, G-BRTM and V-BRTM, a constrained non-linear multivariable solver, fmincon, was used. For the special case of V-BRTM with h = N1 /N2 = 1, fminbnd, an in-built solver developed for single-variable functions on fixed intervals was found to yield faster convergence compared to fmincon. The solver fmincon requires lower and upper bounds for the decision variables and suitable initial guesses. For all G-BRTM solutions, values of 0.01 g/L and 200 g/L were used as the lower and upper bounds for the saturation capacities while an initial guess of 50 g/L was used. Various values within the range of the lower and upper bounds were used and all of them yielded identical results, indicating that local minima were perhaps not present. However, the time required to obtain the solution for different sets of initial conditions varied marginally. For the case of the V-BRTM, the bounds for the decision variables were defined by Eqs. (A.13) and (8) for the general case and the special case, respectively. In all the cases reported below, the run time required for the estimation was less than 10 s on a standard desktop computer. Fig. 4. Non-linear isotherms for two enantiomers of 1-phenyl-1-propanol on Chiralcel OD under SFC conditions [14]. Symbols: isotherm using parameters from literature [14]. Solid lines: isotherm using parameters obtained from G-BRTM.

7. Case studies In the previous sections the BRTM was developed. In order to evaluate the efficiency of the method, a case study using two cases from the literature are considered. The first concerns mass overload injections on a column with moderate efficiency, while the second concerns volume overload injections on a column with high efficiency. The examples considered are enantiomeric systems. Although, more complex isotherms, e.g., bi-Langmuir, are used to describe the adsorption of enantiomers on chiral stationary phases, several examples exist in the literature where the simple competitive Langmuir model seems to adequately describe the adsorption equilibria [14,27]. The primary objective of the case studies is to compare the isotherm parameters obtained by the BRTM with a more rigorous technique, e.g., the inverse method. It is worth emphasizing that in doing so, we implicitly assume that the competitive Langmuir

model describes the adsorption equilibria adequately. 7.1. Case 1: mass overload injections The first case study considers mass overload injections of 1phenyl-1-propanol on a Chiralcel-OD column under supercritical fluid chromatography (SFC) conditions [14]. Supercritical carbon dioxide containing methanol as a modifier was used as a mobile phase. Relevant experimental data are listed in Table 3. As can be seen from the table, this case corresponds to a situation where the column had a moderate efficiency (number of plates ≈ 2200) [28]. In order to illustrate the efficiency of the BRTM, a specific set of experimentally measured chromatograms, shown in Fig. 3, are considered. The operating conditions for these chromatograms are

Table 4 Comparison of isotherm parameters from literature and those from the binary retention time method for the cases studied. Case

Solute

Inverse method

G-BRTM

V-BRTM

Hi [–]

Ni [g/L]

Hi [–]

Ni [g/L]

1

R enantiomer S enantiomer

4.20 5.28

113.89 88.62

4.19 5.28

135.31 86.07

2

R enantiomer S enantiomer

7.24 6.04

79.29 78.90

7.39 6.12

69.44 72.20

Hi [–]

Ni [g/L]

7.39 6.12

70.82 71.33

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

351

Fig. 5. Experimental pulse response (symbols) and calculated profiles from equilibrium theory (lines) of a mixture of S-CPP (2.068 g/L) and R-CPP (2.015 g/L) on Chiralcel OB-H. Dotted lines indicate the location of the shock fronts. Injection volume: (a) 0.25 mL, (b) 0.50 mL, (c) 0.75 mL, and (d) 1.0 mL. Experimental profiles were extracted from Cherrak et al. [29]. Other details are given in Table 3.

the following: Pressure = 170 bar; Temperature = 30 ◦ C, Modifier composition = 4.9% (w/w). In the literature, the inverse method was used to estimate the Langmuir parameters that are listed in Table 4. The set of experiments correspond to mass overload conditions and hence only the G-BRTM can be used. The Henry constants were obtained from experiments performed with dilute injections [28]. The retention time of the shocks were calculated as shown in Fig. 3. The Langmuir isotherm parameters estimated using the G-BRTM by considering all the four chromatograms are listed in Table 4. It can be seen that the estimation of the saturation capacity of the R enantiomer is larger than that obtained from the inverse method, while that of the S enantiomer is comparable. Fig. 4 shows the comparison of the isotherms obtained by the two methods. It can be seen that the deviation between the isotherms obtained from the BRTM and inverse methods are minor. This is further illustrated in the comparison of the experimental elution profiles with those obtained from equilibrium theory as shown in Fig. 3. Given the basis of the equilibrium theory, which does not account for mass transfer and axial dispersion effects, the predictions of the experimental profiles are good. It is further worth noting that the plateau on the heav-

ier component in Fig. 3(c) and the “bump” on the tail of the heavier component at  ≈ 2.5 is clearly predicted by the calculations. These clearly highlight the fact that in this case the adsorption of the two components is well represented by the competitive Langmuir isotherm and the parameters estimated by the BRTM are indeed sound. Fig. 3(d) shows an extreme case of an elution profile where the peak corresponding to the 2 shock is just visible. In this case the retention time of the shock can still be estimated and used for isotherm parameter estimation. Should there be an injection in which the shock front of the stronger component is not clearly visible, it will not be possible to apply the BRTM. This can be considered as a shortcoming of the BRTM. However, by lowering the injection volume, it will be possible to obtain experimental chromatograms where the retention time of 2 shock can be clearly identified. 7.2. Case 2: volume overload injections The second case study considers the enantioseparation of 3chloro-1-phenyl-propanol (CPP) on a Chiralcel OB-H column, the

352

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

8. Conclusions

Fig. 6. Non-linear isotherms for two enantiomers of 3-chloro-1-phenyl-propanol on Chiralcel OB-H [29]. Symbols: isotherm using parameters from literature [29]. Dashed lines: isotherm using parameters obtained from V-BRTM. Solid lines: isotherm using parameters obtained from G-BRTM. The isotherms obtained using V-BRTM and G-BRTM are indistinguishable.

operating conditions of which given in Table 3 [29]. It can be seen that the column showed a rather high efficiency (Number of plates ≈ 8000) as compared to the previous case. Cherrak et al. used both frontal analysis and the inverse method to evaluate isotherm parameters that are listed in Table 4 [29]. The inverse method used a set of four volume overload injections shown in Fig. 5. The isotherm parameters were estimated from all the injections by allowing the two enantiomers to have different saturation capacities. Since the experiments correspond to a set of volume overload injections, both G-BRTM and V-BRTM can be used. The retention times of the two shock fronts were extracted as indicated in Fig. 5 and the reported Henry constants were used for the isotherm parameter estimation using all the four experiments together.1 The parameters obtained by BRTM are compared to those obtained from the inverse method in Table 4 and the corresponding isotherms are plotted in Fig. 6. Though there is a minor deviation in the saturation capacities estimated by the BRTM and the inverse method, the isotherms are rather indistinguishable. It can be noticed that the parameters estimated by G-BRTM and V-BRTM are identical which is worth emphasizing since the V-BRTM method requires no user input ranges for the decision variables. Using the estimated parameters, the elution profiles were calculated using the equilibrium theory and are compared with the experimental ones in Fig. 6. It can be seen that although the column had high efficiency, the comparison between the predicted and experiments show deviations. For example, the plateau on the tail of the heavier component are not seen in the experiments. In fact, a similar trend was observed in the literature [29]. In such situations, a more sophisticated isotherm model is required. However, it is worth emphasizing that the results of the inverse method and BRTM are practically similar.

1 It is worth noting that in the inverse method used by Cherrak et al., all the four isotherm parameters were used as decision variables and hence the Henry constants obtained from the inverse method are slightly different from those obtained from dilute injections. For the BRTM, as described in the protocol, we used the Henry constant measured from dilute injections, as reported in Table 1 of Cherrak et al. [29].

The binary retention time method for the measurement of competitive Langmuir isotherm parameters was developed. The method uses the measured retention times of the two shock fronts of binary injections and the Henry constants to estimate the isotherm parameters. This is achieved by comparing the experimental shock retention times with the calculated values from the equilibrium theory of chromatography. The inverse problem of parameter estimation was formulated as an optimization problem and solved using a standard optimization tool. Two schemes, G-BRTM and VBRTM were introduced. The G-BRTM can be used for both mass and volume overload conditions, while the V-BRTM can be used only for volume overload injections. However, the V-BRTM is expected to be more robust as no bounds for the ranges of the decision variables are required and hence is independent of user input. Similar to any other isotherm determination method, the BRTM has its limitations. Firstly, the main limitation is that this method, as described in this work, is applicable strictly to systems that obey the competitive Langmuir isotherm. Applying this method to systems that do not obey the Langmuir isotherm may lead to spurious parameters which could lack physical significance. Secondly, unlike the inverse method which can also be used to measure mass transfer parameters, the BRTM method yields only isotherm parameters. Thirdly as seen in the case studies, the isotherm parameters obtained by the BRTM are in good agreement to those obtained by the inverse method especially when the column is efficient. When column efficiency is low, minor deviations are seen. Several applications of the BRTM can be envisaged. In the recent times, stationary phases with smaller particle sizes less than 10 ␮m are being increasingly used for both analytical and preparative separations, thus resulting in very high efficiencies. In these columns, the influence of mass transfer will be reduced and hence the assumption of local equilibrium can be approached. Under these conditions, the BRTM can be effectively used. Secondly, even for columns with moderate efficiencies, should the user prefer to use the inverse method, the BRTM can be used to obtain good initial guesses for the adsorption isotherm parameters. This can significantly reduce computational time invested in inverse methods. Thirdly, since the BRTM uses only retention times, it is independent of detector calibration and hence helps avoiding an extra experimental step. Finally, since the time required to obtain isotherm parameters are in the order of a few seconds, on a desktop computer, this method can be attractive for rapid estimation of isotherm parameters and can be envisaged to be used for making quick estimates of the process performance at large-scales. Software A free copy of the implementation of the BRTM in the form of a MATLAB code, to be used for non-commercial purposes, can be obtained by writing to the corresponding author. Acknowledgements Financial support from NTU AcRF grant RG38/05 is acknowledged. Initial discussions with Abhijit Tarafder, that motivated the work and the helpful comments of Marco Mazzotti, both from ETH Zurich are appreciated. Appendix A. Brief review of equilibrium theory In this section, for the sake of completeness the equilibrium theory of chromatography, as developed by Rhee et al. [19], is briefly reviewed. The discussion is limited to aspects that concern the

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

development of the BRTM. For a comprehensive discussion of the theory, the reader is referred to more detailed works [3,19,20,22]. Let us consider the transport of two solutes, 1 and 2, along the axial direction, z, of a chromatographic column of length L and void fraction of b . The solute is assumed to be present in a nonadsorbing solvent whose interstitial velocity is v. Within the frame of the equilibrium theory of chromatography, the component mass balance of the solutes can be written as ∂ ∂c (ci + ni ) + i = 0 ∂ ∂x

(i = 1, 2)

(A.1)

where ci is the concentration of solute i in the fluid phase, ni is the concentration on the stationary phase that is in equilibrium with the fluid phase concentration ci and,  = (1 − b )/b is the phase ratio. In the above equation  = t v/L and x = z/L represent the dimensionless time and axial coordinates, respectively. If the adsorption of the two solutes can be represented by the competitive Langmuir adsorption isotherm, the equilibrium relationship is given by Hi ci 1 + K1 c1 + K2 c2

ni =

(A.2)

Hi = Ni Ki

(A.3)

where Ni and Ki are the saturation capacity of the stationary phase and the equilibrium constant of solute i, respectively. The Henry constant, Hi , represents the slope of the isotherm at infinite dilution. For the isotherm to the thermodynamically consistent the saturation capacities of the two components should be the same, i.e., N1 = N2 [1]. However, for many practical cases, experimental data seem to be better described by relaxing this requirement. Having said this, it should be mentioned that the equilibrium theory of chromatography allows for the development of the solution for the general case when N1 = / N2 . The selectivity, , is defined as  = H2 /H1 with H2 ≥ H1 . Eqs. (A.1) and (A.2) together represent set of first-order, homogenous, and reducible partial differential equations. They can be solved with the definition of a suitable initial value problem. Seminal contributions have been made concerning the solution of Eq. (A.1) for problems of practical relevance [19,30,31]. The method of characteristics is one such method in which the solution of Eq. (A.1) is conveniently solved on the hodograph plane (c1 , c2 ) plane. The characteristics on the hodograph plane is obtained from the solution of the following fundamental differential equation: ∂n2 ∂c1

 dc 2 1

dc2



−

∂n2 ∂n1 − ∂c1 ∂c2



dc1 dc2



−

(A.4)

c2  2 − (k + c1 − hc2 ) − hc1 = 0

(A.5)

where  = dc1 /dc2 , h = N1 /N2 and k = ( − 1)/(K1 ). Eq. (A.5) has two real roots with 1 and 2 being the negative and positive roots, respectively: 1 = 2 =

2c2

 1

2c2

(k + c1 − hc2 ) − (k + c1 − hc2 ) +

 

2

(k + c1 − hc2 ) + 4hc1 c2 2

(k + c1 − hc2 ) + 4hc1 c2



c1 = 2 c2 −

k2 2 + h

−h < 1 ≤ 0 ≤ 2 < ∞

(A.10)

ω1 =

1 + h 1 + h

(A.11)

ω2 =

2 + h 2 + h

(A.12)

The parameter ω1 is constant along 2 characteristic while ω2 is constant along the 1 characteristic. The limits of the two parameters follow from Eq. (A.10) as 0 < ω1 ≤

1 ≤ ω2 ≤ 1 

(A.13)

It is worth noting that compared to the pair (1 , 2 ), the transformed variables (ω1 , ω2 ) are bound within finite limits. A direct mapping between the parameters (ω1 , ω2 ) and the concentrations of the two solutes (c1 , c2 ) exists and is given by c1 =

(1 − ω1 )(ω2 − 1) K1 ( − 1)ω1 ω2

(A.14)

c2 =

(1 − ω1 )(1 − ω2 ) K2 ( − 1)ω1 ω2

(A.15)

The fundamental relationships described above have been used to study Riemann problems, where the initial data is piecewise constant. These problems involve the specification of an initial and a final state that are connected by a discontinuity. The transition between the two states on the hodograph plane occurs through an intermediate state, which is obtained by the intersection of the 1 and 2 characteristics that originate from the initial and final states, respectively. These states are separated by two transitions. The first transition separating the initial and intermediate states and the second transition separating the intermediate and final states, respectively. It can be shown that while moving from left to right on the physical plane, the transitions connecting the initial state and the final state is given by 1

Final state→Intermediate state→Inital state

(A.16)

If moving from left to right, on the physical plane, yields a family of characteristics that fan clockwise, we obtain continuous simple wave transitions, for which, i consists of Cj characteristics whose slopes on the physical plane are given by j =

d = 1 + j dx

(j = 1, 2)

(A.17)

where j is the eigen value of the matrix N = [nij ], with nij = ∂ni /∂ci . It can be shown that

(A.7)

and

(A.8)

(A.9)

A useful variable transformation can be performed to introduce two new parameters ω1 and ω2 which are defined as

1 = ω12 ω2



k1 1 + h

with

(A.6)

Each root yields a family of characteristics that are straight lines on the hodograph plane. One family, denoted by 1 , has a positive slope, while the other, denoted by 2 , has negative slopes. The equation of the characteristics on the hodograph plane can be expressed in terms of parameters, 1 and 2 . 1 characteristic :

c1 = 1 c2 −

2

∂n1 =0 ∂c2

In the case of the competitive Langmuir isotherm, the above equation reduces to

 1

2 characteristic :

353

(A.18)

1 = ω1 ω22

(A.19)

If the characteristics do not fan clockwise, then the initial and final states are separated by discontinuities referred to as shocks (denoted by j ). Under such situations the slope of the characteristic on the physical plane is given by j =

d = 1 + ˜ j dx

(j = 1, 2)

(A.20)

354

A. Rajendran, W. Chen / Separation and Purification Technology 67 (2009) 344–354

where ˜ j =

 [n ]  1

[c1 ]

=

j

 [n ]  2

[c2 ]

(A.21)

j

The terms within the square brackets indicates the jump in the respective quantity across the discontinuity. In such situations, the transition separating the initial and final states, as represented in Eq. (A.16), consists of two discontinuities S1 and S2 corresponding to 1 and 2 characteristics, respectively. It can be shown that ˜ 1 = ω1L ω1R ω2

(A.22)

and ˜ 2 = ω1 ω2L ω2R

(A.23)

where the superscripts L and R refer to the states on the left and right of the characteristic, respectively. In the case of systems following the Langmuir isotherm, since the characteristics on the hodograph plane are straight lines, it is possible to show that, on the hodograph plane, the images of 1 and 2 coincide with 1 and 2 , respectively. References [1] D.M. Ruthven, Principles of Adsorption and Adsorption Processes, John Wiley, New York, 1984. [2] D.M. Ruthven, C.B. Ching, Chem. Eng. Sci. 44 (5) (1989) 1011–1038. [3] G. Guiochon, S.G. Shirazi, A.M. Katti, Fundamentals of Preparative and Nonlinear Chromatography, Academic Press, Boston, 1994. [4] M. Mazzotti, G. Storti, M. Morbidelli, J. Chromatogr. A 769 (1) (1997) 3–24. [5] M. Juza, M. Mazzotti, M. Morbidelli, Trends Biotechnol. 18 (3) (2000) 108–118.

[6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

G. Guiochon, J. Chromatogr. A 965 (1–2) (2002) 129–161. A. Rajendran, G. Paredes, M. Mazzotti, J. Chromatogr. A 1216 (4) (2009) 709–738. A. Seidel-Morgenstern, J. Chromatogr. A 1037 (1–2) (2004) 255–272. A. Felinger, A. Cavazzini, G. Guiochon, J. Chromatogr. A 986 (2) (2003) 207–225. E.V. Dose, S. Jacobson, G. Guiochon, Anal. Chem. 63 (8) (1991) 833–839. O. Lisec, P. Hugo, A. Seidel-Morgenstern, J. Chromatogr. A 908 (1–2) (2001) 19–34. C. Blümel, P. Hugo, A. Seidel-Morgenstern, J. Chromatogr. A 865 (1–2) (1999) 51–71. L. Asnin, G. Guiochon, J. Chromatogr. A 1089 (1–2) (2005) 105–110. S. Ottiger, J. Kluge, A. Rajendran, M. Mazzotti, J. Chromatogr. A 1162 (1) (2007) 74–82. G. Yang, Chromatographia 50 (9–10) (1999) 621–624. D. Burger, R. Neumueller, G.L. Yang, H. Engelhardt, I. Quinones, G. Guiochon, Chromatographia 51 (9–10) (2000) 517–525. T.-W. Chen, N.G. Pinto, L. Van Brocklin, J. Chromatogr. 484 (1989) 167–185. S. Golshan-Shirazi, G. Guiochon, Anal. Chem. 61 (5) (1989) 462–467. H.K. Rhee, R. Aris, N.R. Amundson, First-order Partial Differential Equations, vol. II, Dover Publications, Mineola, NY, 2001. H.K. Rhee, R. Aris, N.R. Amundson, First-order Partial Differential Equations, vol. I, Dover Publications, Mineola, NY, 2001. G. Storti, M. Mazzotti, M. Morbidelli, S. Carra, AIChE J. 39 (3) (1993) 471–492. M. Mazzotti, Ind. Eng. Chem. Res. 45 (15) (2006) 5332–5350. C.J. Welch, P. Sajonz, G. Spencer, W. Leonard, D. Henderson, W. Schafer, F. Bernardoni, Org. Process Res. Dev. 12 (4) (2008) 674–677. S. Golshan-Shirazi, G. Guiochon, J. Phys. Chem. A 93 (10) (1989) 4143–4157. M. Mazzotti, J. Chromatogr. A 1126 (1–2) (2006) 311–322. Z. Zhang, M. Mazzotti, M. Morbidelli, J. Chromatogr. A 989 (1) (2003) 95–108. S. Khattabi, D.E. Cherrak, J. Fischer, P. Jandera, G. Guiochon, J. Chromatogr. A 877 (1–2) (2000) 95–107. A. Rajendran, M. Mazzotti, M. Morbidelli, J. Chromatogr. A 1076 (1–2) (2005) 183–188. D.E. Cherrak, S. Khattabi, G. Guiochon, J. Chromatogr. A 877 (1–2) (2000) 109–122. E. Glueckauf, Discussions of the Faraday Society 7 (1949) 12. F. Helfferich, G. Klein, Multicomponent Chromatography, in: Theory of Interference, Marcel Dekker, New York, 1970.