Ind. Eng. Chem. Res. 2009, 48, 7733–7752
7733
Nonclassical Composition Fronts in Nonlinear Chromatography: Delta-Shock Marco Mazzotti* Institute of Process Engineering, ETH Zurich, Sonneggstrasse 3, CH-8092 Zurich, Switzerland
The local equilibrium theory of chromatography is applied to study binary systems subject to the mixed generalized Langmuir isotherm, ni ) Hici/(1 - K1c1 + K2c2) (i ) 1, 2), where components 1 and 2 are anti-Langmuir-like and Langmuir-like, respectively. Riemann problems, where a stream of a given composition is fed to a column initially saturated at a different composition, are analyzed. Within the hyperbolic region of the hodograph plane close to the origin, it is shown that such Riemann problems admit not only classical solutions, consisting of constant states separated by simple wave and shock wave transitions, but also nonclassical solutions, namely including a continuous nonsimple wave transition and a delta-shock. The latter can be viewed as a traveling spike of growing strength superimposed to the discontinuity separating the initial and the feed states. Exact conditions for the occurrence of classical and nonclassical solutions are derived. Thus, this work adds two new types of composition fronts to the classical ones already known in nonlinear chromatography. In the case of the delta-shock, the most striking new transition, it is possible to derive exact relationships for its properties, including the slope of its image in the physical plane, i.e. the reciprocal of its propagation speed, dt/dz ) {1/W}{1 + ν((H1K2[n2] - H2K1[n1])/(H1K2[c2] - H2K1[c2]))}, where the symbol [ · ] denotes the jump across the discontinuity of the quantity enclosed. These analytical expressions match the results obtained through numerical simulations. The physical legitimacy of the model considered in this work is demonstrated theoretically and is further supported by the experimental evidence of the existence of delta-shocks that has been recently provided [Mazzotti et al. J. Chromatogr., submitted]. 1. Introduction Operating a chromatographic column typically involves adsorbing a mixture of more and less retained chemical species onto an empty column, i.e. saturated with the mobile phase only, and then, after the column becomes saturated with the feed mixture, i.e. the feed mixture reaches the column outlet, regenerating the column by desorbing the feed mixture using the same initial mobile phase. During both adsorption and desorption composition fronts develop and propagate along the column, according to a pattern that is determined by the adsorption isotherm that characterizes the retention behavior of the species involved. The adsorption isotherm establishes a relation between fluid phase and adsorbed phase concentrations and reflects the phase equilibrium thermodynamics of the system. It is worth noting that pulse injection, i.e. an operation mode where a relatively small amount of the feed mixture is injected before switching back to the initial mobile phase, is a very short adsorption/desorption cycle where composition fronts interact before reaching the column end. In this paper, we will not consider situations where such interactions between composition fronts occur. In the limit of an extremely small injected amount of feed, pulse injection is the operation mode of analytical chromatography, where chemical species are so diluted that their retention is subject to a simple linear isotherm. In the following, we will consider only the much richer behavior of systems subject to nonlinear adsorption isotherms. The local equilibrium theory of nonlinear chromatography is the best modeling approach to understand and predict the dynamics of composition fronts in chromatographic operations. This is a model where mass transfer resistance and axial dispersion are neglected, and the behavior of the system is controlled by convection and partition between fluid and adsorbed phases, according to the adsorption isotherm; ther* To whom correspondence should be addressed. E-mail: marcom@ ethz.ch. Tel.: +41-44-6322456. Fax: +41-44-6321141.
modynamic equilibrium is assumed to be attained at any time and position in the column. In the frame of a mathematical description of chromatographic operations, adsorption and desorption are special cases of a more general situation whereby the column is at an initial state, corresponding to an initial composition, and a new state, corresponding to the feed composition, is fed to the column starting at time zero. This is a so-called Riemann problem and will be the subject of this paper. Three types of composition fronts, or waves, have been so far observed in nonlinear chromatography and described by the local equilibrium theory, namely a rarefaction (continuous) wave, a shock (discontinuous) wave and a semishock, which consists of a combination of the previous two types. These are observed in single component, as well as in multicomponent chromatography, and the conditions for the occurrence of one or the other are well understood; in the following, we will call these classical composition fronts.2,3 In this paper, we aim at showing that there are also nonclassical, or degenerate, or singular, composition fronts that may occur in nonlinear chromatography. We will discuss their prerequisites and behavior in the case of a binary system subject to a rather simple isotherm, i.e. the mixed Langmuir isotherm M2, that has been analyzed in the frame of a recent study applying equilibrium theory to a family of generalized Langmuir isotherms.4 These are given by ni )
NiKici 1 + p1K1c1 + p2K2c2
(i ) 1, 2)
(1)
where ci and ni are fluid and adsorbed phase concentrations of species i, respectively; Ki and Ni are the equilibrium constant and the saturation loading capacity of the ith component, respectively; Hi ) NiKi is its Henry’s constant, i.e., the infinite dilution slope of its isotherm, where H2 > H1; p1 and p2 can take the values (1 and define the sign of the corresponding
10.1021/ie9001537 CCC: $40.75 2009 American Chemical Society Published on Web 07/15/2009
7734
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
Figure 1. (a) Characteristic fields in the hodograph plane (u1, u2), together with the parabola, which is their envelope, for the system where a2 ) 2a1 ) 3. The dotted line is the locus where d ) 0. For the sake of simplicity, only regions 1 and 4 are shown (the complete hodograph plane is shown in an earlier paper4). (b) Characteristic fields in the characteristic parameter plane (ω1, ω2). Note that the point of coordinates (a1, a2) corresponds to the origin of the hodograph plane where u1 ) u2 ) 0; the segments ω1 ) a1 and ω2 ) a2 map onto u1 ) 0 and u2 ) 0, respectively; the vertical and the horizontal dotted lines correspond to Γ2 and to Γ1 characteristics, respectively; the portion of the diagonal between a1 and a2 maps onto the portion of the envelope belonging to region 1 of the hodograph plane.
term in the denominator. The mixed Langmuir isotherm M2 is the case where p1 ) -1, i.e. the less retained component exhibits an anti-Langmuir character, and p2 ) 1, i.e. the more retained species exhibits a Langmuir character. Therefore, the adsorption of species 2 is competitive, whereas that of species 1 is cooperative. We will consider only this case in this work. Solutions similar to the traVeling spike that will be presented here for the first time have been obtained in very few cases arising in fields other than chromatography.5-7 These have been called delta-shock or singular shock, and they are not amenable to a general mathematical treatment.8-10 In this paper, we will exploit the methods of the local equilibrium theory of chromatography2-4 to prove the existence of the nonclassical composition fronts and to provide exact conditions for their occurrence in the case of the mixed Langmuir isotherm M2. Then, we will focus on the delta-shock solution and derive exact formulas for its speed and strength by using and extending a mathematical technique proposed earlier.8 2. Background 2.1. Model Equations. In the frame of the local equilibrium theory, assuming isothermal conditions, diluted streams (constant velocities), and homogeneous packing (constant bed void fractions) with no internal particle porosity, the dimensionless mass balance equations for a fixed bed chromatographic column are as follows: ∂ci ∂ (ci + νni) + W ) 0 (i ) 1, 2) ∂t ∂z
(2)
where t and z are the time and space coordinates (with L being the bed length); W is the interstitial fluid velocity; ν is the phase ratio, i.e. solid phase volume to fluid phase volume of the column; and the adsorbed phase concentrations ni are given in terms of the fluid phase composition through the mixed Langmuir isotherm M2, i.e. the adsorption isotherm of eq 1, with p1 ) -1 and p2 ) 1.
For the sake of generality the following dimensionless compact form of the equations above will be used in the following: ∂ui ∂ (ui + Vi) + ) 0 (i ) 1, 2) ∂τ ∂x
(3)
with Vi )
aiui aiui ) d 1 - u 1 + u2
(i ) 1, 2)
(4)
where only the pairs (u1, u2) that make the denominator d ) 1 - u1 + u2 be positive are allowed, i.e. the points in the first quadrant of the hodograph plane (u1, u2) that are below the straight line of unitary slope that crosses the axis u2 ) 0 at u1 ) 1 (see Figure 1a). In the equations above ui ) Kici and Vi ) νKini are the ith component’s dimensionless fluid and solid phase concentrations, respectively; ai ) νNiKi (with a2 > a1) is its Henry’s constant multiplied by the phase ratio, i.e. what is called the retention factor and indicated as k′i in chromatography; τ ) tW/L and x ) z/L are the dimensionless time and space coordinates, respectively. Wherever possible eq 3 will be solved using the method of characteristics. Such equilibrium theory solution will be compared with the numerical solution of a regularization of eq 3 obtained by adding a small dispersion term on its right-hand side, which is proportional to the second-order space derivative of the unknown: ∂ui ∂2ui ∂ (ui + Vi) + ) εi 2 ∂τ ∂x ∂x
(i ) 1, 2)
(5)
where the small parameter εi is the reciprocal of a properly defined Peclet number. This is the equilibrium dispersive model of chromatography.11 2.2. Characteristics in the Physical and in the Hodograph Plane. The mass balance eqs 3 with eqs 4 constitute a system of two first-order, homogeneous, and reducible partial
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
differential equations, which has to be completed by defining an initial value problem, e.g. a Riemann problem. The solution of such a problem can be obtained through the method of characteristics, whose application to the problem at hand will only be summarized here because all the details can be found elsewhere.3,4 It can be shown that there are two corresponding sets of characteristics, each consisting of two families identified by the index j ) 1, 2, i.e. the characteristics Cj in the physical plane (x, τ), which are defined as Cj:
dτ ) σj ) 1 + θj (j ) 1, 2) dx
(6)
and the characteristics Γj in the hodograph plane (u1, u2), which are defined by Γj:
du1 θj - V22 ) du2 V21
(j ) 1, 2)
(7)
In the equations above, the variables θj are the eigenvalues of the matrix V ) [Vij], where Vij ) ∂Vi/∂uj is the partial derivative of the adsorption isotherm 4, i.e.: V11 ) a1(1 + u2)/d2
(8)
V22 ) a2(1 - u1)/d2
(9)
V12 ) -a1u1 /d2
(10)
V21 ) a2u2 /d2
(11)
Hence, the eigenvalues θj are the roots of the quadratic equation (with θ1 < θ2): θ - (V11 + V22)θ + V11V22 - V12V21 ) 0 2
7735
the point of tangency. The points belonging to the convex partition of the hodograph plane defined by the envelope (13) are physically meaningful, since they lie in the first quadrant and do not violate the constraint d > 0. However, in these points, the system of partial differential eqs 3 with 4 is not hyperbolic and the method of characteristics cannot be applied. It is worth noting that region 1 is the region to which the origin belongs, as well as a portion of the diagonal u1 ) u2. For most applications, this is the domain of interest, and this is the one that we will consider in the following. In this region characteristics of both families, i.e. Γ1 and Γ2, have negative slopes as shown in Figure 1a. Characteristics Γ1 are steeper (in absolute terms), cross the horizontal axis, and include the vertical axis. Characteristics Γ2 are less steep, cross the vertical axis, and include the horizontal axis. Note that the two characteristics coincide only for points on the envelope. 2.3. Characteristic Parameters. It can be demonstrated that within the hyperbolic part of the hodograph plane there exists a one-to-one mapping between (u1, u2) and the pair of characteristic parameters (ω1, ω2), which are obtained as roots of the following quadratic equation:3,4 (1 - u1 + u2)ω2 - [a1(1 + u2) + a2(1 - u1)]ω + a1a2 ) 0 (14) The discriminant of this equation is the same as that of eq 12; hence, there are two real and distinct roots when the system of partial differential eqs 3 with 4 is hyperbolic.4 The four hyperbolic regions in the hodograph plane correspond to four different ranges of ω1 and ω2 values, namely: region 1: a1 e ω1 < ω2 e a2
(15)
region 2: 0 e ω1 < ω2 e a1
(16)
(12)
If they are real and distinct, then system 3 is hyperbolic. This is true when the discriminant of eq 12 is positive, i.e. when: ∆ ) [a1(1 + u2) + a2(1 - u1)]2 - 4a1a2(1 - u1 + u2) > 0 (13) The families of characteristics Γj (j ) 1, 2) in the hodograph plane form a network of straight lines.3,4 Since the coefficients of eq 12 depend only on u1 and u2, then such network does not depend on the specific initial value problem and can be drawn once and for all for a given adsorption isotherm 4, as shown in Figure 1a. On the contrary, the characteristics Cj in the physical plane are different for different initial value problems; along Cj characteristics the composition state (u1, u2) remains constant, and therefore, they are also straight lines. The characteristics Γj are tangent to the envelope defined by eq 13. This is a parabola that lies in the first quadrant of the hodograph plane and is tangent to the vertical and the horizontal axes in (0, 1 a1/a2) and (a2/a1 - 1, 0), respectively, and to the d ) 0 line in u1 ) a2/(a2 - a1), u2 ) a1/(a2 - a1) (see Figure 1a). Therefore, there are always four disjoint regions where the system we are considering is hyperbolic. Region 1 is the finite region between the envelope, the vertical axis for 0 e u1 e (1 - a1/a2), and the horizontal axis for 0 e u2 e (a2/a1 - 1). Region 2 is the unbounded region between the parabola and the horizontal axis for (a2/a1 - 1) e u2 < +∞; region 3 is the infinite region between the parabola and the line d ) 0 above the point of tangency; region 4 is the region between the vertical axis for (1 - a1/a2) e u1 < 1, the parabola, and the d ) 0 line below
-∞ e ω1 < ω2 e 0
(17)
region 4: a2 < ω1 < ω2 e +∞
(18)
region 3:
Note that this is different from all other generalized Langmuir cases, where the intervals to which ω1 and ω2 belong are disjoint.4 The triangular region defined by eq 15 that maps onto region 1 in the hodograph plane is shown in the so-called characteristic plane with coordinates (ω1, ω2) in Figure 1b. Characteristics Γj in the hodograph plane map on segments parallel to the axes in the (ω1, ω2) plane, namely characteristics Γ1 on segments with ω2 ) const, and characteristics Γ2 on segments with ω1 ) const. The segment where ω2 ) a2 corresponds to u2 ) 0, whereas that where ω1 ) a1 maps onto u1 ) 0; therefore, the origin of the hodograph plane maps onto (ω1 ) a1, ω2 ) a2); the parabolic envelope in region 1 maps onto a1 e ω1 ) ω2 e a2. The inverse relationships of eq 14 yield the fluid and adsorbed phase concentrations in terms of ω1 and ω2:4 u1 )
a2(ω1 - a1)(ω2 - a1) ω1ω2(a2 - a1)
(19)
7736
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
a1(ω1 - a2)(ω2 - a2) ω1ω2(a2 - a1)
(20)
a1a2 ω1ω2
(21)
V1 )
(ω1 - a1)(ω2 - a1) a2 - a1
Figure 2. Propagation of a shock between states L and R. The solid and the dashed lines represent the concentration profiles of either species at time τ and at time τ + ∆τ, respectively.
(22)
V2 )
(ω1 - a2)(ω2 - a2) a2 - a1
(23)
whose slopes are given by eq 6 where θj ) ωj2ω3-j/(a1a2), with ωRj e ωj e ωLj and ω3-j constant and equal to the value that is common to states L and R. When the conditions above are not fulfilled, i.e. when the Cj characteristics fan counterclockwise when going from left to right in the physical plane, or equivalently when ωj increases moving from left to right through the jth transition, then the transition is a discontinuous shock indicated as Σj in the hodograph plane. The image of Σj in the physical plane is a single straight line through the origin, called Sj, that propagates at constant velocity along the column and separates the two constant states. The dimensionless propagation velocity of a shock of type j, i.e. λ˜ j, can be obtained by solving a material balance across the shock assuming that the discontinuity itself has no volume and no capacity as shown in Figure 2.3,4 This yields the following equation:
u2 )
d)
It can also be shown that the eigenvalues θj of the matrix V ) [Vij] can be recast in terms of ω1 and ω2 as follows:4 ω1 ω12ω2 ) θ1 ) a1a2 d
(24)
ω2 ω1ω22 ) a1a2 d
(25)
θ2 )
Since d is always positive and ω1 < ω2, the sign of ωj and θj is the same; hence, θ1 < θ2, as required. Substituting the last two equations into eq 6, the slope of the characteristics Cj in the physical plane are obtained in terms of ω1 and ω2, and one sees that σ1 < σ2 always. Since along Cj the composition state is constant, so are the two characteristic parameters and the slope of the characteristics, which are in fact straight lines. 2.4. Riemann Problems. Let us consider the typical Riemann problem for eq 3, where there are assigned an initial state (u1B, u2B) for x g 0 and τ ) 0 and an inlet state (u1A, u2A) for τ g 0 and x ) 0, which are separated in the origin by a discontinuity. The classical solutions of such problems consist of three constant states, i.e. regions of the physical plane where the solution is constant, separated by two simple wave transitions, i.e. composition fronts traveling along the column.3,4 Two constant states coincide with the initial and the inlet states, whereas the intermediate state is defined by the intercept (u1I , u2I ) in the hodograph plane of the Γ1 and Γ2 characteristics originated in the initial and inlet state, respectively. Note that, due to the properties of the one-to-one mapping of the hodograph plane onto the (ω1, ω2) plane, if the initial state B and the inlet state A map onto (ω1B, ω2B) and (ω1A, ω2A), respectively, then the intermediate state I maps onto (ω1A, ω2B) and the corresponding composition can be calculated using eqs 19 and 20. The two transitions in the physical plane have their image in the hodograph plane along the segment of Γ2 characteristic connecting the inlet state and the intermediate state and along the segment of Γ1 between the intermediate state and the initial state. Each transition can be either continuous or discontinuous. The Γj simple wave transition is continuous and consists of straight Cj characteristics emanating from the origin with slope σj given by eq 6 if and only if the Cj characteristics fan clockwise when going from the constant state on the left-hand side of the transition (indicated as state L in the following) to that on the right-hand side of it (state R). Due to eqs 24 and 25, this condition holds when ωj decreases from left to right through the Γj transition. Therefore, the second transition is a continuous Γ2 simple wave if the inlet state is above the intermediate state in the (ω1, ω2) plane, whereas there is a continuous Γ1 simple wave if the initial state is on the left with respect to the intermediate state (refer to Figure 1b). Summarizing, a Γj simple wave maps into a pencil of characteristics in the physical plane
[ui + Vi]∆x ) [ui]∆τ (i ) 1, 2)
(26)
where the symbol [ · ] denotes the jump across the discontinuity of the quantity enclosed (left quantity minus right one). By taking the limit for ∆τ f 0, ∆x/∆τ f λ˜ j, the following equation is obtained: σ˜ j )
[V2] [V1] 1 )1+ ) 1 + θ˜ j ) 1 + [u1] [u2] λ˜
(27)
j
which incorporates the so-called Rankine-Hugoniot condition. The quantity θ˜ j can be given in terms of ω1 and ω2 as follows, where the superscripts L and R refer to the states on the left and on the right-hand side of the discontinuity in the physical plane, i.e. either the inlet and intermediate state in the case of S2, or the intermediate and the initial state in the case of S1: θ˜ 1 )
ωL1 ωL1 ωR1 ω2 ωR1 ) R ) L a1a2 d d
(28)
θ˜ 2 )
ωR2 ω1ωL2 ωR2 ωL2 ) R ) L a1a2 d d
(29)
It is worth reminding that ω1 is constant along a Γ2 characteristic and, therefore, across a S2 shock; therefore, ω1 carries no superscript in eq 29. Likewise ω2 is constant on Γ1 and across a S1 shock; hence, it has no superscript in eq 28. A shock of type j can be described by the following Heavyside generalized function (the Heavyside function H(x) is 0 for negative and 1 for positive values of x): ui(τ, x) ) uLi - [ui]H(x - λ˜ jτ) (i ) 1, 2)
(30)
Note that across an S1 shock, due to the condition ω1L < ω1R, the following inequalities apply: σL1 < σ˜ 1 < σR1 < σR2
(31)
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
σ˜ 1
T and for |x| > X. For smooth solutions, the weak formulation is equivalent to the original PDE (3). Note also that for the sake of generality, the last equation has been written for a Riemann problem defined on the whole x axis, i.e.: at τ ) 0, x g 0: ui ) uBi (state B)
(47)
at τ ) 0, x < 0: ui ) uAi (state A)
(48)
where states A and B fulfill the conditions for the occurrence of a delta-shock specified in section 3.2.2. This Riemann problem is equivalent to the one defined earlier, i.e. eqs 35 and 36, since the solution of the two formulations of the Riemann problem in the domain of physical interest, i.e. τ g 0 and 0 e
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
x e 1, is exactly the same. Unfortunately, not even the weak formulation is viable, since the theory of generalized functions is a linear theory that does not allow computing Vi through the nonlinear function (4) when ui is the linear combination of the Heavyside generalized function and of the Dirac’s delta function given by eq 45. One can readily see that the division of ui by d ) 1 - u1 + u2 required by eq 4 is not amenable to any straightforward definition. This implies that even writing Vi formally as Vi(τ, x) ) VLi - [Vi]H(x - sτ) + Viτδ(x - sτ)
(49)
(which would be consistent with the calculated profiles shown in Figure 11d) would not help. We would in fact be able to use this expression in the weak formulation (46), but we would still be unable to relate the quantities Vi with the quantities Ui due to the impossible division operation that this would require. 4.2.2. Background on Delta-Shock. The considerations above demonstrate that the major difficulty is to define a proper notion of weak solution that can apply to the singular solution ui given by eq 45. The approaches proposed in the literature are all based on some kind of regularization, either of the original partial differential equation, or of the solution, or of the notion of equality itself. These three approaches are briefly described in the following, where references providing the necessary details are given as well. The first approach entails studying the regularized eq 5 (with ε ) εi , 1), finding an analytical solution depending on the small parameter ε through for instance perturbation methods, and determining the asymptotic properties of such a solution for ε f 0.6-8 This approach is of course very problem specific. Examples of its successful application to cases yielding deltashock solutions are the Keyfitz-Kranzer equation:6,8 uτ + (u2 - V)x ) 0
Vτ + (u3 /3 - u)x ) 0
(50)
and the following system constituting a one-dimensional zeropressure gas dynamics model:7 uτ + (u2)x ) 0
Vτ + (uV)x ) 0
(51)
where u and V are the unknowns and the subscripts τ and x indicate partial differentiation with respect to time and space. In the second approach, one defines a regularized solution that depends on a small parameter ε and that approaches the singular solution as this parameter approaches zero. As all operations involved in the original PDE are doable for the regularized solution, one can calculate the PDE residues and then enforce the condition that they approach zero as ε f 0. Equation 50 has been solved in this manner using both a smooth approximation and a simple step-function approximation (boxfunction),8 whereas eq 51 has been solved using a similar method called the weak-asymptotics method.10,16 It is worth noting that how to choose the functional form of the regularized solution is indeed problem specific and requires a good understanding of the physics behind the original conservation laws. The third approach is based on Colombeau’s theory of generalized distributions,17 where the concept of equality ()) is generalized to introduce that of association (≈). According to the latter, two generalized functions are associated when the integrals of the product of each function times a smooth test function (as in eq 67 below) are the same. Let us consider as an example the bell function, which is continuous on the real axis with all its derivatives:
B(x) )
{
(
1 1 exp 2 β x -1 0
)
if |x| < 1
7747
(52)
if |x| g 1
∞ B(x) dx where β is such as to fulfill the integral condition ∫-∞ ) 1, i.e. β = 0.444. Then, obviously B(x) * δ(x), but according to Colombeau’s theory, B(x) ≈ δ(x) due to eq 59 below (where η ) 0). Colombeau distinguishes between a microscopic aspect of a function, and its macroscopic aspect.17-19 Functions with different microscopic aspects may have the same macroscopic aspect. The former corresponds to the genuine nature of the physical object described by that function, since any discontinuous variation in a quantity representing a physical property occurs in a very small but nonzero length scale hence it should be described by a smooth function with a very steep variation and not by a step function. On the contrary, the macroscopic aspect of it, as captured for instance by a step function, can be viewed as the mathematical idealization of that physical object, which is valuable as it enables a much simpler mathematical treatment and captures some higher level properties of the physical object.17 In the case of the functions used here as examples, B(x) and δ(x) do not have the same microscopic aspect, hence, B(x) * δ(x), but they have the same macroscopic aspect, thus B(x) ≈ δ(x). Through these new definitions the space of generalized distributions is endowed with an algebra that allows for operations such as the multiplication of distributions. Through this approach and the use of the so-called splitted delta functions, eqs 50 and 51 and extensions thereof have been solved.9,20 Both the second and the third approach, i.e. using either box functions or splitted delta functions, can be applied to determine the unknown quantities in eq 45 as will be discussed elsewhere.21 In this paper, for the sake of brevity, we apply only the second approach, using in particular a smooth approximation of ui, i.e. following a similar approach as Keyfitz and Kranzer used in section 2.2 of their 1995 paper to solve eq 50.8 4.2.3. Definitions. First, some definitions are given that will be useful in the following. Let us consider a smooth function b(x) satisfying b(x) g 0 ∀ ∞ b(x) x ∈ R, b(x) ) 0 ∀ |x| > 1 (i.e., supp(b) ⊆ [ -1,1]), and ∫-∞ dx ) 1. A possible choice is the bell function B(x). Let us also define the one parameter family
1 x bε(x) ) b ε ε
()
(53)
with ε > 0, which as ε approaches zero constitutes a delta sequence, i.e. a sequence of functions that approaches the Dirac’s delta function as shown by eq 59 with η ) 0. Let us then define (see Figure 13): ξ ) x - sτ
(54)
bε0(ξ) ) bε(ξ)
(55)
ε b(ξ) ) bε(ξ + ε)
(56)
ε b+ (ξ) ) bε(ξ - ε)
(57)
hε(ξ) )
∫
ξ
-∞
bε0(y) dy
(58)
Each family bεk (k ) -, 0, +) constitutes a delta sequence as ε f 0, thanks to the following observation (η is -1, 0, +1 for the different k labels):
7748
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
|∫ |∫ 1 ε
ε(1+η)
-ε(1+η)
1
-1
( εx - η)ψ(x) dx - ψ(0)| )
the same, has decisive consequences on the derivation that follows. Using eq 60 to calculate the regularized denominator of the adsorption isotherm yields in fact
b
|
b(y)(ψ(0) + ε(y + η)ψ′(χ)) dy - ψ(0) e 4εsup{bψ′}
ε dε(τ, ξ) ) 1 - uε1 + uε2 ) dL - [d]hε + (uˆ2 - uˆ1)τb-
(59)
(62)
where 0 < χ < ε(y + η) or ε(y + η) < χ < 0. Note that in the frame of Colombeau’s generalized distributions the last result proves that bεk ≈ δ(x), (k ) 0, +, -), as they share the same macroscopic aspect as clarified above, although bε- * b0ε * bε+ because they differ in their microscopic details as shown in Figure 13. Likewise, the one parameter family of functions hε(x) approaches the Heavyside function as ε approaches zero; moreover, dhε/dξ ) b0ε (ξ). 4.2.4. Smooth Approximation of the Delta-Shock Solution. Following a similar approach as Keyfitz and Kranzer,8 and using the definitions of the previous section, the following smooth approximation of ui in eq 45 is proposed:
where its singular part contains only a positive term proportional to the left-hand side of the singular part of uεi . In other words while uεi exhibits two contributions to its singular component that are positioned symmetrically with respect to ξ ) 0, dε is in that respect asymmetric. As both uεi and dε are smooth functions, the regularized adsorbed phase concentration can be calculated using the adsorption isotherm (4), thus resulting in
ε (ξ) + uεi (τ, ξ) ) ujεi (ξ) + u˜εi (τ, ξ) ) uLi - [ui]hε(ξ) + uˆiτbε uˆτb+ (ξ) (60)
where ε is a small positive parameter, uˆi and uˆ are constants, which fulfill the key assumption that uˆ2 > uˆ1 g 0, and the functions ujεi and u˜εi represent the bounded (first and second term) and the singular (third and fourth term) parts of uεi . This definition has one special feature, which is novel with respect to what proposed so far in the literature. The singular part of uεi consists of two terms, which differ in their microscopic details because they have disjoint supports. However, both terms are associated to the Dirac’s delta, i.e. they have the same macroscopic aspect, and therefore, comparing eqs 45 and 60 and considering the behavior of uεi as ε f 0 shows that Ui ) uˆi + uˆ
(61)
Moreover the fact that the first term in the singular part of uεi , i.e. that on the left-hand side of the delta-shock, is different for the two components, whereas that on the right-hand side is
Vεi (τ, ξ) )
{
aiuεi
) Vjεi (ξ) + V˜ εi (τ, ξ) ) {VLi - [Vi]hε} +
dε
ε aiuˆτb+
+
dε
[Vi][d]hε(1 - hε) dε
+
ε τb-
dε
(aiuˆi - (VLi - [Vi]hε))
}
(63) Vjεi
V˜ εi
where the functions and represent the bounded (first and second term) and the singular (remaining terms) part of Vεi . As the denominator dε contains bε- only and not bε+, it will be shown in the next section that only the term in the singular part of Vεi proportional to bε+ contributes to the singular part of the formal solution Vi of eq 49. It can be easily shown that if either both terms proportional to bε- and to bε+ were present in uεi and dε (which would be the case if also the term proportional to bε+ were different for the two species) or only one term were present in uεi and then also in dε, then the singular part of Vi would be zero, which would be in striking contradiction with the numerical evidence provided by the simulations illustrated in Figure 11d. It is also worth noting that the asymmetry exhibited by the singular part of dε is also nicely consistent with its calculated profiles shown in Figure 11c. Thus, summarizing the choice of the regularized function uεi allows accounting for all the important qualitative features of the numerical solution of the problem under consideration. 4.2.5. Conservation Laws in a Moving Coordinate Frame. The proper framework in which to study the delta-shock solution is the moving coordinate frame of reference defined by ξ ) x - sτ, where the conservation eq 3 can be recast as (ui + Vi)τ + ((1 - s)ui - sVi)ξ ) 0
(64)
(the subscripts τ and ξ indicate partial differentiation). The Riemann problem of eqs 47 and 48 is reformulated in the (τ,ξ) space as follows: at τ ) 0, ξ g 0: ui ) uBi (state B)
(65)
at τ ) 0, ξ < 0: ui ) uAi (state A)
(66)
Let us consider now the residues Rεi of the weak formulation (see eq 46) of the PDEs (64) when the regularized solution (60) and the corresponding adsorbed phase concentration (63) are used: Rεi )
∫
∞
0
dτ
∫
∞
-∞
{(uεi + Vεi )ψτ + ((1 - s)uεi - sVεi )ψξ} dξ +
∫
∞
-∞
Figure 13. Plots of the functions defined by eqs 55-58, when calculated with b(x) ) B(x), i.e. the bell function, and ε ) 0.1.
{uεi (0, ξ) + Vεi (0, ξ)}ψ(0, ξ) dξ
(67)
Here, the function ψ(τ,ξ) is a smooth test function having compact support, i.e. which possesses continuous partial deriva-
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
tives of any order, and vanishes for τ > T and for |ξ| > X. Note that in the case of a smooth solution the differential and the integral (weak) formulations are equivalent and that the weak formulation is chosen here for the sake of convenience. Since all functions are smooth, integration by parts can be applied wherever possible; doing this separately for the bounded and the singular parts of uεi and Vεi and considering the compact support of ψ yields Rεi
)-
∫ dτ ∫ ∫ dτ ∫ ∞
0
∞
{(uεi -∞
∞
∞
-∞
0
+
Vεi )τ
+ ((1 -
s)ujεi
-
sVjεi )ξ}ψ
dξ +
((1 - s)u˜εi - sVjεi )ψξ dξ ) -Ψεi + Ξεi
(68)
where the integrand of Ψεi is proportional to ψ and that of Ξεi is proportional to ψξ. Due to this difference, the two integrals have to be analyzed separately in order to find simple relationships (eqs 73 and 75) that guarantee that they vanish as ε approaches zero. In a first reading, one can skip the details that follow and go directly to the beginning of section 4.2.6. Integral Ξεi . Let us use eqs 45 and 49 to rewrite the integral ε Ξi as follows: Ξεi )
∫
∞
0
dτ
∫
∞
-∞
∫
{
∞
0
ε ε (1 - s)(uˆiτb+ uˆτb+ )-s
∫
dτ
∞
-∞
{
ε τbR1 ε
d
ε aiuˆτb+
dε
}
}
ψξ dξ +
+ R2hε(1 - hε) ψξ dξ
(69)
where the symbols Rj represent bounded non-negative functions of ξ and τ, whose detailed expressions can be easily obtained using eqs 45 and 49 (this applies also to eq 74 below). The first double integral consists of terms that approach δ(ξ), i.e. the Dirac’s delta function; hence, their integral in dξ yields finite values. While this is a consequence of eq 59 (with ψξ instead of ψ) for the first two terms multiplied by (1 - s), it has still to be proven for the third term. Let us consider its integral in dξ, which can be written as follows when accounting for the fact that bε+ and bε- have disjoint supports:
∫
∞
ε b+
-∞
dε
ψξ dξ )
∫
ε b+
2ε
0
εf0
dL - [d]hε
ψξ dξ 98
ψξ(0)
) dL - [d]hε(0) 2ψξ(0) (70) dL + dR
where the last limit can be demonstrated following the same approach as in eq 59. The second double integral in eq 69 is O(ε), i.e. it behaves like ε as ε approaches 0. In fact, as all terms are non-negative, the following inequalities hold:
∫
∞
R -∞ 1
∫
∞
-∞
ε τb-
d
ε
ψξ dξ
H2 - H1
(88)
In this very simple case, i.e. component 2 displacing component 1, the initial and feed state lay on the axes of the hodograph plane (see Figure 1a) and are not constrained to be within its region 1 for the delta-shock to be possible (u1 must be smaller than 1, however, to avoid a negative denominator in the adsorption isotherm). This can be clarified by at the same time giving an interesting physical interpretation how the mixed isotherm (4) leads to the delta-shock. First, let us note that the slopes in the physical plane of the shock transitions corresponding to the adsorption of component 2 and to the desorption of component 1 are linear functions of the ratios n2A/c2A and n1B/c1B, respectively. From section 3.2.2, it is known that the condition for the occurrence of a delta-shock is that the desorption front of species 1 should be steeper than the adsorption front of species 2, i.e. that n2A/c2A < n1B/c1B. The last inequality is exactly eq 88 and holds independently of whether states A and B belong to the hyperbolic region 1 of the hodograph plane as long as they contain only species 2 and 1, respectively. Then, note that the ratios n1B/c1B and n2A/c2A increase and decrease, respectively, as c1B and c2A increase. If concentrations are small and eq 88 does not hold, the desorption front of species 1 is less steep (faster) than the adsorption front of species 2. As the two concentration values increase more and more, finally
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
they reach a level where eq 88 is fulfilled and the front of 1 would like to be steeper (slower) than the front of 2. However, due its competitive role the presence of component 2 interfering with component 1 reduces its adsorption and makes its desorption front faster. In turn, the presence of component 1, which enhances component 2 adsorption, slows down its adsorption front. The outcome of these two simultaneous effects taking place at the point in time and space where the two fronts collide is the accumulation of mass at the discontinuity and the development of the traveling spike, i.e. the occurrence of the delta-shock. These remarks confirm that while simple waves, shocks, and semishocks can occur for both single-component and multicomponent chromatographic systems, the delta-shock is a phenomenon that cannot occur in single-component systems. After showing that a delta-shock can occur in binary systems, it remains an open question whether it can in multicomponent systems. Moreover, as a delta-shock is clearly due to the interaction between a competitive and a cooperative species, there is no reason to exclude that it can occur for isotherms other than the mixed M2 isotherm considered here that exhibit the same type of competitive-cooperative interaction. We have just shown that there exists at least one experimental system leading to a delta-shock.1 Why has not this phenomenon been observed earlier in the lab? First, this type of interaction is rather rare, and even when it is observed at low concentrations, it often evolves into a competitive-competitive behavior at higher concentration levels (see, for instance, ref 22). Second, the type of Riemann problem yielding a delta-shock is rather uncommon in chromatography practice. Finally, it is also clear that the indefinite growth of concentration expected at the deltashock front will in practice be limited by the solubility of the compounds involved, thus masking the occurrence of what theory predicts. Further investigation will tell us how often systems with the right characteristics leading to a delta-shock exist in practice, and whether any of the newly discovered nonclassical composition fronts may have practical, possibly useful implications in preparative chromatography. Acknowledgment The author thanks Giuseppe Storti (ETH Zurich and Politecnico di Milano, Italy), Hyun-Ku Rhee (Seoul National University, Korea), Thes Ba¨bler (ETH Zurich), Daniel Tondeur (University of Nancy, France), and Barbara Lee Keyfitz (Ohio State University, OH) for very helpful discussions. Notation ai ) retention factor of component i, ai ) νNiKi b(x) ) function B(x) ) bell function, eq 52 ci ) fluid phase concentration of component i Ce ) continuous nonsimple wave transition in the physical plane Cj ) jth characteristics in the physical plane d ) denominator in the adsorption isotherm, d ) 1 - u1 + u2 dε ) smooth approximation of the adsorption isotherm denominator in the delta-shock, eq 62 hε(ξ) ) regularized Heavyside function; see Figure 13 H(x) ) Heavyside function Hi ) Henry’s constant of component i, Hi ) NiKi k ) very large parameter, eq 44 ki′ ) retention factor of component i, ki′ ) νNiKi Ki ) adsorption equilibrium constant of component i L ) length of the chromatographic column
7751
n ) number of theoretical plates ni ) adsorbed phase concentration of component i Ni ) adsorbed phase saturation concentration of component i pi ) parameter in the generalized Langmuir isotherm, pi ) (1 Riε ) residue of equation i, eq 46 s ) nondimensional speed of propagation of a delta-shock Sj ) jth shock in the physical plane t ) time T ) finite parameter ui ) dimensionless fluid phase concentration of component i, ui ) Kici uˆ ) component of the delta-shock strength uˆi ) component of the delta-shock strength of species i uiε ) smooth approximation of the fluid phase delta-shock solution for component i, eq 60 ujiε ) bounded part of uiε, eq 60 u˜iε ) singular part of uiε, eq 60 Ui ) delta-shock strength of the dimensionless fluid phase concentration of component i Vi ) dimensionless adsorbed phase concentration of component i, Vi ) νKini Viε ) smooth approximation of the adsorption phase delta-shock solution for component i, eq 63 Vjiε ) bounded part of Viε, eq 63 V˜ iε ) singular part of Viε, eq 63 Vij ) partial derivative of Vi with respect to uj V ) delta-shock strength of the dimensionless adsorbed phase concentration, V ) Vi/ai Vi ) delta-shock strength of the dimensionless adsorbed phase concentration of component i V ) matrix defined as V ) [Vij] wi ) hold-up in the delta-shock of component i W ) interstitial fluid phase velocity x ) dimensionless axial coordinate, x ) z/L X ) finite parameter z ) axial coordinate Greek Letters Rj ) finite parameter or bounded function β ) finite parameter eq 52 Γe ) image of the continuous nonsimple wave transition in the hodograph plane Γj ) jth characteristics in the hodograph plane δ(x) ) Dirac’s delta function ∆ ) discriminant of a quadratic equation, eq 13 ε ) small regularization parameter εi ) reciprocal of Peclet number of component i θj ) jth eigenvalue of the matrix V θ˜ j ) equilibrium theory parameter for a Sj shock, eq 27 θ˜ δ ) equilibrium theory parameter for a delta-shock, eq 85 η ) parameter λ˜ j ) nondimensional speed of propagation of a Sj shock λ˜ δ ) nondimensional speed of propagation of a delta-shock ν ) phase ratio ξ ) dimensionless moving coordinate, ξ ) x - st Ξiε ) double integral in equation i, eq 67 σj ) slope of the jth characteristic in the physical plane, eq 6 σ˜ j ) slope of a Sj shock in the physical plane σ˜ δ ) slope of a delta-shock in the physical plane Σj ) image of a Sj shock in the hodograph plane τ ) dimensionless time, τ ) tW/L χ ) parameter ψ ) 1D or 2D test function Ψiε ) double integral in equation i, eq 67 ω ) characteristic parameter, eq 14
7752
Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009
Subscripts and Superscripts 0 ) function centered with respect to the delta-shock; see Figure 13 - ) function on the left-hand side (lhs) of the delta-shock; see Figure 13 + ) function on the right-hand side (rhs) of the delta-shock; see Figure 13 * ) modified isotherm A ) feed state B ) initial state e ) nonsimple wave transition i ) component index I ) intermediate state j ) eigenvalue and eigenvector index L ) state on the lhs of a transition R ) state on the rhs of a transition δ ) delta-shock ε ) regularized function of the small parameter ε
Literature Cited (1) Mazzotti, M.; Tarafder, A.; Gritti, F.; Guiochon, G. Experimental evidence of a delta-shock in nonlinear chromatography. J. Chromatogr. A, 2009, submitted. (2) Rhee, H.-K.; Aris, R.; Amundson, N. R. First order partial differential equations; Prentice-Hall: Englewood Cliffs, NJ, 1986; Vol. I. (3) Rhee, H.-K.; Aris, R.; Amundson, N. R. First order partial differential equations; Prentice-Hall: Englewood Cliffs, NJ, 1989; Vol. II. (4) Mazzotti, M. Local equilibrium theory for the binary chromatography of species subject to a generalized Langmuir isotherm. Ind. Eng. Chem. Res. 2006, 45, 5332. (5) Korchinski, D. J. Solution of a Riemann problem for a 2 × 2 system of conservation laws possessing no classical weak solutions. Ph.D. thesis, Adelphi University, New York, 1977. (6) Keyfitz, B. L.; Kranzer, H. C., A viscous approximation to a system of conservation laws with no classical Riemann solution. In Nonlinear Hyperbolic Problems, Lecture Notes in Mathematics; Carasso, C. Ed.; Springer-Verlag: Berlin/New York, 1989; Vol. 1402, pp 185-197. (7) Tan, D.; Zhang, T.; Zheng, Y. Delta-shock waves as limits of vanishing viscosity for hyperbolic systems of conservation laws. J. Differen. Equat. 1994, 112, 1.
(8) Keyfitz, B. L.; Kranzer, H. C. Spaces of weighted measures for conservation laws with singular shock solutions. J. Differen. Equat. 1995, 118, 420. (9) Nedeljkov, M. Delta and singular delta locus for one-dimensional systems of conservation laws. Math. Meth. Appl. Sci. 2004, 27, 931. (10) Shelkovich, V. M. δ- and δ′-shock wave types of singular solutions of systems of conservation laws and transport and concentration processes. Russ. Math. SurV. 2008, 63, 473. (11) Migliorini, C.; Gentilini, A.; Mazzotti, M.; Morbidelli, M. Design of simulated moving bed units under non-ideal conditions. Ind. Eng. Chem. Res. 1999, 38, 2400. (12) Liu, T.-P. Existence and uniqueness theorems for Riemann problems. Trans. Am. Math. Soc. 1975, 212, 375. (13) Kvaalen, E.; Neel, L.; Tondeur, D. Directions of quasi-static mass and energy transfer between phases in multicomponent open systems. Chem. Eng. Sci. 1985, 40, 1191. (14) Jedrzejak, A.; Gorius, A.; Tondeur, D. Some consistency problems arising in the description of equilibrium adsorption of mixtures. Chem. Eng. Sci. 1989, 44, 1315. (15) Lax, P. Weak solutions of nonlinear hyperbolic equations and their numerical computation. Commun. Pure Appl. Math. 1954, 7, 159. (16) Danilov, V. G.; Shelkovich, V. M. Dynamics of propagation and interaction of δ-shock waves in conservation law systems. J. Different. Equat. 2005, 211, 333. (17) Colombeau, J. F., Multiplication of distributions - A tool in mathematics, numerical engineering and theoretical physics, Lecture Notes in Mathematics; Springer-Verlag: Berlin/Heidelberg, 1992; Vol. 1532. (18) Colombeau, J. F.; Le Roux, A. Y. Multiplications of distributions in elasticity and hydrodynamics. J. Math. Phys. 1987, 29, 315. (19) Colombeau, J. F.; Le Roux, A. Y.; Noussair, A.; Perrot, B. Microscopic profiles of shock waves and ambiguities in multiplications of distributions. SIAM J. Numer. Anal. 1989, 26, 871. (20) Krejic´, N.; Krunic´, T.; Nedeljkov, M. Numerical verification of delta shock waves for pressureless gas dynamics. J. Math. Anal. Appl. 2008, 345, 243. (21) Mazzotti, M. Delta-shock waves in nonlinear chromatography. J. Differen. Equat. 2009, in preparation. (22) Piatkowski, W.; Antos, D.; Gritti, F.; Guiochon, G. Study of the competitive isotherm model and the mass transfer kinetics for a BET binary system. J. Chromatogr. A 2003, 1003, 73.
ReceiVed for reView January 28, 2009 ReVised manuscript receiVed April 26, 2009 Accepted June 24, 2009 IE9001537
