Two-Phase Flow in Liquid Chromatography, Part 2: Modeling

Feb 14, 2018 - Two-Phase Flow in Liquid Chromatography, Part 2: Modeling ... This work considers an equilibrium theory model and a lumped kinetic mode...
1 downloads 8 Views 4MB Size
Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Two-Phase Flow in Liquid Chromatography, Part 2: Modeling Franziska Ortner and Marco Mazzotti* Institute of Process Engineering, ETH Zurich, 8092 Zurich, Switzerland S Supporting Information *

ABSTRACT: This contribution presents the derivation and solution of model equations for a chromatographic system in the presence of multiple fluid phases. While chromatographic models commonly account for adsorption and one single convective phase, macroscopic models for applications in natural reservoirs deal with multiphase flow but usually neglect adsorption effects. This work considers an equilibrium theory model and a lumped kinetic model which account for both multiphase flow and adsorption phenomena. Analytical and numerical solutions of these models are derived and applied to a specific chromatographic system, which was investigated in detail in a companion paper. The analysis of equilibrium theory solutions provides a deep insight in the mathematical and physical implications of two-phase flow and adsorption. The comparison of analytical and numerical solutions validates both models and highlights advantages and drawbacks. The presented model equations and solutions are of considerable interest not only in the context of liquid chromatography but also for several applications involving natural reservoirs.

1. INTRODUCTION Chromatographic model equations are commonly based on component mass balances which account for the adsorption of one or multiple species, as well as for the convection of one single phase.1 Since chromatographic processes are often conducted in a dilute regime, the established models are usually based on the assumption of a dilute, ideal behavior in liquid and adsorbed phases. However, under certain conditions, components can become enriched in the liquid phase, which may result in supersaturation and a subsequent liquid−liquid phase split. Reasons for this enrichment are interactions between different adsorbing components, as observed previously,2,3 but also the impact of modifiers or chemical reactions occurring in chromatographic reactors. Under these conditions, standard chromatographic models fail for the following reasons. First, the concentrated liquid phase(s) feature a highly nonideal behavior, which has to be taken into account when describing the thermodynamic equilibrium between liquid phase(s) and adsorbed phase. Second, in the presence of a phase split and subsequent twophase flow, the flow behavior can be much more complex, as different convective phases can move with different velocities. In this contribution, we set up model equations (based on component mass balances), which account both for adsorption effects and for two-phase flow. This is an equilibrium theory model, assuming thermodynamic equilibrium between all convective phases and between the convective phases and the adsorbed phase, which can be solved analytically with the method of characteristics.4−6 We will also consider a modification of this model (accounting for kinetic limitations in one lumped term) to be solved numerically and compared with the analytical model. © XXXX American Chemical Society

