Modelling and simulation of multicomponent nonlinear chromatography

Computerschem. Engng Vol. 22, No. I-2, pp. 259-281. 1998

Pergamon

Copyright O 1997 Elsevier Science Lid Printed in Great Britain. All rights reserved

PH: S0098-1354(96)00364-X

oo98-1354/97$17.00+0.00

Modelling and simulation of multicomponent nonlinear chromatography Reto Zenh~iusern* and David W. T. Rippin EidgenOssische Technische Hochschule (ETH), Laboratorium for Technische Chemie, 8092 Ztirich, Switzerland

(Received 25 October 1994; revised 27 November 1995) Abstract The equilibrium relationship between the stationary and the mobile phase plays a leading role in chromatography; especially for higher concentrations where the equilibrium relationship becomes nonlinear. Based on local equilibrium theory, a simple numerical algorithm is developed by considering the multicomponent nonlinear equilibrium relationship and its Jacobian matrix to quickly obtain a first prediction of the chromatographic separation without the need of providing mass transfer data. The impact of several favorable muiticomponent equilibrium isotherms on chromatographic separations is studied: The multicomponent Langmuir isotherm, a polynomial and an exponential extension of the multicomponent Langmuir isotherm. © 1997 Elsevier Science Ltd

Introduction Chromatography is being used increasingly for production scale separations. It can provide high selectivity while avoiding severe operating conditions which may be damaging to sensitive molecules. To support the design and optimal operation of such separations, available theoretical models must be more effectively applied. Efficient procedures have to be developed for characterization of the multicomponent mixtures. Today, modelling of multicomponent chromatography has reached a high degree of complexity. Many factors such as micropore and macropore diffusion, temperature dependencies, and complex equilibrium relationships obtained from thermodynamics can be modelled. However, due to the freedom of choosing the operating conditions in chromatography (mobile phase, stationary phase, temperature, pressure, concentration range by choosing the dilution of the components in the feed) it is often not possible to estimate all parameters needed for the simulation. Common to all chromatographic simulations is the estimation of the equilibrium relationship which describes the concentrations on the stationary phase as a function of the concentrations in the mobile phase. Since the right choice of the stationary and the mobile phase is one of the most important steps during * Present address: TeranolAG/Lalden,Postfaeh310, 3930 Visp, Switzerland.

the design of a chromatographic separation, quick qualitative answers are needed during the initial stage of a chromatgraphic separation project. Ideal nonlinear chromatography, based on local equilibrium theory, assumes immediate equilibrium between the mobile phase and the stationary phase. Any contributions from hydrodynamic effects or mass transfer phenomena are neglected. They only have a smoothing effect on the chromatographic peaks obtained from local equilibrium theory. The resulting system of coupled hyperbolic partial differential equations belongs to the well-known family of conservation laws. Solution of these systems takes advantage of wave theory (Courant and Friedrichs, 1952; Courant et al., 1952; Courant and Friedfichs, 1976). A good overview of this subject, with more emphasis on fixed-bed processes is given by Tondeur (1987a) Tondeur and Bailly (1987b) and Schlechter et al. (1987). Other systems in process engineering such as e.g., fixed-bed sorption, fluid-fluid extraction, sedimentation and catalytic fixed-bed reactors are also found to belong to the same family (Loureiro et al., 1983; Marquardt, 1989). Numerical solution methods for such systems recently have been summarized by Loureiro and Rodrigues (1991). DeVault (1943) first formulated the qualitatively correct form of the basic equations of single and multicomponent chromatography. Based on this work Glueckauf (1949) developed an analytical solution for binary systems using the Longmuir isotherm. The

259

260

R. ZENH.~USERNand D. W. T. RIPPIN

solution can generally be found by constructing the trajectories in concentration space (Glueckauf, 1949). Analytical solutions for the continuous case (representing elution) and the discontinuous case (representing loading) have been found. Model simplification was always a main goal of local equilibrium theory. However, all early attempts suffered from either being limited to the binary Longmuir isotherm or they were only of a qualitative nature. In the remarkable monograph by Helfferich and Klein (1970) multicomponent chromatography was discussed for the first time in a rigorous manner. Of fundamental importance within this context is the concept of coherence. Coherence is considered as the state that a transient, nonequilibrium, non-steady state perturbation seeks to assume (Helfferich, 1986) i.e., a distinct concentration path is followed during sorption. The idea of regarding coherence as a steady state in addition to equilibrium is rather unfamiliar to chemists and chemical engineers. The theory of Heifferich and Klein seems to be very intuitive in large parts. Rhee et al., 1968, 1970, 1986, 1989 followed a mathematical approach. As in the work of Helfferich and Klein the theory of chromatography has also been extended to multicomponent systems. Unfortunately, only the case for the multicomponent Langmuir isotherm has been solved explicitly which is equivalent to constant separation factor systems. For this type of equilibrimn isotherm an analytical solution could be found. The influence of mass transfer on the shape of the profiles has also been discussed. It has been shown (Rhee et al., 1971; Rhee and Amundson, 1972; GolshanShirazi et al., 1989b) that the form of the equilibrium isotherm is still the major factor in determining the peak shape. Increasing influence of mass transfer will result in smoother peaks compared to the ideal shape. Based on the local equilibrium theory of Rbee, Aris and Amundson a numerical algorithm has been delveloped (Zenh~iusern, 1987, 1993) with which the key examples could be successfully reproduced. Successful experimental verification for the binary Langmuir isotherm has been recently reported (Golshan-Shirazi and Guiochon, 1988, 1989a: Katti et al., 1990). Following the trend of simplification Basmadjian and Coroyannakis, 1987Basmadjian et al., 1987a,b) presented a combined mathematical-graphical approach for binary systems. By taking advantage of the different representation in the concentration space diagram (hodograph), isotherm graph and concentration history graph, simplified rules could be derived for a graphical construction of the solution path in either diagram. Furthermore, Basmadjian's method relieves the user of solving complicated mathematics. His method is also valid for other types of isotherms e.g., solutions for unfavorable, sigmoidal and even azeotropic systems are presented. The aim of this work is to show how the theory of local equilibrium theory could support experiments to aid the design of preparative chromatographic systems. The fundamental mathematical equations presented by

Rhee (1968) will be introduced in the first part. The numerical algorithm developed in this work is discussed in the second part. The main focus in this part is on the detection of sharp fronts (discontinuous solutions) and component interactions that occur for multicomponent nonlinear equilibrium isotherms. Only systems without selectivity reversal are considered. Examples in the last part demonstrate the algorithm for extensions of the multicomponent Langmuir isotherm and a simple form of the Fritz-Schltinder isotherm.

Mathematical modelling

By neglecting mass transfer effects a system of hyperbolic partial differential equations is obtained. The mathematical model is based on analysis of the eigenvalues and eigenvectors of the gradient matrix of the multicomponent equilibrium isotherm. The solution in concentration space is obtained by combined integration along the eigenvectors and analysis of this concentration path with respect to its physical relevance. Its corresponding image in the physical plane diagram is represented by wave trajectories. The slope of these trajectories, which actually represent reciprocal concentration velocities, are found to be directly proportional to the eigenvalue of the gradient matrix of the equilibrium isotherm.

Fundamental equations A liquid isothermal system with a uniformly packed bed with constant void fraction e, constant flow rate, constant feed and initial conditions is considered. Equilibrium is established immediately at every location and instant. Any influence on adsorption caused by mass transfer effects in either radial or axial directions is neglected. Generic favorable and unfavorable multicomponent isotherms without selectivity reversal are considered. The individual components are classified by increasing respectively decreasing adsorption affinity. For gas systems an additional term including the socalled sorption effect which influences the flow rate due to volume change has to be considered. The fundamental mass balance on any component is obtained by considering a differential element of the column. Introduction of dimensionless time (t=u~Z) and distance (x=z/Z) as the independent variables and a phase factor (7/=(1 - ¢)/e) finally leads to the system of fundamental hyperbolic partial differential equations

Oc, (

-- +

l+'q

~n~)Oc,

=0,

(1)

where the change in total adsorbed amount of component i is represented by the directional derivative (~n~)l (~c,)

~n~ --= ~c~

Ont dc~ + . . . + -Onl - dcN __-_ ~ Oni dc~ Ocl dc~ OCN dci jffit Ocj dG"

Multicomponent nonlinear chromatography The components are arranged in such a way that l
nj/ct

(2)

where aj~ is defined as the relative adsorptivity or the separation factor. When the sequence of components defined by (2) is maintained throughout the concentration region considered this ensures that the order of adsorption and desorption remains unchanged during the movement of the individual components through the column. Consequently, component 1 is eluted first while component N is eluted last. (1) can be rewritten in the form of the well-known conservation laws

8 e(x,t)+ ~8 f(ct(x,t) ..... cN(x,t))=O 8x

(3)

where e is an N-dimensional vector of conserved quantities and f represents the flux function. More properly, the flux function can be interpreted as the column isotherm where the ith element is defined as f~=f~(cl.....c~)=c,+ ~ln,(cl..... cN).

(4)

Thus the column isotherm frepresents the total concentration in the differential element expressed per volume of fluid in the fixed bed.

Continuous solution In the following, solutions with piecewise constant initial and boundary conditions for the system of (1) are considered where pulse injections of duration ti into an initially and finally empty column are made c~=

initial condition: c°=0, feed condition: ci=cYi, final condition: c,'=0,

O<--x<-L,t=0 O<--t<-ti, x=0 (5) t>-t~, x=0.

This will define a Riemann problem. The iniection of a pulse into a column can therefore be assumed as two initially independent Riemann problems: an adsorption and a desorption step. The two will begin to interact after a while resulting in the desired separation of the components. Both steps taken together are often referred to as elution chromatography. It has been shown in the literature (Lax, 1957; Courant and Friedrichs, 1976) that for hyperbolic partial differential equations of this type the problem solution for each component can be expressed as a function of (xl t) i.e. the slope o"kof the trajectory of the corresponding concentration in the physical plane. On the other hand it is evident from the inverse function of ct that each concentration on the trajectory may also be expressed in terms of any other concentration on the same trajectory. Consequently the problem can be formulated in a new parametric representation, e.g., to

261

The independent parameter to then runs along the solution curve in the concentration domain or hodograph space while continuously changing from the initial state toi to the final state a/. (1) now can be rewritten as

(,o[ on].o),c -~x +

l+r/~c

~

~---~=0.

(7)

Since c~ is changing along the solution curve in the hodograph space one can immediately see that dc--2# 0 dto and consequently the first part of (7) must be zero giving (dt)

o'*= ~

8mlOx [

~n,)

,~=- o~/at = ~ 1 + ~ ¢

.

(8)

(8) which is also known as the kinematic wave equation is the key to the whole solution. Since (8) applies to all components the following condition must be met ~)nl ~c~

~n2 ~c2

~n N ~cN

tr k - 1 - -

= A k.

(9)

(9) requires the directional derivative of the equilibrium isotherm to be equal for all components at the actual state. The matrix representation of (9) de ( V n - Aq) ~ =0

(10)

with I representing the identity matrix, 0 the null matrix, e the column vector of concentrations finally reveals that Ak represents the kth eigenvalue of the gradient matrix Vn of the multicomponent equilibrium isotherm whereas (dc)l(dto k) corresponds to the kth eigenvector of the system. Since the system of PDE's is fully hyperbolic one obtains N distinct eigenvalues which are different from zero. Since the mathematical form of the equilibrium isotherm in general may be very complex it can easily be seen that no general analytic solution of the problem is available. Analytical solutions for (10) have only been found for a few special cases such as e.g., the multicomponent Langmuir isotherm. Classical numerical schemes lead to catastrophic results since unusual phenomena (shock waves, simple waves, wave interaction) occur depending on the initial and boundary conditions. In general, only special numerical procedures like e.g., high resolution methods (Harten, 1983) or modified orthogonal collocation schemes (Chiang and Hwang, 1989) can be applied. A different approach will be followed here. The gradient in the hodograph space with respect to cN can be calculated by the elements of the kth right eigenvector

dc__L = 8c/Oto

r~

dcN ac~ato r~ (6)

where r~k represents the ith element of the kth eigenvector.

262

R. ZENHAt:SERN and D. W. T. RIPPIN

The entire continuous solution will then consist of N pieces of concentration paths / " s which connect the initial and the feed state. They are in a descending order f r o m / ~ originating at the feed state to _F~ ending at the initial state which is in agreement with the adsorptivity order. The vertices at which the individual /-'s meet represent plateau zones. Furthermore it has been shown (Liu, 1975; Rhee et al., 1989) that there exists such a path and that it is unique. In parallel to the integration over the hodograph space one can readily construct the solution in the physical plane. To each point on t h e / ~ there exists a corresponding trajectory in the physical plane. For simplicity, however, a discrete number of trajectories will provide sufficient accuracy. The image of t h e / ~ in the physical plane will therefore be represented by a bundle of m straight lines forming the kth simple wave with the given slope o'k(j)= 1 + r/A k~/)(j= ! ..... m), each beginning at the origin. Mathematical solutions with multiple concentrations at the same time and location are not physically realistic. Consequently, the following condition must be fulfilled

= ~ - >0,

(11)

which simply means that the trajectories in the physical plane must fan clockwise when proceeding from the time axis to the distance axis. Equivalent to condition (11) is the relation At(j- I)< At(j)

(12)

when j is numbered in the direction of time. Smaller eigenvalues of the gradient matrix therefore correspond to larger concentration velocities Uc. One can further observe, that the entire solution can be constructed by means of the gradient matrix and the initial and feed conditions alone.

trY=

At) Ani Azz = l + r / A c ~ '

(13)

The A-terms refer to the difference between the left and right side of the discontinuity where the left side is closer to the feed state and the right side closer to the initial state. As with the kinematic wave equation for the continuous case (8) one can readily see that (13) also applies to all components. One therefore obtains the compatibility relation An I _ An2 . . . . . ~n~, Acz Ac~ '

Ac~

(14)

which is a set of ( N - 1) algebraic equations. (13) and (14) are also known as the R a n k i n e - H u g o n i o t relationships in the theory of conservation laws. The discontinuous solution then can be obtained by knowing the state on one side of the shock and either the reciprocal velocity of the shock wave o "~ or one concentration on the other side of the shock. A discontinuous change in concentration of one component has a direct effect on other components through the multicomponent equilibrium relationship. Consequently discontinuities occur at the same time for all components. The image of the shock path in the hodograph space can be represented by a ~k which is similar to the continuous case. For infinitely small differences at a point the shock path ~k becomes tangential to the continuous path/~. An ambiguity still remains concerning the direction of jump. The solution to (13) and (14) is totally symmetric with respect to the correct sign. Therefore, an additional condition has to be formulated to find the reasonable physical solution. The extended entropy condition (Liu, 1975: Rhee et al., 1989) helps find the correct physical solution. Extended Entropy Condition: Let ~.o k be the curve .vk emanating f r o m the point c t in the hodograph space and consider another point c r on .~o k. I f the two states are to be connected by a k-discontinuity with c t on the left-hand side c r on the right, we must have

Discontinuous solutions

One could also imagine a case where the characteristics would fan counterclockwise or partially clockwise and partially counterclockwise in the physical plane when proceeding from the time axis to the distance axis. Due to the different velocity of the characteristics the shape of the peak becomes asymmetric with time. After some time the characteristics begin to cross. Since each characteristic is associated with a composition one would then obtain overhanging concentration profiles having more than one composition at the same location. This has of course no physical meaning. To find the physically realistic solution one has to introduce a discontinuity which is referred to as a shock wave. Across the shock wave a mass balance is formulated by applying the same assumptions as in the continuous case. This finally leads to an alternative form of the kinematical wave equation which is also known as the jump condition

o'~(c ~, e ' ) - > ~ ( c ~, c)

(15)