Equilibrium theory models are commonly used in liquid chromatography to understand and to describe the behavior of experimental systems3,5,7−11 and to design chromatographic processes.12,13 While all these models account for adsorption but do not consider multiple convective phases, equilibrium theory models in the context of different applications in natural reservoirs describe two- or multiphase flow through a porous medium but most often neglect adsorption effects.6,14−16 To the best of our knowledge, there exists only one contribution accounting for both multiphase flow and adsorption effects in the context of enhanced coal-bed methane recovery and CO2 sequestration.17,18 Focused on the mathematical derivation of model solutions, such study assumes simple relationships to describe the thermodynamic equilibria between gas and liquid phase (constant distribution coefficients) and between convective and solid phases (Langmuir isotherm as a function of the partial pressures in the gas phase). Accounting for the physically realistic thermodynamic properties of a specific experimental system (described in the companion paper19) and combining adsorption and two-phase flow, the model presented in the current contribution constitutes a combination and extension of equilibrium theory models used in liquid chromatography and in applications to natural reservoirs. The paper is structured as follows: In section 2, model assumptions are specified and the equilibrium theory model, as well as a lumped kinetic model, is established. Subsequently, the model system (investigated in detail in the companion Received: Revised: Accepted: Published: A

December 14, 2017 February 13, 2018 February 14, 2018 February 14, 2018 DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research paper19 is introduced in section 3, providing information about thermodynamic liquid−liquid and adsorption behavior and about the two-phase flow behavior. In sections 4 and 5, the equilibrium theory equations are solved, assuming equal or different interstitial velocities of the convective phases, respectively. A general solution of the model equations is provided, which is subsequently applied to the model system, and characteristics in the hodograph plane, as well as elution profiles for specific initial and feed states are derived. On the basis of these results, mathematical and physical implications of the equilibrium theory model are highlighted. Finally, the lumped kinetic model is solved numerically and simulations are compared to the equilibrium theory results in section 6. Conclusions are drawn in section 7.

is described with thermodynamic consistency via an adsorption isotherm, which gives the adsorbed phase concentrations, ni, as a function of the liquid phase activities, ai, that are identical for all liquid phases in equilibrium. Dealing with a nondilute system, volumetric changes due to adsorption/desorption cannot be neglected. Thus, we define the stationary phase to be made up of the adsorbent (with reference volume (1 − ϵref)Vc and reference porosity ϵref at the completely regenerated state, i.e., in the absence of any adsorbed component) and the adsorbed phase, and we assume every component to occupy the same volume in the adsorbed as in the liquid phase (defined by the component densities ρi), thus yielding the following mathematical relationship: 1 − ϵref 1−ϵ

2. MODEL EQUATIONS 2.1. Equilibrium Theory Model. The component mass balances for a multicomponent system in the presence of multiple convective (liquid) phases, accounting for adsorption and neglecting dispersive effects and mass transfer resistances, is given as ∂(ϵCi) ∂n ∂(uFi ) + (1 − ϵref ) i + = 0, ∂t ∂t ∂z

NP

NP

NC i=1

j=1

j=1

∑ xijρj̃ f j j=1

(5)

(6)

Thus, under the previous assumptions of variable porosity, volume additivity in the liquid and the adsorbed phase, and equal component densities in the liquid as in the adsorbed phase, the superficial velocity (overall flow rate divided by the cross-sectional area of the column) remains constant. Finally, time and space coordinates are made dimensionless:

∑ cijf j j=1

ξ=

(2)

z ; L

τ=

ut L

(7)

where L is the column length. The final form of the component mass balances is then

where NP denotes the number of phases, xij and cij are the mass fraction and concentration of component i in phase j, respectively, and ρ̃j is the density of phase j (note that ρi corresponds to the density of pure component i). The volumetric fraction of phase j, also called phase saturation, is indicated by Sj, whereas f j is the volumetric fractional flow of phase j. Accounting for the fact that different convective phases can move with different interstitial velocities, the fractional flow f j is not necessarily equal to Sj. Commonly, f j is expressed as an empirical (rational) function of Sj, as discussed in detail below. As a consequence of f j being in general not equal to Sj, the values of Ci and Fi can differ, since they are averages of the phase concentrations cij, weighed by the phase saturations Sj or the fractional flow f j, respectively. The phase densities ρ̃j are given as follows, as a function of the phase composition and of the component densities, by assuming additivity of volumes: ⎛ NC xij ⎞ ρj̃ = ⎜⎜∑ ⎟⎟ ⎝ i = 1 ρi ⎠

ni ρi

∂u =0 ∂z

NP

=

(4)

By normalizing each component mass balance by the component density ρi and summing up over all the component mass balances, the overall mass balance reduces to

i = 1, ..., NC

NP

Fi =

i=1

ni ⎞ ⎟=1 ρi ⎟⎠

ϵ = ϵref − (1 − ϵref ) ∑

where t and z are time and space variables, respectively, NC is the number of components, and ni is the adsorbed phase concentration of component i in thermodynamic equilibrium with the liquid phase(s) (mass of component i per volume of adsorbent in the completely regenerated state). The overall liquid phase concentrations, Ci, and overall fractional flows, Fi, are defined as

∑ xijρj̃ Sj = ∑ cijSj ,

NC



which can be recast as

(1)

Ci =

⎛ ⎜⎜1 + ⎝

∂(ϵCi) ∂n ∂F + (1 − ϵref ) i + i = 0, ∂τ ∂τ ∂ξ i = 1, ..., (NC − 1)

(8)

The NC − 1 independent concentrations Ci are uniquely determined by the system of NC − 1 component mass balances, while the concentration of component NC is a function of the other component concentrations, based on the assumption of volume additivity: NC

∑ i=1

Ci = ρi

NC

∑ i=1

cij ρi

=1 (9)

Since the conservation eqs 8 form a system of first-order, homogeneous, partial differential equations with unknowns Ci (i = 1, ..., (NC − 1)), the model can be solved by applying the method of characteristics. 2.2. Lumped Kinetic Model. We further consider a modification to the above conservation laws, to be solved numerically, which lumps dispersion and mass transfer effects into a simple empirical mass transfer term. It should be pointed out that, while this additional term accounts for possible dispersive effects and mass transfer resistances between liquid

−1

(3)

We further assume thermodynamic equilibrium between the liquid phases (i.e., equal chemical potentials of each component in all liquid phases), which allows us to calculate phase split and phase compositions based on a thermodynamic model. The equilibrium between the adsorbed phase and the liquid phases B

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Table 1. UNIQUAC Parameters for the Exemplary Systema

a

components

1−2

2−1

1−3

3−1

2−3

3−2

Δuij/R [K]

1.097 × 103

−76.94

1.165 × 103

270.7

−83.02

201.0

Adopted from the companion paper.19

3. MODEL SYSTEM Let us consider a ternary system with one adsorbing component (component 1) and two inert components (components 2 and 3) at a constant temperature of 296.15 K. System properties are those of the experimental system studied in the companion paper19 with phenetole, methanol, and water being components 1−3 and a reference porosity of the packed column of ϵref = 0.615. Note that for this system with only one adsorbing component eq 5 reduces to n ϵ = ϵref − (1 − ϵref ) 1 ρ1 (15)

phases and adsorbed phase, the assumption of thermodynamic equilibrium between the liquid phases is nevertheless maintained. The resulting modified conservation equations read: ∂F ∂(ϵCi) ∂n ̃ + (1 − ϵref ) i + i = 0, ∂τ ∂ξ ∂τ i = 1, ..., (NC − 1)

∂nĩ = St i(ni − nĩ ), ∂τ

(10)

i = 1, ..., NC

(11)

3.1. Thermodynamic Equilibria. The nonideal features of the liquid phase(s) can be described by the UNIQUAC model in the companion paper19 and ref 20, with the parameters given in Table 1. This model allows for the determination of liquid phase activities for any composition in the ternary diagram, and, for a composition located in the immiscible region, it allows calculating the phase split (phase saturations Sj) and the compositions xij of the two liquid phases in thermodynamic equilibrium. The binodal curve (black line), the plate point (empty circle), and a few exemplary tie-lines are illustrated in Figure 1. (For the calculation procedure, see the companion

The variable ni denotes the adsorbed phase concentration in equilibrium with the liquid phase(s) throughout the entire work, whereas its actual, nonequilibrium value, in the context of eq 11 is called ñi. The dimensionless Stanton number Sti is defined as St i =

akiL u

(12)

where a is the specific area of the adsorbent and ki is the mass transfer coefficient of component i. The introduction of the mass transfer term results in band broadening effects in the elution profiles, reducing the sharpness of fronts, which is physically more realistic than the very sharp transitions obtained under the assumption of instantaneous thermodynamic equilibrium. However, this lumped kinetic term was not only introduced for physical but mainly for numerical reasons, as it circumvents the calculation of derivatives of ni with respect to C, which would have been necessary for a numerical solution of the equilibrium theory model. With n1 being a function of the liquid phase activities, analytical expressions of such derivatives would have been rather complex in the miscible region (involving a differentiation of the UNIQUAC equation) and impossible in the immiscible region. The introduction of the mass transfer term instead of a numerical calculation of the derivatives of n1 allowed a reduction in computational effort. A rather high Stanton number, Sti = 200 (corresponding to minor dispersion effects), was assumed for all components, so the solution of the lumped kinetic model is expected not to deviate excessively from that obtained using the equilibrium theory model. Since in this work we are dealing with Riemann problems (i.e., stepwise constant initial value problems), the initial condition is defined as Ci(ξ , 0) = Ci0 ;

ni(ξ , 0) = ni0

Figure 1. Ternary diagram (mass fractions) of the model system (phenetole (1), methanol (2), and water (3)), with binodal curve (continuous black line), plait point (open circle) and a selection of tielines (thin gray lines) predicted by the established UNIQUAC model.

paper.19 To avoid computationally expensive equilibrium calculations and thus to increase the simulation speed, 5000 calculated tie-lines spanning the entire two-phase region were used as a look-up table for solving both the equilibrium theory model and the lumped kinetic model. The implementation of the look-up table is model-specific and is explained in detail in the Supporting Information. The adsorption equilibrium of species 1 is described by an anti-Langmuir isotherm, defined as a function of the liquid phase activity (which is obviously equal for all the liquid phases in thermodynamic equilibrium), a1:

(13)

where C0i and n0i define the initial state of the column at τ = 0. For the first-order lumped kinetic model used in this work, a Dirichlet boundary condition can be applied:

Fi(0, τ ) = FiF

(14)

where FFi denotes the feed overall fractional flow. C

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

n1 =

H1a1 1 − K1a1

Table 2. Parameters of the Relative Permeability Functions for the Exemplary Systema

(16)

The adsorption isotherm parameters are H1 = 41.5 g L−1 and K1 = 0.862 (from the companion paper19). Obviously, n2 = n3 = 0. 3.2. Two-Phase Flow Behavior. For single-phase flow in a porous medium, Darcy’s law establishes that the superficial velocity u depends linearly on the pressure gradient dP/dξ (assuming linear, one-dimensional flow while neglecting gravitational effects).21 The constant of proportionality is the permeability K, an intrinsic property of the porous medium, divided by the dynamic viscosity μ of the fluid phase. In the presence of multiple liquid phases, the extended Darcy’s law relates the superficial velocity of phase j, uj , to the corresponding pressure gradient dP/dξ through the relative permeability krj of the phase j itself, as follows:22

uj = f j u = −

parameters values a

(17)

The relative permeability, ranging between 0 and 1, accounts for the fact that the phase j cannot occupy the entire void space within the porous medium, but has to share that space with the other fluid phases present. Thus, the superficial velocity of that phase at the same pressure gradient is smaller than it would be in the case of single phase flow. Assuming the same pressure for all fluid phases, Pj = P, at any position in the column, i.e., neglecting the effect of the capillary pressures due to curved interfaces between phases, one can derive the following relationship for the fractional flow f j: k rj

fj =

N

k ri μi

(18)

Under the assumption of equal velocities of all the convective phases, it is obvious that f j = Sj. However, it has been observed that the relative permeabilities of many experimental systems can be accurately described by a power law relationship, given as a function of the phase saturations S.23 For our system containing a maximum of two liquid phases, we assume that the following relationships hold true (superscripts L and R denote the component 1 lean and rich phase, respectively): krL = 0; krR = 1

L S L ≤ S lim

(19a)

L 2 krL = krL,max(Seff ) ; L 3 L R ≤ S L ≤ (1 − S lim krR = krR,max(1 − Seff ) S lim ) (19b)

krL = 1; krR = 0

R (1 − S lim ) ≤ SL

(19c)

L Seff =

S −

Adopted from the companion paper.19

⎛ ⎛ ∂C1 C ⎞⎞ ∂C + ⎜⎜(1 − ϵref )n12⎜⎜1 − 1 ⎟⎟⎟⎟ 2 + =0 ∂ ∂ξ ρ τ ⎝ ⎝ 1 ⎠⎠

L S lim

L R − S lim 1 − S lim

SlimR 0.32

⎛ ⎛ C ⎞⎞ ∂C ⎜ϵ + (1 − ϵref )n11⎜⎜1 − 1 ⎟⎟⎟ 1 ⎜ ρ1 ⎠⎟⎠ ∂τ ⎝ ⎝

with L

SlimL 0.00

4. SOLUTION OF THE EQUILIBRIUM THEORY MODEL ASSUMING EQUAL INTERSTITIAL VELOCITIES In a first step, we analyze the dynamic behavior of a chromatographic system where the two liquid phases move at the same interstitial velocity, i.e., when assuming f j = Sj and thus Fi = Ci in the model equations above. We first provide the general solution of the equilibrium theory model for equal velocities, using the method of characteristics.5 Then, we apply the general solution to the model system of interest here. We calculate the network of characteristics in the hodograph plane, i.e., the (C2, C1) plane, which provides images of the solutions independent of initial and feed states, and we explore characteristic features and physical implications of the solutions. Finally, the chromatographic behavior is illustrated by elution profiles for specific initial and feed states. 4.1. General Solution. Under the assumption that f j = Sj and thus Fi = Ci and keeping in mind that n2 = n3 = 0, eq 8 for components 1 and 2 can be recast as

μj

∑i =P1

kr,maxR 0.77

pressure drop, in the presence of the other phase Y being completely trapped at a saturation SYlim. The assumed parameters are independent of the interfacial tension between the phases, i.e., of the position of the tie-lines with respect to the plait point. In principle, as the properties of the two phases in equilibrium approach each other when approaching the plait point, the interfacial tension between the phases decreases, hence the two-phase flow behavior should approach a single-phase flow behavior ( f j = Sj). For the sake of simplicity of the equilibrium theory derivations, however, such dependence is neglected in this work. To fully determine the fractional flow functions (eq 18), we assume a constant viscosity μR = 1.17 mPas for the component 1 rich phase. The viscosity of the component 1 lean phase is estimated for the relevant solvent ratio cL2 /cL3 , assuming that the ternary mixture behaves like the binary mixture methanol (component 2)−water (component 3); for details, see the companion paper.19 An estimation procedure proposed in the literature was applied,24 with parameters determined from viscosity data provided elsewhere.25 The assumed relative permeability functions and the resulting relationship f X−SL, assuming viscosities μR = 1.17 mPas and μL = 1.42 mPas, are illustrated in Figure 2.

k rjK dPj μj dξ

kr,maxL 0.94

(21a)

⎛ ⎞ ∂C ⎛C ⎞ ∂C ∂C2 C ⎜⎜ϵ − 2 (1 − ϵref )n12⎟⎟ 2 − ⎜⎜ 2 (1 − ϵref )n11⎟⎟ 1 + (21b) ∂ξ ρ1 ⎝ ⎠ ∂τ ⎝ ρ1 ⎠ ∂τ =0

(20)

The model parameters in the relative permeability functions are reported in Table 2, as obtained for the experimental system investigated in the companion paper.19 The parameters SXlim indicate the volume fraction of phase X which is completely trapped in the porous medium. Maximum relative permeabilities kXr,max determine the flow velocity of phase X at a certain

Here, C1 and C2 are the unknown variables, and n1j denotes ∂n1/∂Cj (j = 1, 2). Note that assuming an isotherm depending on the liquid phase activity a1, and not directly on the liquid D

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

The Γi characteristics, resulting from the solution of eqs 22 and 23, describe the change in composition (C1, C2) that occurs when connecting an initial state to a feed state. Equation 22 can be written as ∂n1 ∂n1 dC1 + dC2 = 0 = dn1 ∂C1 ∂C2

(24)

which allows to conclude that the adsorbed phase concentration n1 remains constant along the Γ1 characteristics. Substituting eq 9 both in finite and in differential form into eq 23 yields the following differential equation in C2 and C3

dC 3 C = 3 dC 2 C2

(25)

Such equation holds along Γ2 and implies that the ratio of the concentrations of the two inert components (in the following called solvent ratio) remains constant along the Γ2 characteristics. Each composition on the solution path mapping onto Γ1 (Γ2) characteristics moves through the column at a constant propagation velocity λ1 (λ2). The reciprocals of the propagation velocities, σj = λ−1 j , determine the slopes of the characteristics in the physical plane (ξ, τ); for the investigated system of eqs 8, the slopes are σ1 = ϵ (26) ⎛ C ⎞ C σ2 = ϵ + (1 − ϵref )n11⎜⎜1 − 1 ⎟⎟ − (1 − ϵref )n12 2 ρ ρ1 ⎝ 1⎠

Figure 2. (a) Relative permeability function assumed for the exemplary system. (b) Fractional flow function based on the relative permeability functions presented in (a) and assuming μR = 1.17 mPas and μL = 1.42 mPas.

(27)

Since both n11 and n12 can be positive or negative, it is difficult to make a general statement about whether the value of σ2 is larger than ϵ or not; however, it should be pointed out that this is always the case for our system (thus, σ2 > σ1). We consider a Riemann problem with an initial (downstream) state A (CAi for ξ ≥ 0 and τ = 0) and a feed (upstream) state B (CBi for ξ = 0 and τ > 0). The solution in the physical plane (ξ,τ) consists of three states, namely, the initial state A, an intermediate state I, and the feed state B. These three states are connected by two transitions, each of which maps onto a Γj (j = 1, 2) characteristic. Depending on the directional derivative of σj along the Γj characteristic, the transition can be a contact discontinuity (directional derivative equals zero), a simple wave (positive directional derivative moving from downstream states D to upstream states U), a shock (negative directional derivative moving upstream), or a semishock (directional derivative changing sign). Since σ1 remains constant along Γ1, the directional derivative is 0, and consequently, a transition mapping on a Γ1 characteristic is always a contact discontinuity. In turn, the directional derivative of σ2 along Γ2 can be positive or negative, allowing for different types of transition. If σ2 increases when moving from downstream states D to upstream states U, then the transition is a simple wave; in the opposite case, it is a selfsharpening shock, i.e., a “weak” solution obtained from the integral form of the mass balance and with a slope in the physical plane σ̃2 given by

phase concentration C1, the dependence of n1 on Cj (j = 1, 2) is not trivial. The isotherm equation itself gives n1 as a function of a1, which is a function of the mole fractions of the three components, i.e., a1 = a1(Z1, Z2, Z3), through the UNIQUAC equation; mole fractions Zi in turn depend on the liquid phase concentrations, i.e., Zi = Zi(C1, C2, C3), as specified in the Supporting Information. Finally, based on the assumption of volume additivity, the number of unknown variables is reduced to 2, by expressing C3 as a function of C1 and C2 using eq 9 (for the details of such calculation, see the Supporting Information). Applying the method of characteristics,5 one can derive equations for the characteristics, i.e., the images of the solutions, in the hodograph (C2, C1)-plane and in the physical (ξ, τ)-plane, based on which it is then possible to determine elution profiles for specific initial and feed conditions. The procedure of deriving these analytical solutions has been explained in great detail in the literature.3,5,8,10,11 For the sake of brevity, we keep explanations concerning the solution procedure to a minimum and focus on the specific results obtained with this model. In the case of hyperbolic systems of two equations, such as that given by eqs 21, there are two sets of characteristics, Γ1 and Γ2, with slopes ζ1,2 = dC2/dC1, in the hodograph plane given by n Γ1: ζ1 = − 11 n12 (22)

C2 Γ2: ζ2 = − ρ1 − C1

σ̃2 =

[ϵC1] + [(1 − ϵref )n1] [ϵC2] = [C1] [C2]

(28)

The sign [·] denotes a jump of the enclosed quantity across the discontinuity. Note that the conditions given in eq 28 for

(23) E

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research components 1 and 2, which both have to be fulfilled for a shock transition, also provide the possible states which, assuming a certain downstream state D (a certain upstream state U), can be connected to D (U) through a shock. These states are located on a shock path Σ2, which is tangent to the Γ2 characteristic in D (U). Since in our case Γ2 characteristics are straight lines, a shock path Σ2 coincides with a characteristic Γ2. A further requirement for a shock connecting a downstream state D to an upstream state U is that σU2 < σ̃2 < σD2 .4 If the directional derivative of σ2 changes in sign along the Γ2 characteristic connecting the intermediate to the feed state, then the transition will be a semishock (shock-wave or waveshock). A detailed explanation of the necessary conditions and the behavior of semishock transitions is provided in the literature.3−5 The two transitions in the solution of the Riemann problem map on two different types of characteristics Γ1 and Γ2, which intersect at the intermediate state I and connect I to initial state A and feed state B. In a physically meaningful solution, downstream states have to propagate through the column faster than upstream states. Thus, with σ1 < σ2 and hence λ1 > λ2, the intermediate state has to be connected to the initial (downstream) state A by a Γ1 characteristic, and to the feed (upstream) state B by a Γ2 characteristic. A physical interpretation of the derived mathematical solution is as follows. Along a Γ1 characteristic, the adsorbed phase concentration n1 remains constant. Since neither adsorption nor desorption occur along such a characteristic, all states along it move at the same propagation velocity, namely, the interstitial velocity ϵ−1. In turn, along a Γ2 characteristic, the concentration ratio of inert components 2 and 3 remains constant. Since neither inert component interacts with the adsorbent, and all convective phases move with the same interstitial velocity, these components propagate at the same velocity. Therefore, unless the concentration ratio of inert components is different in the initial and in the feed state, i.e., these two states are not located on the same Γ2 characteristic, there is no driving force for the ratio C2/C3 to change. As the concentration and activity of component 1 vary along these characteristics, adsorption or desorption occur and propagation velocities vary as well, leading to the formation of waves, shocks, or semishocks. 4.2. Hodograph Plane. In the following, the equations derived above are applied to determine characteristics in the hodograph plane for the model system introduced in section 3, while still assuming equal velocities of the convective phases ( f j = Sj). Γ1 and Γ2 characteristics in the hodograph plane are plotted in Figure 3a,b, respectively. Γ1 characteristics were calculated by solving the ordinary differential eq 22 with a builtin Matlab routine (ode113). In turn, the calculation of Γ2 characteristics is trivial, as they feature a constant solvent ratio C2/C3 and thus are straight lines in the ternary diagram. With n1 being constant along blue Γ1 characteristics (see eq 22), based on eq 16, Γ1 characteristics constitute paths of constant a1. The qualitative behavior of the Γ1 characteristics is illustrated in Figure 4, where the miscible region is enlarged for better visibility. All Γ1 characteristics connect a point on the x2 = 0 axis with a point on the x3 = 0 axis. In the soluble region, these characteristics are arc-shaped, while in the immiscible region they map on the tie-lines. Characteristics emanating from a point M located on the x2 = 0 axis and close to the x3 = 1 corner remain entirely in the miscible region. Increasing the mass fraction of component 1 (xM 1 ) of this starting point, a Γ1

Figure 3. Characteristics in the Hodograph plane (mass fractions) for the ternary system with one adsorbing component introduced in section 3, assuming equal interstitial velocities of all convective phases. (a) Γ1 characteristics and (b) Γ2 characteristics. The black curve corresponds to the binodal curve. Red arrows indicate the direction in which σ2 increases along a Γ2 characteristic. Red squares denote conditions at which the directional derivative of σ2 changes from positive to negative.

characteristic, located entirely in the miscible region and tangent to the plait point, is attained. A further increase of xM 1 results in Γ1 characteristics passing through the immiscible region along tie-lines and extending into the miscible region to reach the x2 = 0 and x3 = 0 axes. Once xM 1 reaches the solubility limit of component 1 in component 3, the corresponding Γ1 characteristic maps onto the tie-line located on the x2 = 0 axis. Finally, if M is located in the miscible region on the component 1 rich side, then the corresponding Γ1 characteristic again connects directly to the x3 = 0 axis, without passing through the immiscible region. It should be kept in mind that along Γ1 characteristics, σ1 remains constant. Thus, the directional derivative equals 0, and transitions along these paths are always contact discontinuities. As illustrated in Figure 3b, Γ2 characteristics are straight lines along which the solvent ratio C2/C3 remains constant, i.e., connecting a point on the x1 = 0 axis to the x1 = 1 corner. Red arrows in Figure 3 indicate the direction in which σ2 increases, while red squares denote locations at which the directional F

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 4. Illustration of the qualitative behavior of Γ1 characteristics in the hodograph plane.

Figure 5. Paths in the hodograph plane (calculated based on eqs 22 and 23) for exemplary cases with an initial state A and two different feed states in the soluble region. Paths describe the solution for the entire chromatographic cycle, i.e., adsorption and desorption. Feed state B1: all paths are located in the soluble region; feed state B2: paths lead through the immiscible region during desorption. States Ia and Id denote intermediate states formed during adsorption and desorption, respectively. Note that blue and red paths correspond to Γ1 and Γ2 characteristics, as they are illustrated in Figure 3a,b, respectively.

derivative of σ2 equals zero and changes sign. As discussed above, the nonmonotonic behavior of σ2 allows the formation of wave, shock, and semishock transitions. Note that at the intersection of Γ2 characteristics with the binodal curve, σ2 is discontinuous, so concentration plateaus can form at these locations in the corresponding elution profiles, although the solution path remains on the same characteristic. 4.3. Solution of Exemplary Cases. The behavior of the exemplary system shall now be illustrated for two different cases. As the equilibrium theory solution under the current assumption of equal interstitial velocities of the convective phases is rather simple, the examples that follow can be easily exploited to analyze solutions for different initial and feed states (also for those located in the immiscible region). With reference to Figure 5, the initial state A is located on the x1 = 0 axis, i.e., a mixture of the inert components 2 and 3. Let us consider two feed states, B1 and B2, which are located in the soluble region on the x3 = 0 axis, i.e., two different binary mixtures of components 1 and 2. In both cases, we consider a complete chromatographic cycle, consisting of the adsorption step displacing the initial state A by the feed state Bi, and of the desorption step eluting the feed state by the initial state. Accordingly, for the derivation of the solution of the desorption step, the initial state and the feed state are swapped. The corresponding paths in the hodograph plane are shown in Figure 5, whereas elution profiles are illustrated in Figure 6 with adsorption profiles on the left and desorption profiles on the right. Black thin lines in the concentration profiles show the solution obtained by solving numerically the lumped kinetic model of eqs 10 and 11 (see the discussion in section 6 below). As described in the general solution above, all paths connecting initial to feed states emanate from the initial state along a Γ1 characteristic until reaching an intermediate state I, which is then connected to the feed state by a Γ2 characteristic. As a consequence, the first transition connecting the initial state to I is always a contact discontinuity propagating with the interstitial velocity, i.e., appearing at the column outlet at the void time. The second transition in most cases (and in our two examples) is a shock, although under certain conditions it can also be a wave or a semishock. Whenever an initial state located on the x1 = 0 axis is connected to a feed state located on the x3 = 0 axis, the intermediate state of the adsorption step Ia is located in the x2 =

1 corner. The physical reason for this is that the inert component 3 (absent in the feed state) propagates at the interstitial velocity and is thus eluted within the void time. In turn, the adsorbing component 1 (absent in the initial state) propagates more slowly and arrives at the column outlet later; hence, at the intermediate state, only component 2 is present. In general, for all initial states located on the x1 = 0 axis connected to any feed state in the miscible region, the image of the adsorption step belongs entirely to the miscible region. The desorption step for initial and feed states located in the miscible region can pass through the two-phase region, as seen for the example with feed state B2, where the two-phase region is crossed along a Γ1 characteristic. The intermediate states of the desorption steps Id of the two exemplary cases are both located in the miscible region and on the same Γ2 characteristic (determined by state A) but differ due to the different Γ1 characteristics to be considered (determined by state Bi). The resulting elution profiles (see Figure 6, right column) are qualitatively similar, with both intermediate states located in the miscible region, connected to states Bi through a contact discontinuity and to A through a shock. A potential crossing of the two-phase region (as is the case for example 2) occurs during the contact discontinuity. We want to point out that an analogue solution can be derived for most initial and feed states located in the soluble region. An intermediate state located in the two-phase region can occur under specific and rare initial and feed states, such as in the case that during desorption, the component 1 lean phase of the tie-line involved in the solution path has a solvent ratio C2/C3, which is lower than the solvent ratio of state A. This is only the case for feed states B located very close to the x1 = 1 corner, and is not considered any further in this work. G

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

5.1. General Solution. Solving the equilibrium theory model with different interstitial velocities based on the dependent variables C1 and C2 is possible but results in rather complex expressions for the characteristics in the hodograph and the physical plane. Instead, we perform the variable transformation suggested in ref 6, by defining two new variables, namely, the saturation of the phase rich in component 1, SR, and a tie-line indicator η, uniquely defining the tie-line on which a state is located. Here, we choose η to be the logarithm of the mass fraction of component 1 in the lean phase η = ln(xL1 ). The original system of eqs 8 in the unknowns C1 and C2 is transformed into the following system in the unknowns SR and η: Ai

∂η ∂η ∂SR ∂SR + Bi + Ci + Di = 0, ∂ξ ∂τ ∂ξ ∂τ

i = 1, 2

(29)

with Ai = f R

dciR dc L ∂f R + (1 − f R ) i + (ciR − ciL) , dη dη ∂η

i = 1, 2 (30a)

⎛ ⎛ dc R dc L ⎞ C ⎞ dn B1 = ϵ⎜SR 1 + (1 − SR ) 1 ⎟ + (1 − ϵref )⎜⎜1 − 1 ⎟⎟ 1 dη ⎠ ρ1 ⎠ dη ⎝ dη ⎝ (30b)

⎛ dc R dc L ⎞ C dn B2 = ϵ⎜SR 2 + (1 − SR ) 2 ⎟ − 2 (1 − ϵref ) 1 d η d η ρ dη ⎝ ⎠ 1 Ci = (ciR − ciL)

∂f R ∂SR

Di = ϵ(ciR − ciL),

,

i = 1,2

(30c)

(30d)

i = 1,2

(30e)

Note that with n1 being a function of the liquid phase activity of component 1 (a1) and ϵ being defined as in eq 5, both ni and ϵ remain constant along the same tie-line; thus, ∂n1/∂SR = ∂ϵ/∂SR = 0. This variable transformation results in simpler expressions for the characteristics, allowing to draw interesting physical conclusions, as seen below. We now define two sets of characteristics Γt and Γnt in the new hodograph plane with coordinates SR and η, with slopes ψt/nt = dη/dSR. These characteristics are the image of the solution in the hodograph plane. Substituting ψ in eqs 29 yields two equations in the unknown SR, which must be identical. This is the case if and only if ψ is a solution of the following quadratic equation:5

Figure 6. Concentration profiles of the chromatographic cycles connecting the initial state A to the feed state (a) B1 and (b) B2. Left: adsorption, right: desorption. Colored thicker lines correspond to equilibrium theory solutions, thin black lines to solutions of the lumped kinetic model solved with a number of discretization steps N = 1000. The corresponding images of the solution in the hodograph plane are shown in Figure 5.

5. SOLUTION OF THE EQUILIBRIUM THEORY MODEL ASSUMING DIFFERENT INTERSTITIAL VELOCITIES In this section, we derive the equilibrium theory solution, considering the more realistic case of different interstitial velocities of the convective phases, i.e., f j = f j(Sj) ≠ Sj and Fi ≠ Ci. The solution described in this section applies only within the immiscible region, while in the single-phase region, the solution derived in section 4 remains valid. We proceed as in the previous section, by first providing the general solution for the equilibrium theory model, before applying it to the model system. Here the fractional flow function introduced in eq 18 is utilized. Again, we investigate the hodograph plane, as well as elution profiles derived for specific initial and feed conditions. Note that in the main text we focus on the most important mathematical and physical attributes of the equilibrium theory solution. A detailed explanation of the construction of solutions is provided in the Supporting Information.

(A1B2 − A 2 B1)ψ 2 + (B2 C1 − B1C2 + A1D2 − A 2 D1)ψ + C1D2 − C2D1 = 0

(31)

Since C1D2 = C2D1, the following two expressions are the two exact solutions of the last equation, which in fact define the two families of characteristics:

Γt: ψt = 0 Γnt: ψnt =

(32)

B2 C1 − B1C2 + A1D2 − A 2 D1 A 2 B1 − A1B2

(33)

From eq 32 and the definition of η, it is clear that Γt characteristics map on the tie-lines. We therefore refer to Γt characteristics as tie-line characteristics and to Γnt characteristics as non-tie-line characteristics (hence the label used as subscript). H

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research The propagation velocities in the physical plane λx = dξ/dτ (x = t, nt), associated to states on the corresponding characteristics Γx, are derived by the method of characteristics as λt =

R C 1 1 ∂f = 1 = σt D1 ϵ ∂S R

(34)

λnt =

A1ψnt + C1 1 = σnt B1ψnt + D1

(35)

The directional derivatives of both λt and λnt can be negative or positive; thus, no conclusions about the type of transition can be drawn a priori. In the case where the directional derivative of λx increases along a Γx characteristic connecting a downstream state to an upstream state, a “weak” solution is required, and the propagation velocity of the corresponding shock transition can be derived from the integral mass balance as [ϵC1 + (1 − ϵref )n1] [ϵC2] 1 = σ̃ = = [F1] [F2] λ̃

(36)

The Rankine−Hugoniot condition given by eq 36 also determines the shock paths Σx in the hodograph plane, since all conditions being connected to a specified initial condition through a shock are uniquely identified by the two expressions of σ̃. It is worth recalling the existence of a further condition to be fulfilled by a shock transition connecting a downstream state D to an upstream state U, namely, σUx < σ̃ < σDx (see section 4.1). As discussed below, states in the immiscible region are most often connected to states in the miscible region through shocks. Since eq 36 is derived from the integral component mass balances, it remains valid also for these shocks crossing the binodal, whereas for the single phase state it is obvious that Ci = Fi = ci. 5.2. Hodograph Plane. We now apply the general solution to the model system to investigate the behavior of the characteristics in the hodograph plane and of the propagation velocities associated to the corresponding liquid phase compositions. Figure 7a presents the Γt and Γnt characteristics in blue and red, respectively, in the (ln(xL1 ), SR)-plane. As expected, Γt characteristics map on the tie-lines and are thus straight horizontal (blue) lines in Figure 7a. Each Γ t characteristic is tangent to two Γnt characteristics; the points of tangency where they meet are singular points where the discriminant of the quadratic eq 31 is zero. Hence, ψt = ψnt and λt = λnt. These points are called watershed points (or, in some contributions, equal eigenvalue points). For the discussed system, there is an infinite number of watershed points, which are all located in the physically relevant domain of (η, SR) values. Let us now consider the propagation velocities λt and λnt along the characteristics, which are illustrated for two selected tie-line characteristics in Figure 8a and for four selected non-tieline characteristics in Figure 8b. The selected characteristics are the highest and the lowest tie-line characteristic in Figure 7a, as well as the two pairs of non-tie-line characteristics with watershed points located on the highest and lowest tie-line characteristic. Considering λt and λnt along their corresponding characteristics (Γt and Γnt, respectively), it can be observed that the propagation velocities do not exhibit a strictly monotonic behavior. The velocity λt exhibits a maximum at a saturation SR

Figure 7. Characteristics in the hodograph plane, constructed based on the transformed variables SR and η = ln(xL1 ) and assuming the hydrodynamic behavior described by eqs 18. (a) Tie-line characteristics (Γt, blue) and non-tie-line characteristics (Γnt, red), as defined in eqs 32 and 33. (b) Γ1 and Γ2 characteristics, constructed by combining parts of Γt and Γnt such that λ1 > λ2 always. Circles (in the color of the characteristics) indicate a change in the directional derivative of σx along a Γx characteristic (x = t, nt, 1, 2) from negative to positive, squares reflect a change from positive to negative.

ranging between 0.6 and 0.8 and decreases to 0 when approaching 0 and 1. The slopes of the propagation velocities λnt along non-tie-line characteristics exhibit multiple changes in sign. With directional derivatives of λx along Γx characteristics changing in sign, complex semishock transitions, which possibly exhibit multiple wave and shock parts, can be formed. Points where the directional derivatives are zero (and change sign) are indicated in Figure 7a as open circles and squares. With reference to Figure 8, there are two points along each tie-line characteristic (one point along each non-tie-line characteristic) where λt = λnt, which are indicated in Figure 8 as open black circles and correspond to the watershed points. I

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

A final observation concerning Figure 8 is that unlike in the case of equal interstitial velocities treated in section 4 neither of the two propagation velocities λt and λnt is always larger than the other. Focusing on the behavior along tie-line characteristics (Figure 8a), it can be noted that between the two watershed points λt > λnt. In turn, beyond the watershed points toward SR → 0 or SR → 1, the propagation velocities λnt are larger than λt. This behavior complicates the derivation of solution paths and elution profiles for specific initial and feed states. As already discussed for the case of equal interstitial velocities (section 4.1) in a physically meaningful solution, downstream states have to propagate faster than upstream states. Hence, a physically meaningful intermediate state has to be reached from upstream states with a propagation velocity higher than the propagation velocity with which it is connected to downstream states. For the solution discussed in section 4.1, with σ1 < σ2 and hence λ1 > λ2, this consideration leads to the simple rule that the solution path would always map on a sequence of Γ1 characteristic (connecting a downstream state A to the intermediate state I) and Γ2 characteristic (connecting the intermediate state I to an upstream state B), intersecting at I. In the case of different interstitial velocities discussed in this section, with no defined order of λt and λnt, the derivation of solution paths is therefore more complex. In order to simplify the derivation of solution paths and elution profiles for specific initial and feed states, and following the suggestion of Helfferich,14 we define two new sets of paths, namely, “fast” paths and “slow” paths, with propagation velocities along the “fast” paths being always larger (or equal) to the propagation velocities along the “slow” paths. With reference to the discussed behavior of Γt and Γnt characteristics and their propagation velocities, this corresponds to reconnecting segments of Γt and Γnt characteristics at the watershed points in a different manner. The resulting “fast” and “slow” paths are presented in Figure 7b. With the propagation velocities of “fast” paths being larger than those of “slow” paths, it follows that solution paths always map on a sequence of “fast” and “slow” path, intersecting at the intermediate state I and connecting I to the downstream state A and the upstream state B, respectively. We have observed that solving eq 8 with respect to the original variables C1 and C2 (which is more complicated from a mathematical and a computational point of view) the calculated characteristics map on the “fast” and “slow” paths such that for the remainder of this work we call these paths Γ1 and Γ2 characteristics with the corresponding propagation velocities λ1 and λ2. The (ln(xL1 ), SR)-plane presented in Figure 7 only shows the immiscible region. In order to assess the location of the characteristics with respect to the fluid phase composition, and to take the soluble region into account, Γ1 and Γ2 characteristics determined in Figure 7b are mapped into the ternary diagram in Figure 9. In this diagram, it can be noted that all Γ2 characteristics and a few Γ1 characteristics intersect with the binodal curve. Γ1 characteristics which do not intersect with the binodal can only be accessed from the miscible region through shocks fulfilling eq 36. In turn, characteristics in the immiscible region intersecting with the binodal are continued in the soluble region along the Γ1 and Γ2 characteristics derived in section 4. However, it should be noted that the slopes of Γi characteristics as well as the corresponding propagation velocities λi are discontinuous at the binodal curve. As a consequence, crossing the binodal curve along a characteristic either occurs through a shock transition or it produces a plateau

Figure 8. Propagation velocities λt (blue) and λnt (red) (a) along tieline characteristics Γt, (b) along non-tie-line characteristics Γnt. Continuous and dashed lines in (a) correspond to the highest and lowest tie-line (highest and lowest ln(xL1 )) illustrated in Figure 7. In (b), continuous and dashed lines correspond to the two pairs of nontie-line characteristics with watershed points located on the highest and lowest tie-line (two watershed points per tie-line). Open black circles indicate watershed points. Note that the non-tie-line characteristic with a watershed point at SR = 0.94 features a very narrow range of SR, such that only the watershed point and no curve is visible.

The states associated to these points are special because their propagation velocity in the physical plane is the same whether they occur within a transition whose image is a Γt or a Γnt characteristic. Therefore, such points can be reached through a path along a Γt (Γnt) characteristic and departed in the direction of a Γnt (Γt) characteristic without its propagation velocities changing. As a consequence and contrary to what happens at all other points, a change of path through a watershed point can occur without any intermediate constant state emerging in the solution. J

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 10. Initial and feed conditions (mole fractions) of the five different exemplary cases solved by equilibrium theory. All cases constitute displacement cycles with the same initial state A (mixture of inert components 2 and 3) and different feed states B1−B5 (B1−B3 being mixtures of adsorbing component 1 and inert component 2, B4 and B5 being located in the immiscible region). Connecting lines are for illustrative purposes only.

Figure 9. Mapping of Γ1 and Γ2 characteristics (compare Figure 7b) in the ternary diagram.

in the concentration profiles which does not correspond to a classic intermediate state, as the solution path remains on the same characteristic. For detailed explanations on the crossing of the binodal curve, see the Supporting Information. It is worth making one last comment on the impact of adsorption on the shape of the characteristics in the hodograph plane. With the S-shaped fractional flow functions (see 18), there is one point on every tie-line at which f X = SX (apart from SX = 0 or 1); thus, Fi = Ci. At this point, both phases travel with the same interstitial velocity. Connecting all these points over all tie-lines, one obtains the “equivelocity” curve, which intersects all tie-lines and ends in the plait point. In the absence of adsorption effects, it can be shown that both the equivelocity curve, as well as the binodal curve are non-tie-line characteristics.6 In contrast, in our case, accounting for adsorption effects and for the resulting volume (porosity) changes, it can be shown that neither the equivelocity curve nor the binodal curve constitute non-tie-line paths (proof is provided in the Supporting Information and in the literature).17 5.3. Solution of Exemplary Cases. The results discussed in the previous section emphasize the complexity and diversity of possible solutions for different initial and feed states, featuring (i) different possibilities of crossing the binodal, (ii) characteristics with peculiar shape, and (iii) directional derivatives of propagation velocities changing sign along these characteristics, possibly several times, thus producing complex semishock transitions. We want to provide a few exemplary elution profiles and their mapping in the hodograph plane, but for the sake of brevity, we only highlight the most relevant mathematical and physical aspects. A detailed explanation is provided in the Supporting Information. With reference to Figure 10, we consider 5 exemplary cases, all of them sharing the same initial state A (mixture of inert components 2 and 3), but featuring different feed states, namely B1−B5. The first three feed states are located in the miscible region, more precisely they are binary mixtures of adsorbing component 1 and inert component 2. Feed states B4 and B5 are located in the immiscible region. Paths in the hodograph plane, concentration and flow profiles (entire chromatographic cycles) for examples 1−3 (miscible feed states) are provided in Figures 11−13, those for examples 4 and 5 are in Figures 14 and 15, respectively. Note that, apart from the equilibrium theory solutions, concentration and flow

profiles of Figures 11−15 also exhibit thin black lines, corresponding to simulations based on the lumped kinetic model, which will be discussed in detail in section 6. First of all, it is worth noting that, with the assumption of different interstitial velocities of the convective phases, Fi differs from Ci and thus concentration and flow profiles differ. While concentration profiles describe the overall composition of the phase mixture at a certain position in the column (here at the column end where ξ = 1), flow profiles describe the composition of the overall flux (in this case again at ξ = 1). Experimentally measured elution profiles (e.g., obtained by gathering the eluate in fractions and analyzing these fractions, as it is done in the companion paper19) thus correspond to the flow profiles. Let us first focus on the cases with feed states in the miscible region. Comparing Figures 11−13, it can be noted that, although changes in the feed composition are rather small, the paths in the hodograph planes differ considerably in the three cases. Yet, all cases exhibit one unique solution, which in the three examples cross the immiscible region, during adsorption and desorption in cases 1 and 2, and during desorption in case 3. Concerning the adsorption step, one possible solution path entirely located in the soluble region is always available for initial states located on the x1 = 0 axis, connected to feed states located on the x3 = 0 axis. This path is the same solution that would be obtained using the method described in section 4.3, with an intermediate state located in the x2 = 1 corner and being reached through a contact discontinuity propagating at the interstitial velocity. Under certain conditions, a second path, crossing the immiscible region, should be considered, where the initial state is connected to a state in the immiscible region through a shock (or semishock). This second solution is the physically valid solution if the shock (semishock) emanates from the initial state with a propagation velocity faster than the interstitial velocity (ϵ−1) and enables a further valid solution path connecting to the feed state. In examples 1 and 2, such transitions propagating faster than the interstitial velocity exist; hence, we observe adsorption steps leading through the immiscible region. On the contrary, for example 3, the solution through the miscible region is valid. K

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 12. As in Figure 11, but with a feed state B2 (located on the c3 = 0 axis).

Figure 11. Solutions for exemplary case 1 with initial state A (located on the c1 = 0 axis) and feed state B1 (located in the x1 = 1 corner). Top: paths in the hodograph plane resulting from the equilibrium theory derivation (mass fractions). Center: concentration profiles. Bottom: flow profiles. Thick, colored lines in the center and bottom graph provide the solution of the equilibrium theory model, thin black lines correspond to numerical simulations with a number of discretization steps N = 200.

Two remarks are worth making. First, due to the different interstitial velocities of the convective phases, states are able to propagate faster than the interstitial velocity and can thus appear at the column outlet before the void time (which is impossible when a single-phase flow is present, as it is L

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 13. As in Figure 11, but with a feed state B3 (located on the c3 = 0 axis).

Figure 14. As in Figure 11, but with a feed state B4 (located in the immiscible region).

commonly the case). This property can also be concluded from Figure 8a, as propagation velocities of both tie-line and non-tieline characteristics feature values well above the interstitial velocity. Second, certain states in the miscible region can be connected through two possible paths, one leading entirely

through the miscible region and one crossing the immiscible region. Neither of these two paths leads to an apparent physical inconsistency, but the physically correct path is only one, namely, that featuring the faster transition emanating from the initial state. M

DOI: 10.1021/acs.iecr.7b05154 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

intermediate states when assuming equal velocities of the convective phases (see examples 1 and 2 in section 4.3). The mathematical reason for this can again be found in the propagation velocities λt along tie-line characteristics (compare Figure 8a), which approach 0 when approaching the binodal (i.e., as SR → 0 or SR → 1). Physically, this is explained by the fact that liquid phases at a low phase saturation (