f o r every c on ~o k between e t and c r. It has been shown for continuous solutions that the eigenvalue of the gradient matrix of the equilibrium isotherm is continuously decreasing when proceeding from the feed state to the initial state. The occurrence of discontinuities must therefore be related to the reverse behavior of the eigenvalue, i.e. an increasing eigenvalue along/*. The kth eigenvalue is monotonically increasing when following the whole concentration path /~ between two succeeding states.

dA t d---ff>0.

Consequently, the corresponding characteristics in the physical plane would fan counterclockwise. For this case a shock path ,vt has to be constructed starting at the feed state with the concentration c I. Since the discontinuous

Multicomponent nonlinear chromatography path £k proceeds in the same general direction as the continuous p a t h / ~ the eigenvalue is also monotonically increasing and thus the extended entropy condition is satisfiedwith c I f c t and c t f c : Both statesare therefore

263

Constructing the solution

The computational work basically consists of the integrations over the N state transitions. The algorithm for the kth integration is shown schematically in Fig. 2. It consists of the following steps:

connected by a k-shock where the shock slopefulfillsthe shock conditionas listedin Table I. On the other hand one may also have cases with partially continuous solutions during the transition 1. Integration over the transition k, between two states.Such cases are summarized in Table 2. Analysis of the concentration path and the characterI and Fig. I. istic slopes, Table 1. Conditions for different shock panerns

Pauem

Shock condition

Reference

Shocks Left semishock Right semishock Contact discontinuity

O'/(C~>O'k'(Cl,C~>O'k~(C~ o'/(¢b=o'k'(c~,c~>o'/(e~)

Fig. l(a) Fig. l(b)

O'/(Cl)>Crk'(Cl,C~=Cr/(C~ Fig.l(c) cvkC(c~ffio'tJ(c~cJ+t)fo'/~(e j+i) Fig.l(d)

c~

t sk

(b) (a)

--~il°a(a)

ts

•

(a)

X

t

t.

:Y sk

V

C

(a)

ct

(b) - - ~

(a)

• x

ct (b) (e) ts

(a)

X

ci

(d)

t. (a) x

X

Fig. 1, Summary of possible shock patterns: (a) shock waves, (b) left-side semishock, (c) right-side semi-shock, (d) two-sided contact discontinuities.

R. ZENH.~USERNand D. W. T. RIPPIN

264 3. Solution of wave interactions.

Each of these topics will now be discussed. Integration over the transition The integration has to be performed in the direction of decrease of the component k since no initial guess of the coherent path can be made at zero concentration of component k. For the integration over the transition k the component k is chosen preferably as the independent parameter since its values are known on both sides of the transition, namely ckt=O for the lower boundary and ck"= c J respectively cku= cki for the upper boundary. The integration then can be formulated as follows:

r4,

c= / Fdck Jq where the N - 1 elements of the integration function F' are obtained by the elements of the kth eigenvector r k of the gradient matrix of the equilibrium isotherms F'(i)= dc-2= rk(i) . dck rk(k) Consequently, one has to find the set of eigenvalues A and eigenvectors r of the gradient matrix along the

integration path. However, only the kth eigenvalue and eigenvectors will be used since the integration is performed in the kth direction of the hodograph space. With the calculation of the coherent concentration path one immediately finds the slopes of associated characteristics in the physical plane by inserting the kth eigenvalue into the kinematic wave equation o~=I+v/A k at each integration point. Since the individual characteristics in the physical plane are needed for subsequent calculation of chromatographic profiles or wave interactions a distinct number of integration points has to be chosen. Consequently the integration can be performed by any integrator which calculates the solution at a given number of integration points. The integration is repeated N times, i.e. for each component until the final state is reached. Analysis of the coherent concentration path The coherent concentration path P then has to be analyzed in order to find the physically correct solution. This can be solved by two simple methods • Analysis of the sequence of the slopes ~'s, • Find the coherent concentration path in the n-c diagram.

creed ' k = ko [

+

-

Integration over I ck untilto ck = o

+

l

Concentration path F k with eigenvalues ~k l i

increasing @

~

decreasing

I increasing and decreasing [Shock wave I

I Mixed shocks

+

I

I Simple wave[

ewk

- onstruct solution in physical plane Obtain profiles from physical plane

I

Fig. 2. Basic tasks for the solution of single step inputs.

I

Multicomponent nonlinear chromatography

265

The simplest method is only to analyze the sequence of the slopes of the characteristics in the physical plane. From the requirement that no more than one concentration of the same component can exist at the same time and location, one can conclude that the characteristics of the transitions k must fulfil the following condition to be physically reasonable and consequently form a simple wave C k. Rule 1: Existence of simple waves The slope of the characteristics in the physical plane must fan clockwise when proceeding from the time axis to the distance axis. If Rule 1 is fulfilled then the coherent concentration path obtained represents the physically reasonable solution. If Rule 1 is not satisfied then Rule 2 will check for the presence of a shock wave. Rule 2: Existence of shock waves The slope of the characteristics in the physical plane must fan counterclockwise when proceeding from the time axis to the distance axis.

1. Perform same integration over the wave transition of interest as in the fully continuous case. 2. Choose a component, e.g., k at the kth wave transition and find the stationary phase concentration nk through the isothermal relationship. 3. Construct the envelope between the start and end point of the integration by satisfying the condition that the slope of the envelope must either increase or remain constant. 4. Portions of the envelope which coincide with the previously calculated coherent concentration path of the integration are then represented by continuous changes whereas portions of the envelope not following the path represent discontinuities between the start and end point of this part of the path. 5. Construct the corresponding wave pattern in the physical plane from the above analysis. Continuous portions are represented by simple waves, discontinuous portions by shock waves.

Both rules have also been formulated by Glueckanf (1949) and Basmadjian and Coroyannakis (1987). For the multicomponent Langmuir isotherm Rule 1 and 2 prove to be enough since only shock waves and simple waves are obtained. For other isotherms, however, more sophisticated checking needs to be done in order to find the correct coherent wave path with its corresponding wave pattern in the physical plane. To determine continuous and discontinuous portions in the solution the complete coherent wave path of the actual transition has to be analyzed. From a physical point of view one can readily agree that the concentration velocity of a coherent composition within a wave transition must either decrease or remain constant during its movement through the column, i.e., the shape of the actual wave cannot become overhanging with time or distance. On the other hand the velocity is obtained by the reciprocal value of the slope of the trajectory in the physical plane diagram

The procedure is shown in Fig. 3 for an artificial isotherm. As soon as the slope of the integrated path in the n-c diagram does not fulfil the criteria stated in paragraph 3 of Rule 3, a shock must be introduced which continues until Rule 3 is satisfied again. From the resulting final coherent path one then can immediately derive the image in the physical plane diagram. The combination of these N integrations will now complete the solution of either of the step inputs in the physical plane diagram. Since the slope and composition of each characteristic is known one can easily construct the corresponding chromatographic profiles.

~nk where the directional derivative must be equal for each component. This is required by the coherence condition (9) formulated earlier. Consequently the directional derivative must increase or remain constant along the coherent concentration path. This can be verified by analyzing the coherent path of the mobile and stationary phase concentrations in the two-dimensional n--c diagram for any of the components. By definition the directional derivative evaluates the slope in the n--c diagram in the direction of integration which in this case results in the coherent concentration path. Projecting this path into the two-dimensional plane of e.g. component k, one can verify the slope of the coherent curve between the two states which are connected by the actual wave transition. The following rule can therefore be formulated Rule 3: Determination of the physically correct coherent concentration path in the n--c diagram:

Solution of wave interaction

In the modelling of elution chromatography the rectangular input pulse of a sample can be represented by two initially independent Riemann problems. The first Riemann problem describes the loading of the column to the feed concentration while the second Riemann problem represents the elution of the loaded column. Each of these partial problems may involve both shock waves and simple waves or a combination thereof. Each wave is travelling at a different speed through the column. To combine the two Riemann problems one has to consider wave interactions between the different components. At first sight the calculation of such wave interactions seems to be rather difficult. In the previous treatment it has been shown that the coherent concentration path in the hodograph space is obtained by the successive integration along the eigenvectors of each component. By simple coordinate transformation of the hodograph space into the eigenvector space ~bi=Rrci, where R contains the N eigenvectors, the concentration path becomes orthogonal, i.e., the kth wave only affects the kth concentration along the kth eigenvector in the eigenvector space. However, no analytical solution is available for the transformation into the eigenvectors and consequently the transformation has to be applied to

266

R. ZENH,~USERNand D. W. T. RIPPIN

,,F

nk

a

f,i

d:," J/j t

/ ~ /

(b) f

-- -- shocks

/~h

,

'Ck (c)

(a)

Fig. 3. Determination of the coherent concentration path and the physical plane patterns in the n--c diagram (a). Corresponding image in the physical plane for loading (b) and elution (c) of component k. a few selected points along the coherent concentration path. On the other hand the same may also be true for the waves of another Riemann problem. Since in our case both have the same multicomponent equilibrium isotherm they will also share the same eigenvector space. Consequently, the wave interactions can be calculated in terms of the orthogonal intersection points of the corresponding waves in the eigenvector space. The backtransformation into the hodograph space is possible if the transformation matrix consisting of the eigenvectors is available. Since the eigenvectors remain constant for the multicomponent Langmuir isotherm the backtransformation is straightforward for this type of isotherm. For the case where the backtransformation matrix from the eigenvector space into the hodograph space may be not available one would have to perform an integration along each trajectory in the direction of the corresponding eigenvectors to construct the grid of coherent paths in the hodograph from which the wave pattern in the physical plane can be deduced. Simple wave transmission by a shock wave

In Fig. 4(a) the simple wave C mof the mth component is overtaking the shock wave S k of the kth component. The simple wave C" in this case is represented by five characteristics. At point A the shock wave S k starts to

t

interact with the simple wave C m. During the transmission the velocity of the shock wave changes continuously until point E is reached. The characteristics of the simple wave, however, experience a discontinuous change in velocity. This can be observed by the discontinuous change in velocity over the shock wave resulting in the characteristics l', 2', 3', 4', 5'. The corresponding image in the eigenvector space is shown in Fig. 4(b). Discontinuous changes are indicated as dashed lines whereas continuous changes are drawn in solid lines. The construction of the solution consists of three tasks • Determine the intersection points A, B, C, D, E, • Determine concentrations on the right side of the shock, • Find the slopes of the characteristics after the transmission by the shock wave. The first task is relatively simple. Since all characteristics i through 5 and both sides of the shock S k are known initially one can easily construct the first intersection point A. For the case where c(k)=0 on the right side of the shock the slope of the shock wave that corresponds to each of the remaining characteristics can immediately be obtained. Thus the remaining inter-

U s ' E

4'

1'

2'

3'

4'

I

I

I I I

I I I

5'

2'

' T

i

1 x

(a)

sk

i

3

4

em (b)

Fig. 4. Transmission of a simple wave by a shock wave in the physical plane (a) and eigenvector space (b) representation.

Multicomponent nonlinear chromatography section points are readily available and since sufficient conditions are met one can also construct the full composition on the fight side which corresponds to each characteristic. A distinct number of characteristics is chosen initially to construct the simple wave. Consequently, one can assign a slope and a set of concentrations to each characteristic in the simple wave. This also applies to the shock waves. The procedure that can be applied to each of the wave interactions shown above can be outlined as follows: 1. Ensure that all the states are known that precede the zone of interaction. 2. Calculate the point of intersection of the straight shock lines and simple wave trajectories. For the part of the shock wave that transmits or absorbs the simple wave the term "straight" applies to the piece between the succeeding simple wave characteristics. However, the average of the shock slope corresponding to the concentrations on the current and next trajectory of the simple wave should be taken. In fact, this corresponds to an integration of the shock path. More accurate results may therefore be obtained by using a higher order integrator. But this would require the calculation of additional internal points. 3. Determine the set of concentrations on the right side by taking advantage of the already known states and the properties of the shock wave. This will be needed in order to calculate the new slope of the simple wave characteristic or part of the shock wave. 4. Calculate the new slopes of simple or shock waves by using the kinematical wave equations that apply to simple waves and shock waves. .,

c

267

Alternatively the problem can also be solved by first constructing the solution in the eigenvector space, backtransforming into the hodograph space, calculating the corresponding slopes of the shock wave and characteristic of the simple wave and finally obtaining the intersection points. This procedure has already been pointed out.

Simple wave absorption by a shock wave Simple wave absorption always occurs between a simple wave and a shock wave of the same component. The graphical representation of this phenomenon is shown in Fig. 5(a). The corresponding image in the eigenvector space in Fig. 5(b) also illustrates that this behavior only affects the kth eigenvector axis in the eigenvector space. At the beginning of the interference of the two types of waves the constant plateau zone F disappears. Furthermore the magnitude of the chromatographic peak is continuously changing. Each characteristic that intersects the shock wave is immediately absorbed while changing the path of the shock wave. The computational task is the same as already described except that each characteristic is absorbed at the point of intersection. Simulation e x a m p l e s

The multicomponent Langmuir isotherm The success of the multicomponent Langmuir equilibrium isotherm can be explained by its practical and mathematical simplicity, namely

N,K#i hi=

.

One can immediately draw several fundamental conclusions: • The equilibrium parameters are associated with pure component favourable behaviour. For the determination of these parameters one only needs to find the parameters of the individual components. No extensive measurements are needed to obtain the set of 2N parameters. • Calculating the selectivity order,

~nm

Simple wave C :O'm=1 + rt~c, "

Shock wave Sk: o-~= 1 + r/Ank

ACk

5. Continue with step 2 until the interaction disappears or the desired range of interest has been covered.

I

F

C

E

L _ _ L Sk 1 2 3 4 5

/ / / ~

(I) (a)

(16)

j=l

x ~v (b)

Fig. 5. Absorption of a simple wave by a shock wave in the physical plane (a) and eigenvector space (b) representation.

268

R. ZENH~USERNand D. W. T. P.JppIN Table2. Parametersfor the simulationof the 4-componentLangmuir example Component

1

Saturation valueN~ Equilibriumvalue K~ Concentrationc~fe~ Columndiameterd Void fraction ~ Fluid velocityu Cross sectionalareaA Sample volume V

mol/L3 0.5 L3/mol 1 mollL 3 0.05 L

2

3

4

0.5

0.5

0.5

2

3

4

0 . 0 5 0 . 0 5 0.05 2 0.4

LIT

1

L2

1 1.25

L3

NzKz N3K3 <...< NNK~ 1 < NiKI < ~ NtKi ,

one finds that it does not depend on the actual composition but only on the equilibrium parameters of the isotherm. Therefore the selectivity order remains constant for the whole separation process. Since one can express the selectivity order by constant values regardless of the actual composition, considerable simplifications are possible. For such systems analytical solutions for multicomponent systems have been found in the past (Helfferich and Klein, 1970; Rhee et al., 1970). This not only includes the computation of the coherent path but also the whole interaction analysis in the hodograph and physical plane representation. Nevertheless the application of such analytical expressions to systems with more than two components can become very cumbersome. This includes the determination of the eigenvalues and the eigenvectors of systems of higher orders and especially the calculation of wave interactions. The algorithm in this paper, however, provides an excellent and rapid means for the calculation of multicomponent Langmuir equilibrium systems. For

the subsequent simulations the parameters chosen are summarized in Table 2. The dimensions herein are only specified by generic units of time T and length L. In practice, however, appropriate dimensions may be assigned. The chosen values of the sample volume V, cross sectional area A, and fluid velocity will result in an initial pulse length of t~= 1 on the time axis. By combination of the loading and elution process one is able to compute the separation of multicomponent systems by chromatography. This step can be done by constructing the grid of simple waves and shocks waves in the physical plane diagram (see Fig. 6). One preferably starts with computation of the integration of the shock path S N with the simple waves C~-C N-I (transmission) and C N(absorption). At each intersection point the state on the right side of the shock is immediately available since the concentration cN=0 on the right. This process is repeated for each shock wave. The solution for a 4-component mixture is illustrated in Fig. 6. The following observations can be made: • The pulse injection is indicated by the interval between 0 and 1 on the time axis. • The solution is further characterised by four shock

14 12 10 8

6 4

2 C 0

0.5

1

1.5

2

2.5

3

3.5

4

4.3

5

Distance x Fig. 6. Physical plane diagram for the separation of a four componentmixture using the multicomponentLangmuir isotherm: Simple waves (--), shock waves (- - -).

Multicomponent nonlinear chromatography

•

•

•

•

waves and four simple waves which start at the beginning and end of the pulse injection respectively. Since the slopes of the waves represent inverse concentration velocities one can easily observe that the simple waves will tend to overtake the slower shock waves. Two types of wave interactions can be observed: simple wave transmission between the shock wave of component k and the simple wave m and simple wave absorption between the shock and simple wave of component k. Each wave type will become slower during wave transmission. This is indicated by steeper slopes after the interaction. Two points in the physical plane deserve special attention • 1. Point (0.7, 1.9) indicates the location where wave interaction begins. Before this point the concentration profiles in either direction can be constructed by the solutions of loading and elution. • 2. At point (3.2, 10.6) the transmission of shock wave S 4 through simple wave C 3 is complete. Although other components may already be separated earlier, this represents the minimal theoretical column length required in order to completely separate all the componcnts. This point will subsequently be called the characteristic separation point.

• The intermediate plateau zones will disappear as time and distance evolve. Finally there will only remain a shock front immediately followed by the corresponding simple wave for each component. The effect of variation of the equilibrium parameters on the location of the characteristic separation point is summarized in Table 3. The results are consistent with physical intuition. Increasing the equilibrium constants K~therefore causes stronger sorption than in the standard case and consequently decreases the required column length. On the other hand, longer columns will be needed to achieve the same degree of separation if smaller values for K~ are considered. One could also expect variation of Ni, which actually represents the "number" of adsorption sites to have a direct influence on the required column length. For the cases studied a directly proportional influence has been found. By cutting the physical plane diagram in the direction of distance at several times, column profiles are obtained. This is illustrated in Fig. 7 at four selected

269

times. At t= 1 the pulse injection has just finished. The wave pattern at this point consists of four shock waves while the simple waves have not yet started to form. At succeeding times t=4 and t=6 ongoing separation of the components can be observed while at t= 10 the mixture has completely separated.

Polynomial extension to the Langmuir isotherm The multicomponent Langrnuir isotherm can be regarded as the ratio of two polynomials of the first order. Adding product terms of different concentrations simply represents the continuation of the polynomial series in the numerator and denominator of the fraction giving

ni =

+ d f , cj u I+ E bFj+ Eje,c,cj j=l

Distance ratio Time ratio

0.75 2.6

(17)

"

The product terms herein represent additional contribution of interactions between different components to the stationary phase concentration. For dj=0 and ej=0 (17) reduces to the classical multicomponent Langmuir isotherm. Various cases have been studied for this type of isotherm. The equilibrium parameters used for the simulation runs are summarized in Table 4. The parameters for ai and b i are the same as chosen in the simulation of the pure multicomponent Langmuir case (Table 2). For case (a) the resulting physical plane diagram is illustrated in Fig. 9. As could be expected the trajectories are much flatter than in the pure multicomponent Langmuir case. Due to the additional interaction represented by the product terms of concentration the efficiency of the separation has decreased. The range of concentration velocities has become smaller. This can be observed by the simple and shock waves which move closer together and more nearly parallel. As a result, the required minimal column length and minimal time to achieve complete separation is increased. The various cases are compared in terms of the characteristic separation points in Table 4. While nearly double of the column length is needed only a modest increase in time is required to achieve the complete separation of component 4. An interesting behavior in case (a) can be recognized from the breakthrough histories in Fig. 8. The elution of the four component mixture beginning from the feed state first leads to a modest increase in concentrations while component 1 is eluted from the column. This may

Table 3. Comparison of distance and time ratio of the characteristic separation point with reference example (Fig. 6) through variation of equilibrium isotherm parameters

10XK,

"

0.1XK, 2XN, 2.6 1.2

0.5 0.9

0.5xN~ 2.0 1.3

R. Z~'~H~USEa~and D. W. T. RIPPIN

270

(a), simple wave C 1 in (b) and simple wave C 3 in (c) of each figure. The results of these figures can be summarized as follows:

be explained by the higher influence of dr in the equilibrium relationship of components 2-3 which has a slight displacement effect on component 1. The rear front of component 1 is consequently very steep compared to the rear boundaries of the other components. More information on the development of the composition with time is obtained by looking at the hodograph diagram. In Fig. 10 the hodographs of the undisturbed simple waves (i.e. before any of the simple waves meets the shock S 4) are illustrated, namely simple wave C ~ in

•

The nonlinear shape of the curves increases as contributions due to dj and e1 are increasing (see Fig.

iO). • A strong influence on the nonlinear shape of the concentration path is caused by the parameter d~

(case (d)). • The parameter e/alone seems not to have a great

0.08 0.0'7 0.0~

t=l

0.0~

c4

--

0.04 0.03 0.02 0+01 i

2

3

4

s

0.08 0.07

~

0.0($

t=4

0.05

c4

0.04 0.03

/ /

f-!

!/

0.02

__L

0.01

2

3

4

5

O.OI 0.07 0.06

i

0.05

-! ~]

t=O

, / ....... 0.04 0.03

,../,! I/

0.02

o

I I

i/i I/ i l I

0.01

~

4

S

6

5

0

0.08 0.0"/

/I /'

0.0~ 0.05

c4

0.04 0.03 0,02

0.01 01

/ ........]/ /

i

,/ !,/, Distance x

Fig. 7, Column profiles for the multicomponenlLangmuirisotherm at different times. Components: I (--), 2 (- - -), 3 (...), 4 (...).

Multicomponent nonlinear chromatography

271

Table 4. Equilibrium parameters and characteristic point for the case study of the four component extended Langmuir system Case

Components

ai bi

a b c d e

dj ej ~ ej dj ej ~ ej dj ej

1

2

3

4

0.5 1 0.2 0.5 0.5 0.5 1 1 1 0 0 1

1 2 0.4 0.5 0.5 0.5 1 1 1 0 0 1

1.5 3 0.6 0.5 0.5 0.5 1 1 1 0 0 1

2 4 0.8 0.5 0.5 0.5 1 1 1 0 0 1

impact on the form of the hodograph representation (case (e)). This can be observed from the linear shape of the concentration paths in Fig. 10. Since the classical Langmuir relationship also leads to straight lines it is difficult to distinguish the behaviour from that of a classical multicomponent Langmuir system. These qualitative observations may provide useful insights for the determination of an appropriate equilibrium isotherm of a given system. Furthermore it could also be shown that linear curves in the hodograph space must not necessarily be related uniquely to Langmuir isotherms. The range of application of the polynomial extensions may therefore be considerably increased since Langmuir behaviour may be partly satisfied, i.e., over a limited concentration range, while it may provide a useful representation in other regions where the classical multicomponent Langmuir isotherm seems to be inappropriate.

Exponential extension of the Langmuir isotherm A convenient way to extend the classical multicomponent Langmuir isotherm is to add exponents to the concentrations. A simplified form of the Fritz-Schltinder isotherm has recently been proposed by Basmadjian et al. (1987a) giving a~oci ni= di+ .l~,'=aoc~'J"

Characteristic point

(18)

Only exponents in the denominator are considered. Since (18) satisfies the Gibbs-Duhem requirement (Glueckauf, 1949)

Xchar /char

3.4

11.5

3.2

10.9

3.3

11.3

3.3

11.2

3.2

10.7

and Henry's Law at low coverages, basic thermodynamical requirements are met. Multicomponent equilibrium isotherms of the form in (18) are able to represent multicomponent systems with favorable equilibrium isotherms (Type I of the BET classification). However, it can also be extended with modest changes to Type I, II and III systems (Basmadjian et al., 1987b). For this illustration the following form of (18) will be used

ni= 1+

a~£i N j= I

•

(19)

Three different cases will be presented: (a) Increasing values of b~ where 0
272

R. ZENH#kUSERNand D. W. T. RIPPIN

waves. This result is also shown by means of the characteristic separation point in Table 5, where this

effect is expressed mainly by the considerable differences in time.

O.1 0.09 0,08

L=0.5

0.07 0.06

c~

0,05 0.04. 0.03 0.02

i!i L~ i

0.01 0

0

~

2

~

1~

1'2

lg

14

l'S

20

18

20

18

20

0.1 0.09 0.08

L=l.5

0.07 0.06

c/

0.05 0.040.03 0.02 0,01 0

I

2

ili

4

,

6

8

10

12

14

16

0.1 0.09 0.08

L=2.0

0.07 0.06

c/

i i~.,-.~

0.05 0.0.4 0.03 0.02 0.01 0

2

4

6

8

10

12

14

16

0.1 0.09 0.08

L = 4.0

0.07 0.06

ci

0.0~; 0.04 0.03 0.02 0.01 0

0

2

4

6

V1 8

10

i., i \.

i

!....... "L

i

.\ 12

14

16

20

Time t Fig. 8. Breakthrough history profiles for the extended multicomponent Langmuir isotherm for case (a) at different column lengths. Components: I (--), 2 (- - -), 3 (.. 9, 4 (-. -).

Multicomponent nonlinear chromatography

273

14 12 10 8 6 4 2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Distance x Fig. 9. Physical plane representations of case (a). Simple waves (--), shock waves (- - -). The Fritz-Schliinder isotherm A popular form where exponents are added to the numerator and denominator has been proposed by Fritz and Schliinder (1974, 1981a,b) nl =

a/oc,%

•

(20)

di+ j =~l aoc~'J

Considerable drawbacks are associated with multicomponent systems described by (20). Although it may well represent the equilibrium situation over a given range of concentration 1 serious problems arise especially when moving towards low concentrations. Since the velocity of a coherent composition is proportional to the inverse of the directional derivative (5~ni)

-J

Uci OC ~C i

one can readily see that for values of the exponents b~0 and b o less than 1 an undefined situation may develop if the concentrations move towards zero. The concentration velocity uc, moves towards zero for this case resulting in infinite long concentration peaks in the chromatograms. Since NX(2N+3) parameters can be specified it is adaptable to a wide range of multicomponent equilibria. For low concentrations, however, it must be replaced by a thermodynamically consistent isotherm. It will be demonstrated here that the numerical algorithm proposed can also handle these complicated forms of nonlinear isotherms. The simulations are based

on work by Balzli et al. (Balzli, 1977; Balzli et al., 1978). For the subsequent simulation a ternary mixture of Phenyl/t-Amylalcohol/2-Butanol is considered. A pulse of 6.4 t duration and a concentration of 1000 mg/1 for each component is simulated. The measurements by Baizli have led to the FritzSchliinder parameters shown in Table 6. One can readily see that the exponential parameters b are less than zero. This represents a major problem and can be demonstrated by means of component 3 whose equation in the absence of the other components becomes 5~

0.22

n~= ~ ,L C 3

~ =~z: c~022 .

(21)

1 c3

Calculating the derivative with respect to c3 which is needed to obtain the slope or the simple wave characteristics one finds that the stationary phase concentration tends to infinity as its corresponding mobile phase concentration becomes zero. For the modelling of such isotherms this behaviour has to be treated separately. For the current simulation this case has been handled by applying the concept of coherence. By this means the following wave types could be found: 1. Loading step: - - Introduction of component 1: shock wave, - - Introduction of component 2: left-side semishock, - - Introduction of component 3: shock wave. 2. Elution step: m Elution of component 1: right-side semi-shock, Elution of component 2: right-side semi-shock, Elution of component 3: simple wave.

--

274

R. ZE~;~USERN and D. W. T. RIPPIN

The left-side semi-shock of the loading step, however, must be a consequence of the inconsistency of the isotherm. Regardless of how many integration points are chosen the last trajectory always introduces this artificial effect just at the point where c2 becomes zero. On the other hand the beginning of the elution of components 2 and 3 is always started with a step decrease. This can also be observed in the breakthrough histories in Fig. 15. The integration has been carded out with forty wave trajectories. The physical plane diagram (Fig. 14) shows

a significant difference between the components. The adsorptivity of Component 3 is considerably higher requiring high desorption times. On the other hand the low adsorptivity of components 1 and 2 will also demand higher column lengths. Discussion

Simulations of pulse chromatography for different multicomponent nonlinear isotherms have been presented.For the range of concentrationsand equilibrium

0.053 .,~÷÷÷++ ..... . . . . . . ~ ÷ ÷ + 0.052

.........

~ . ~ ÷ +

0.051 ..............................

0.05

C1

÷,-~.+

............

ooo ° 000

0.049

O0

00

0 0 0 00

0.048

O0

0

0

0

0.04"7 0

0 0

0.046 0.045

O 0

0 0 0 0

0

0

0 0

0

C1

0 0 0

o.~os

0

0.51 0.6,5 0.~

i

0.625 0.53 0.635

0.04

0.(}45

0.05

C4 0.054

iiiii?/.........

0.052 0.05 . .+ .°~ "~ ÷ * +' ÷ "1"

0.048

.--

-"

0.046

C1

0.044

..o-

.--"

..-"

0.042

° ooo o

ooo

oooo

.--" 0.04 0.038

C 2

0.03~

0.034

o.51

0

0.52

0.53

o.~

0.55

0.56

0.07

C3 0.05

O.045

÷ ~..~.°÷+÷

+.~.~""

0.04 .~o,+--o~'°'°

C1

0.035

0.03

~

o

+ "~ "+"" _ - ~ o

o

.

. - -

.~.O,,¢.:~"'~" ....ooo°°oq ÷.::° * ....oooo ~.:~.~+" ..--" oOO o ..~...~" . o . o - " oOOo° .~..b..'*" ' t . - - " o o o o o °

z

z°°

.~:.--~ 0.025

0.02

Ca 0.1500

i

0.01

i

0.015 0.~

i

i

0.025 o.b3 0.035 0.~,, 0.~45 0.05

Fig. 10. Hodograph representation (a) Complete elution of component 1 (C~). (b) Complete elution of component 2 (C2). (c) Complete elution of component 3 (Ca). Cases: (a) (--), (b) (- - -), (c) (-. -), (d) (+), • (0).

Multicomponent nonlinear chromatography

275

Table 5. Equilibriumparameters for the ease study of the four componentextendedLangmuirsystem Case

Components

a~ ~ngmuir a b c

b~ ~ ~ ~

Characteristic point

1

2

3

0.5 1 l.O 0.2 1.2 0.5

1 2 1.0 0.4 1.4 1.7

1.5 3 1.0 0.6 1.6 1.3

parameters it has been found that the loading step is represented by shock waves whereas the elution step is described by simple waves. Differences, however, are observed in the hodograph representations of the simple waves. Simulations with the polynomial extension have shown only minor differences with respect to the physical plane representation. The characteristic separation point does not change significantly. The short increase at the start of simple wave C ~ seems to be particularly coupled with the parameter d,. which serves as weight in the numerator whereas parameter e~ supports classical Langmuir behavior in the hodograph representations. An interesting effect has been observed in case (b) where each of the two different effects becomes dominant during the elution of component 1 leading to the parabolic behavior shown in Fig. 10. Considerable differences from the classical multicomponent Langmuir isotherm have also been found for its exponential extension in the hodograph space, namely • bj< 1: Convex towards the concentration of cN,

4

Xch,~

tc~

3.2 14.0 1.6 2.7

10.6 43.8 6.2 9.0

2 4

1.0 0.8 1.8 0.9

• bi= 1: Linear towards the concentration of cN,

• bi> 1: Concave towards the concentration of cN. For b~<0 simple waves tend to be broader. A significant change of the slopes of the last trajectories within the eluting simple waves is observed (see Fig. 11). This can only be explained by a stiff behaviour of the multicomponent equilibrium isotherm along the corresponding coherent path. On the other hand for bj> 1 the simple wave is found to be significantly more dense leading to quicker lower values for the characteristic separation point. Combining case (a) and case (b) the same behavior can be found for the individual components. All simulations have been carried out using different numbers of trajectories to represent the simple wave. At these points the corresponding coherent composition has to be calculated in order to find the slope of the associated trajectory. This will be needed further for the calculation of wave interaction in the physical plane domain. The integration over the simple wave has been made by a conventional 4th order Runge--Kutta algorithm. The calculation of the shock path through the simple wave, i.e., the transmission or adsorption, in fact

14 12 10 8 6 4 2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Distance x Fig. 1 !. Physical plane diagram of case (a) of the exponentialextended Langmuirisotherm: Simple waves (--), shock waves (---). CACE22-1/2-J

276

R. ZEN.£USEP~ and D. W. T. RIPPIN 14 12 10 8

6 4 2 0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Distance x Fig. 12. Physical plane diagram of case (b) of the exponential extended Langmuir isotherm: Simple waves (--), shock waves (---). corresponds to another integration. For the present case the intersection points of the shock wave with the characteristics of the simple waves have been determined by taking the average shock slope corresponding to both characteristics. To obtain an estimation for the error made, comparisons of the coordinates of the time of the characteristic separation point have been made as a function of the number of characteristics used to represent the simple waves. Similar behaviour is obtained for the location of the characteristic point. The results are summarized in Fig. 16. As can bee seen at least ten characteristics will be needed in order to have an error in the region of 5 percent. The worst case herein is seen to be case (a) of the exponential extension. This is the effect discussed above where the slopes of the characteristics are asymmetrically distributed across the simple wave. A considerable number of trajectories will be needed to achieve acceptable error limits. For this case 60 trajectories have been chosen giving an error of less than 5 percent. The remaining simulations have been carded out with 40 trajectories on a VAX 6620. The computation time varies linearly between 5 and 20 minutes for 20 respectively 100 trajectories. For the error it has been found that error ocnpts J.2~ represented by the lines drawn in Fig. 16 for the corresponding cases. Conclusions

A theoretical study has been presented modelling multicomponent nonlinear equilibrium chromatogrnphy. Since the influence of nonlinear multicomponent equilibrium isotherms was of primary interest in this paper local equilibrium theory has been studied. The neglect of

additional mass transfer due to diffusion means that the solution provided by local equilibrium theory represents the theoretical optimum which could be achieved if perfect ideal conditions are established for the separation. In reality, however, mass transfer effects could not be eliminated completely from the separation process. This effect can be observed in chromatographic profiles which may be considerably smoothed compared to the ideal case. Therefore local equilibrium chromatography cannot exactly represent the real situation. Nevertheless local equilibrium theory has advantages which can favour its use in the design and development stage of a chromatographic separation: • The ideal situation shows the essentials of the chromatographic separation, the approximate peak shape, interactions between components etc., which are points of central interest in the design stage. • Fewer parameters are needed for the simulations compared to models claiming to represent the "real" situation by including mass transfer effects. • Evaluation of other multicomponent equilibrium isotherms can easily be carried out since the algorithm could be formulated in terms of generic stationary phase concentrations and gradient matrices which are the only functions to be specified for the new isotherm. • The transparent algorithm presented can easily be implemented. • The algorithm presented is found to be sufficiently fast for the types of multicomponent nonlinear equilibrium isotherms studied. Consequently, this provides a rapid and effective mean for the evaluation of different multicomponent equilibrium isotherms. The simulations have been carried out for three different favorable multicomponent equilibrium iso-

Multicomponent nonlinearchromatography

behavior of the concentration path for certain cues. This, however, is believed to be caused by its mathematical properties rather than being based on physical grounds. The results obtained from the simulations of the multicomponent Langmuir isotherm and its two proposed extensions, i.e. polynomial and exponential, provide an interesting starting point for further work in

therms. A major concern which especiallyapplies to modelling and simulation of multicomponent chromatography is to work with thermodynamically consistent multicomponent equilibrium isotherms.Although some of the multicomponent isothermspresentedare empirical they fulfilat leastbasic thermodynamic requirements. The use of the Fritz-Schl|InderIsotherm should be discouraged in this context since it shows unusual

0.048 F ........................... ~ ' " ' " ' " " 0'046 /

/ J ,

277

"'"

, f , ,,.-'"'°'"

C1

0.0440.042/ " "

0.04

0.03%

""'"'Y'""'"

/." ....

C1

o.&s

o.bl 0.615 0.~

o.&s

o.b3 0.635 0.~

0:;45 0.05

c4

0.05 0.045 ................................................................................ 0.04 ¢1

0.035 0.03 0.025 C2

o.o2 o.o~s

o.~os o.b~ 0.615 o.G2 o.&s

o.b3 0.635

o.G4 o.&5

0.05

¢3

0.05 [

.

.

.

.

.

.

.

.

0.045 V ..................................................... 0.04 ............................

........................................................

0.035 f

0.03 ........................................... Cl

0.025 0.02 0.015 0.01 0.oo%

C3

O.()05

i 0.01

t 0.015

o.b2

0.625

0.~3

0.635

i 0.04

0.045

C2 Fig. 13. Case (a): Hodograph representation of the undisturbed simple waves of the exponential extension of the multicomponent Langrnuir isotherm. (a) Complete elution of component 1 (C J). (b) Complete elmion of component 2 (C~). (c) Complete elution of componem 3 (C3). Components: 1 (--), 2 (- - -), 3 (-. -).

278

R. ZENHA,USERNand D. W. T. RIppIN Table6. Equilibrium parameters of the FritzSchlUnderisothermat 293 K 2-Butanol ao b0 atj a~ a~j b~j b2j b3j

5.36 1.13 1 0.25 0.03 0.73 0.88 0.29

a,

o

t-Amylalkohol Phenol 12.16 1.18 0.6 1 0.02 0.79 0.83 0.24 o

this field. From the behavior of the concentration path in the hodograph space initial guesses about the appropriate form of isotherm are possible. This is of particular interest in the estimation of multicomponentequilibrium isotherms where sets of mobile phase and stationary phase concentrations are used for the estimation of equilibrium parameters of an assumed equilibrium relationship. This work has presented theoretical concepts of multicomponent chromatography. Consequently the experimental verification will be the very next step in this field. The following topics would be of special interest • Development of criteria under which conditions of local equilibrium theory could be applied. This includes questions regarding operating conditions, requirements for the hardware, achievable accuracy, etc. • Since the approach is based on the availability of the individual peak shape it has to be verified under which conditions one can measure these peaks in practice.

52 0.22 20.2 17.8 1 0.48 0.54 0 o

Sometimes, however, only the summation curve of the individual components is available. For such cases methods of peak recognition as developed, e.g., in chemometrics could be investigated. Most of these methods, however, have not yet been applied to nonlinear peak shapes. The practical realization of preparative chromatography need not be implemented only by pulse chromatography as presented in this study. Other processes such as displacement chromatography, gradient chromatography, countercurrent chromatography and the recycling of unseparated portions could represent more efficient applications. The numerical algorithm for various multicomponent nonlinear equilibrium isotherms can also be applied to such systems with only minor changes (Rhee et al., 1989). A more conceptual scope for the future would be to define a common set of multicomponent nonlinear equilibrium isotherms on which modelling and simulation of such systems can be based. These isotherms should meet at least basic thermodynamic and mathe-

60 50 40 30 20 10

0

0

5

1'o

1;

2'0

2'5

30

Distance x Fig. 14. Physicalplane diagramfor the separationof a three componentmixtureusingthe Fritz-SchlUnderisotherm:Simplew a v e s (--), s h o c k w a v e s ( - - -).

Multicomponent nonlinear chromatography

279

2500

L=25 2000

1500 c-4 1000

500

0

0

20

40

60

80

100

120

140

Time t Fig. 15. Breakthrough histories for the Fritz-Schltlnder isotherm Breakthrough history profiles. Cases: (a) (--), (b) (---), (c) ("3. matical requirements as introduced in this paper. By these means the comparison between different separation processes as well as the availability o f sufficient equilibrium data could be considerably improved.

List of Symbols

A c C K L

cross sectional area fluid phase concentration simple wave equilibrium constant o f sorption column length

n N r t S u x z Z a

stationary phase concentration total number o f components right eigenvector dimensionless time defined as t=(u¢)l(Z) shock wave velocity dimensionless distance defined as z=(z)/(Z) distance characteristic length selectivity factor void fraction continuous concentration path in the hodograph space

F

25

20 .............. ~.............. + ...............i................~................i................~ ..............!.......... i................i..............

15

10

0

0

10

20

30

40

50

60

70

80

90

I00

N u m b e r o f t r a j e c t o r i e s / s i m p l e wave ( n p t s )

Fig. 16. Error E as function of trajectories assumed in simple waves: classical (*). polynomial extension (O), case (a) (+) and case (b) ( x ) of the exponential form of the multicomponent isotherm.

R. ZENH,KUSERNand D. W. T. RlPPlN

280 A

4, (7"

7" tO

phase constant defined as r/=(1 - e)/(e) eigenvalue concentration in the eigenvector space slope of a characteristic in the physical plane discontinuous concentration path in the hodograph space real time characteristic parameter

Subscripts m Mobile phase s Stationary phase i Component i

Supe~c~p~

f 0 e

k l m r s u

feed state initial state final state wave corresponding to component k lower boundary or left side of a shock wave corresponding to component m right side of a shock shock upper boundary

References Balzli, M. W. (1977) Einsatz von Aktivkohle zur Reinigung einer Mehrkomponenten-Chemieabwassers. Ph.D. thesis, ETH ZUrich. Balzli, M. W., Liapis, A. I. and Rippin, D. W. T. (1978) Applications of mathematical modelling to the simulation of muiticomponent adsorption in activated carbon columns. Trans. IChemE. 56, 145-156. Basmadjian, D. and Coroyannakis, P. (1987) Equilibrium theory revisited. Isothermal fixed-bed sorption of binary systems--I. Solutes obeying the binary Langmuir isotherm. Chem. Eng. Sci. 42, 1723-1736. Basmadjian, D., Coroyannakis, P. and Karayannopoulos, C. (1987) Equilibrium theory revisited. Isothermal fixed-bed sorption of binary systems--II. NonLangmuir solutes with type I parent isotherms azeotropic systems. Chem. Eng. Sci. 42, 1737-1752. Basmadjian, D., Coroyannakis, P. and Karayannopoulos, C. (1987) Equilibrium theory revisited. Isothermal fixed-bed sorption of binary systems--III. Solutes with type I, II and IV parent isotherms phase separation phenomena. Chem. Eng. Sci. 42, 1753-1764. Chiang, A. S. T. and Hwang, G.-W. (1989) Simulation of breakthrough curves by a moving zone collocation method. Computers Chem. Engng 13, 281-290. Courant, R. and Friedrichs, K. O. (1952) Supersonic Flow and Shock Waves. Interscience, New York. Courant, R. and Friedrichs, K. O. (1976) Supersonic Flow and Shock Waves. Springer Verlag, New York. Courant, R., Isaacson, E. and Rees, M. (1952) On the

solution of nonlinear hyperbolic differential equations by finite differences. Comm. Put. App. Math. 5, 243-255. Fritz, W. and SchlUnder, E. U. (1974) Simultaneous adsorption equilibria of organic solutes in dilute aqueous solutions on activated carbon. Chem. Eng. Sci. 29, 127%1282. Fritz, W. and Schltinder, E. U. (1981) Competitive adsorption of two dissolved organics onto activated carbon--I. Chem. Eng. Sci. 36, 721-730. Fritz, W. and Schltlnder, E. U. (1981) Competitive adsorption of two dissolved organics onto activated carbon--II. Chem. Eng. Sci. 36, 731--741. Glueckauf, E. (1949) Theory of chromatography: Part VII. The general theory of two solutes following non-linear isotherms. Discuss. Faraday. Soc. 7, 12-25. Golshan-Shirazi, S. and Guiochon, G. (1988) Comparison between experimental and theoretical band profiles in nonlinear liquid chromatography with a binary mobile phase. Anal. Chem. 60, 2634-2641. Golshan-Shirazi, S. and Guiochon, G. (1989) Experimental characterization of the elution profiles of high concentration chromatographic bands using the analytical solution of the ideal model. Anal. Chem. 61, 462--467. Golshan-Shirazi, S., Lin, B. and Guiochon, G. (1989) Effect of mass-transfer kinetics on the elution of a binary mixture in nonlinear liquid chromatography. J. Phys. Chem 93, 6871-6880. Harten, A. (1983) High resolution schemes for hyperbolic conservation laws. J. Comp. Phys. 49, 357-393. Helfferich, E and Klein, G. (1970) Multicomponent Chromatography. Marcel Dekker, New York. Helfferich, E G. (1986) Theory of multicomponent chromatography--a state-of-the-art report. J. Chromatogr. 373, 45-60. Katti, A. M., Ma, Z. and Guiochon, G. (1990) Prediction of binary, overloaded elution profiles using the simple wave effect. AIChE J. 36, 1722-1730. Lax, P.D. (1957) Hyperbolic systems of conservation laws II. Comm. Pur. App. Math. 10, 537-566. Liu, T.-P. (1975) The Riemann problem for general systems of conservation laws. J. Differenial Equations 18, 218-234. Loureiro, J. M. and Rodrigues, A. E. (1991) Two solution methods for hyperbolic systems of partial differential equations in chemical engineering. Chem. Eng. Sci. 46, 3259-3267. Loureiro, J. M., Costa, C. and Rodrigues, A. E. (1983) Propagation of concentration waves in fixed-bed adsorptive reactors. The Chemical Engineering Journal 27, 135-148. Marquardt, W. (1989) Wellenausbreitung in verfahrenstechnischen Prozessen. Chem. lng Tech. 61, 362-377. Rhee, H. K. (1968) Studies on the Theory of Chromatography. Ph.D. thesis, University of Minnesota. Rhee, H. K. and Amundson, N. R. (1972)A study of the shock layer in nonequilibrium exchange systems. Chem. Eng. Sci. 27, 199-211. Rhee, H. K., Aris, R. and Amundson, N. R. (1970) On the theory of multicomponent chromatography. Phil. Trans. Roy. Soc. A 267, 419-455. Rhee, H. K., Bodin, E and Amundson, N. R. (1971) A

Multicomponent nonlinear chromatography study of the shock layer in equilibrium exchange systems. Chem. Eng. Sci. 26, 1571-1580. Rhee, H. K., Aris, R. and Amundson, N. R. (1986) FirstOrder Partial Differential Equations, Volume I, Theory and Application of Single Equations. Prentice Hall Inc., Englewood Cliffs, NJ. Rhee, H. K., Aris, R. and Amundson, N. R. (1989) FirstOrder Partial Differential Equations, Volume II, Theory and Application of Hyperbolic Systems and Quasilinear Equations. Prentice Hall Inc., Englewood Cliffs, NJ. Schlechter, R. S., Bryant, S. L. and Lake, L. W. (1987) Isotherm-free chromatography propagation of precipitation/dissolution waves. Chem. Eng. Comm. 58, 353-376.

281

Ton~ur, D. (1987) Unifying concepts in non-linear unsteady processes. Part I, Solitary travelling waves. Chem. Eng. Proc. 21, 167-178. Tondeur, D. and Bailly, M. (1987) Unifying concepts in non-linear unsteady processes. Part II, Multicomponent waves, competition and diffusion. Chem. Eng. Proc. 22, 91-105. Vault, D. De (1943) The theory of chromatography. J. Am. Chem. Soc. 65, 532-540. Zenhllusem, R. (1987) Modellierung chromatographischer Trennprozesse. Diplomar~it, ETH ZUrich. Zenhltusem, R. (1993) Modelling and Simulation of Multicomponent Nonlinear Chromatography. Ph.D. thesis, ETH Z{hich.