Bifurcation analysis and steady-state multiplicity of multicomponent, non-equilibrium distillation processes

Chcmicol Engweering Scwnce, Vol. 46, No. 7, pp. 1757-1769, Printed in Great Britain.

1991.

WOY ~2509/91 $3 00 + 0.00 :c 199 1 Pergamon Press plc

BIFURCATION ANALYSIS AND STEADY-STATE MULTIPLICITY OF MULTICOMPONENT, NON-EQUILIBRIUM DISTILLATION PROCESSES

Institut

fiir Systemdynamik

A. KIENLE and W. MARQUARDT+ und Regehmgstechnik, UniversitHt Stuttgart, Pfaffenwaldring X0, Germany

(First received 1 June 1990; acceptedjbr

publication in reoisedform 5 October

9,700O Stuttgart

1990)

AbstractPSteady-state multiplicity of multicomponent distillation processes is investigated by means of a continuous mass transfer model. As a first step towards the analysis of complex columns a single column section is considered. The structural properties of the steady-state concentration profiles in the limit of an infinite spatial domain are closely related to the global material balances of the column section. It is shown that there exists a direct relationship between multiple steady-state concentration profiles and multiple solutions to the global material balances. Multiplicity of the solution of the global material balances is investigated using the Newton polyhedron method. The theory is applied to ternary mixtures with an ideal vapour phase and an ideal as well as a non-ideal liquid phase. For the former system multiple solutions occur only on thz boundaries of the concentration state space; the corresponding concentration profiles are therefore only of theoretical interest. In contrast, mixtures with a nonideal liquid phase homogeneous as well as heterogeneous-are shown to reveal multiple steady states in various situations. The results of the theoretical investigations are confirmed by simulation.

lation

INTRODUCTION

The first

simulation

of multiple

steady

states

in a dis-

tillation column was reported by Magnussen et al. (1979) for azeotrope distillation of some ternary alcohol-water mixtures, which exhibit a miscibility gap outside the concentration domain considered. The same phenomenon has been found by Prokopakis and Seider (1981) for some other ternary mixtures with similar properties. “Multiple steady states” mean, in this context, that the concentration and temperature profiles inside a distillation column arc ambiguous for some fixed system paramctcrs and input variables. This should not be mixed with so-

called design multiplicities, where a certain product specification can he obtained with different plant parameters or operating conditions. Such design multiplicities have been reported for example by Chavez et al. (1986) for a two-column system. The results of Magnussen et al. (1979) have caused significant research effort in recent years. Since the thermodynamic properties of the mixtures play an important role on the qualitative behaviour of separation processes, the main emphasis has been first on the analysis of vapour-liquid equilibria. Lucia (1986) proved the uniqueness of solutions to single-stage isobaric flash processes involving homogeneous mixtures. The behaviour of heterogeneous mixtures was considered by van Dongen et al. (1983). They have found multiple solutions of the equilibrium relations and claimed that this was a possible reason for multiple steady states in heterogeneous azeotropic distil-

‘Author CES 46:7-C

to whom

correspondence

should

be addressed.

columns.

solutions

The

to the phase

practical

relevance

equilibrium

of multiple

of heterogeneous

mixtures was discussed controversially by Michelsen (1984) and Doherty et al. (1984) because multiple solutions occur only in the region of a metastable liquid phase. Other investigations were concerned with the theoretical analysis of complete columd models. A comprehensive review of the literature was given by Doherty and Perkins (1982). They further proved that multiple steady states can not occur in binary homogeneous distillation and multicomponent homogeneous fash distillaton. Their models assumed constant molar overflows and phase equilibrium on the trays. The results for binary homogeneous multistage distillations were confirmed by Sridhar and Lucia (1989). Additionally they proved that uniqueness still

holds for a more detailed model, which includes the energy balance and hence relaxes the assumption of constant molar overflows. A generalization of the techniques and results to multicomponent multistage distillation processes is difficult and only possible under further restrictions (Sridhar et al., 1989). The very recent results of Sridhar et al. (1989) showed some evidence that homogeneity of the liquid phase may be a prerequisite for the uniqueness of the steady state not only in the binary but also in the multicomponent case. So far, these investigations have provided only a little insight into the physical reasons causing steady-state multiplicity in distillation processes. This may be due to the fact that previous approaches have been either based on very simple situations, i.e. a single equilibrium stage, or on rather complex processes such as complete columns or even interlinked 1757

A. KIENLE

1758

and W. MARQUARD?

systems. We suggest here to study a model system, which is of practical relevance on one hand and which allows on the other hand a rigorous (quasi-)analytical treatment of the model equations to yield as much insight as possible. For that purpose we start with a division of the whole distillation column into its interconnected subsystems, which are a condenser, a reboiler, a feedstage and (usually) two column sections as depicted in Fig. 1. The first naturally arising question is whether a single column section is able to exhibit multiple steady states or whether this phenomenon is caused by the interaction between these different subsystems. Such interaction phenomena between different dynamic subsystems are known to be the reasons for instability and steadystate multiplicity in many chemically reacting systems as well as in many other systems in scieuce and engineering. As a first step the analysis of a single column section is carried out infthis paper, since its understanding is a necessary prerequisite for the study of interaction phenomena in columns or even more complex plants. In contrast to previous investigations the column section is modeled by a differential mass transfer model instead of the commonly employed tray-bytray model. Such a model is usually used to describe packed columns but is also convenient for the approximate modelling of tray columns (Mecklenburg and Hartland, 1975). The column section is assumed to be infinitely long in the direction of flow. The mathematical treatment of the system simplifies that way, because phase equilibrium is attained at both ends of the column section. Thus, an analytical treatment of the model equations becomes tractable, which is indispensable for elucidating the possible phenomena in a high-dimensional space of model parameters. As it will turn out, this simplification allows a further analysis of the physical reasons for potential multiplicity closely related to phase equilibrium properties. According to Sridhar and Lucia (1989) and Sridhar et al. (1989) we neglect the energy balance and consider only multicomponent mixtures. Though our focus will be on ternary mixtures here, the methods presented

are also applicable components.

to separations

MODEL

with any number

of

EQUATIONS

Consider a section of a distillation column as presented in Fig. 1. The compositions of the liquid and vapour phase are given by the mole fractions xi and yi, respectively. Let A be the ratio of the molar flow rates in the liquid and the vapour phase. Assuming: l l l l

constant molar flow rates, thermal equilibrium, negligible axial dispersion, and constant pressure,

a dimensionless mathematical modei of the column section is given by the steady-state material balances for N - 1 linearly independent compositions in both phases: dxi -

--jai@,

dZ

i = I(l)N

- 1,

YX

-aI<
where &(x, y) represents the mass transfer tween the two phases. The boundary conditions

rate be-

lim xi(z) = xi+ i-m (3) lim y,(z) = yi_ z---o0 are determined by the given inlet concentrations xi+ and yi_. For an infinitely extended column section the vapour outlet is in equilibrium with the liquid inlet at the top of the section and the vapour inlet is in equilibrium with the liquid outlet at the bottom of the section according to yi+ =A@+) (4) yi_ =fi(x_). Model (lH3) can be simplified by the fact that the concentrations J+(Z) are linearty dependent on xi(z). This is easily shown by integration of the sum of eqs (1) and (2) from the lower boundary of the section to some coordinate z: yi -X(x_)

= A(& - Xi_).

(5)

Thus the system dimension can be reduced for the infinitely long system from 2(N ~ 1) to N - 1 by eliminating the vapour compositions yi in eq. (I) by means of ey. (5) to yield

3_. -.LiCx*

Fig. 1. Infinitely extended

column

section.

f(x-) + A(x - x-)l~ i=

l(1)N - 1.

(6)

Bifurcation

analysis

and steady-state

Integration from z = -co is arbitrary. The choice z = + a, would lead to a similar model, depending on x+ instead of x- _ The multiplicity analysis is carried out with a special and rather simple expression for the mass transfer function: hi

= :

CA(X)

-

Yil,

i =

l(l)N

- 1

(7)

where Bi is a dimensionless mass transfer coefficient. This relation follows from film theory. The assumptions, which results in equal transfer coefficients for all components. were discussed by Krishna and Standart (1979). Serwinski and Gorak (1983) investigated assumptions which result in different Bi. With eq. (7) the more general eq. (6) simplifies to

2 =:

1759

multiplicity

in the ternary case (N - 1 = 2). Rhee and Amundson (1974) provided such a solution for a similar problem in fixed-bed adsorption. The qualitative behaviout of the system, however, can be studied by an analysis of the trajectories in the (N - 1)-dimensional concentration phase space. The features of these trajectories are discussed here for the ternary case only to allow a geometrical representation of the trajectories in the two-dimensional phase plane. An extension to systems with more components is straightforward but limited with respect to geometrical representation. Division of eqs (8) yields the gradient field of the system in the phase plane:

h

4C_W -_Ux-) - Ah - x2-)1

dx,= B,Cf,(x) -fib-1

id4 - 4x1 - Xl-11=iG (10)

[j(x)

-j(x_)

= A($ - Xi_)], i =

l(l)N

- 1.

(8)

The right-hand side of eq. (8) represents the transfer rate jOi, which now depends only on the liquid composition x. The non-linearity of the system is given by the equilibrium relationf;, which is in general a function of all mole fractions xi. Solutions to eqs (8) or (6) with boundary conditions (3) and equilibrium conditions (4) do not exist for arbitrary feed concentrations xi+ and yi_ . As a necessary condition for the existence of such a solution the input concentrations xi+ and yi_ have to satisfy the global material balances under equilibrium conditions Ji‘(x+) -X(x_)

- A&+

-xi-)

= 0, i = l(l)N

- 1.

(9)

For a given mole rate ratio A and given compositions at any of the two boundaries, the compositions at the opposite boundary are determined by eq. (9). EXISTENCE OF MULTIPLE

SOLUTIONS

One potential reason for multiplicity of the concentration profiles could he multiplicity in the boundary concentrations determined by the boundary conditions (3) and the equilibrium conditions (4). This possibility can be ruled out easily for homogeneous mixtures because of the uniqueness of the phase equilibrium (Lucia, 1986). According to the results of van Dongen et al. (1983) we expect the same behaviour for potentially heterogeneous mixtures in the region of a stable liquid phase. Hence, we confine our considerations to multiplicities of the concentration profiles inside the column section with fixed boundary values. This requires the analysis of the two-point boundary value problem (8) and (9). The corresponding initialvalue problem as a model for the dynamic behaviour of the so-called simple distillation was treated intensively by Doherty and Perkins (1978a, b, 1979). Some facts for ternary mixture.5 Only in a few special cases a closed-form solution of the non-linear differential eqs (8) can be found even

The curves jO; = 0 (i = 1,2) play an important role in the system’s behaviour. They divide the concentration triangle defined by x,,xz>O X1 +x2

(Ill < 1

(12)

into regions where the trajectories exhibit qualitatively different slopes. For the intersection of the curves j., = 0 and j,, = 0 both transfer rates vanish simultaneously. Phase equilibrium is established. These points are called the singular points of the system (subscript “s”). They represent all possible boundary concentrations of the infinite column section as they fulfill the global material balances (9). In the remainder we only consider the special case of equal transfer coefficients, which makes the analysis much more easier. It is expected that this is no incisive restriction, because it is known for ideal mixtures that the ratio of the transfer coefficients does not play any role in the existence and number of solutions (Rhee and Ammundson, 1974; Kienle, 1989) and it is expected that this property still holds for non-ideal mixtures. For equal transfer coefficients B1 =E*

= B

(13)

further statements about the nature of the singular points as well as the trajectories in their neighbourhood are easily possible. The eigenvalues ui of the linearized system (8) are given as the solution of u

B afl -A ax2

B& A

ax,

The eigenvalues ui can be obviously related to the mole rate ratio A and the eigenvalues li of the linearized equitibrium relation:

A.

1760

KIENLE

and

W. MARQUARDT I

ni = ;

(ni -

A),

i=

1,2.

The eigenvalues oi are real and distinct provided that the eigenvalues Li of the equilibrium relation are real and distinct. The results of Kvaalen et al. (1985) and Doherty (1985) revealed that all ii are always real and further positive at any point of the concentration triangle, where the two phases are thermodynamically stable. The proofs require the assumption of constant latent heats. Thus the singular points are always nodes or saddles for such mixtures. Kvaalen et al. (1985, p. 1201) claimed that the eigenvalues & inside the composition space are usually distinct and confined multiple eigenvalues to “extremely special cases”. As we will show below, however, multiple eigenvalues may arise quite frequently in vapourliquid-equilibria. At the singular points the right-hand side of eq. (10) becomes indeterminate. Application of L’Hospital’s rule yields

which is a quadratic equation for the unknown slope in the singular points. For equal transfer coefficients this equation has always two solutions, which are real and distinct in all nonsingular cases, as the discriminant of eq. (17) becomes identically to the discriminant of the eigenvalue problem (15), which always has only real solutions (Kvaalen et al., 1985). With the knowledge of the gradient field, the singular points, their eigenvalues and the slope of the trajectories in the singular points, the trajectories can be constructed in the whole concentration triangle. Multiplicities of the global material balances Multiple steady states are marked by the fact that the two singular points at the boundaries of the system, which are fixed by the corresponding input flow rates, are connected by more than one trajectory in the same direction. This is likely to occur if the system exhibits more than two singular points. Then the global material balances are not unique. This is illustrated by a constructed example in Fig. 2. According to Fig. 2 we have to distinguish between two different kinds of multiplicities. On one hand the same two singular points (A and B in Fig. 2) are connected mukiply by different trajectories; on the other hand different solutions may be established with different pairs of singular points (A-B and A-C in Fig. 2). The latter does not represent steady-state multiplicity in the sense of our initial definition, as they require different input specifications. In this context we are therefore interested only in the first kind.

~I C

x2

B

A

0

I

XI

Fig. 2. Trajectories for a system performing multiple steady states (schematic).

As an analytical treatment of the non-linear twopoint boundary value problem is impossible in generaI, we confine ourselves first to the analysis of the global material balances, which form a system of nonlinear equations. Suitable mathematical tools for the analysis of solution multiplicity of such systems are available from bifurcation theory. Afterwards we go back to the phase plane and study the effect of multiple singular points on the trajectories of the distributed parameter system. The bifurcation analysis of the global material balances (9) is c_arried out with the Newton polyhedron method. The next section summarizes some of the results of Kouchnirenko (1976), which are necessary for our application. Arnold (1981) gave a comprehensive review on the development of this theory. The application of this method to chemically reacting systems has been demonstrated in a recent paper by Lyberatos et al. (1984). The main advantage of this theory in comparison to classical singularity and catastrophe theory (Golubitsky and SchBifer, 1985) is its direct application to a set of algebraic equations without the requirement of its transformation or reduction to a single equation. The application of this method to two-dimensional systems, which correspond to the global material balances for ternary mixtures, leads to quite simple criteria for the occurrence of bifurcation points. This facilitates a physical interpretation and finally a clear statement on the feasibility of the resulting criteria for meaningful model parameters. The Newton polyhedron method. For simplicity the method is described here only for a set of two algebraic equations. In this case its application becomes particularly simple and requires only a geometric construction. An extension to higher-dimensional systems is also possible. The interested reader is referred to Kouchnirenko (1976) and Lyberatos et al. (1984).

Bifurcation

analysis

and steady-state

multiplicity

1761

The theory is a local one, which means that any statement on the number of solutions only holds in some neighbourhood of a certain solution x0. Lyberatos et al. (1984) illustrate this fact by the almost classical “pitchfork” presented in Fig. 3: x3 - Px = 0.

(W

Equation (18) shows only one solution branch (x = 0), if the parameter P is less than zero. Two more solution branches (x = i fi) occur for P greater than zero. Except for P = 0, there exists not more than one solution in a local neighborhood of any point of the solution graph. Hence, the so-called multiplicity+ of this solutions is 1. In the neighbourhood of P = 0, there exist at most three solutions for small perturbations of the parameier P. Thus the multiplicity of the solution x = 0 is therefore said to be 3. P = 0 is the bifurcation point of eq. (18). Let us consider a system of two non-linear equations in two unknowns, x1 and x2: gj(xl,X2)

=O,

j=

1,2

Fig. 3. Pitchfork.

i

h

Oi

(1%

to introduce

the basic definitions. Originally, the theory has been developed for polynomial systems only. It can be applied, however, to all kinds of non-linear systems provided that the finctions are analytic and hence can be expressed by a Taylor series near a certain solution x0 (Kouchnirenko, 1976):

i l[ &* (“)&I,,

gj(xl, x2) = &- z

k-i

Fig. 4. Newron x Ax\Ax:-’

, I

j = 1,2

(20)

diagram of a function with non-vanishing derivatives (a = ,O = 1).

where Axi = xi - xi,,.

(21)

The Newton diagram of such a Taylor series is the graph of a function which maps the natural numbers N into a two-dimensonal product space ZEN x N. Z consists of all ordered pairs {i, k - i} of powers corresponding to non-vanishing partial derivatives in eq. (20). If none of the derivatives of gj vanishes, we obtain a Newton diagram as shown in Fig. 4. The convex envelope of this graph (see Fig. 4) forms a polyhedron, to which the name of the method refers. In the two-dimensional case the polyhedron degenerates to a polygon. The Newton number v is in general defined with respect to the Newton polyhedra of the functions gi and those of the composite map of gl.qz. In the two-dimensional case it can be determined directly using only the Newton diagrnms of g1 and g2.

‘This is a different multiplicity definition as usually employed. The usual definition of multiplicity refers to the situation where multiple solutions occur for one set of parameters. Multiplicity of soluti&s in bifurcation theory refers to the local appearance of additional solutions and is therefore related to the bifurcation points.

Let sj and pj be the sections of the axes under the polyhedron in the Newton diagram of gj. Then the Newton number is given as v = min(alD2,azfi1).

(22)

For a system of two equations, where none of the partial derivatives vanishes, both functions reveal a Newton diagram as shown in Fig. 4. Thus application of eq. (22) yields a Newton number v of 1. According to Kouchnirenko (1976) the following central theorem of the theory relates the Newton number to the multiplicity of the solution x0: Theorem 1: In all non-degenerate cases the maximum number of solutions in the uicinicty of x0 for small parameter perturbations is equal to the Newton number v. In all degenerate cases the maximum number of solutions is greater than v. Some final comments on the so-called degenerate cases are appropriate. In non-degenerate cases the occurrence of bifurcation points is revealed by the topological properties of the polyhedron. For vanishing low-order terms the sections of the axes under the

1762

A.

KIENLE

and W.

polyhedron are greater than 1, which leads to a Newton number greater than 1, indicating a bifurcation point. It is possible, however, that a bifurcation point occurs for some parameter values without any noticeable change in the structure of the Taylor series and hence the geometric properties of the polyhedron. An exact statement on the number of solutions is therefore not possible with the polyhedron method. These cases are degenerate in the sense of this theory. For the formulation of the degeneracy conditions we refer to Kouchnirenko (1976) or Lyberatos er al. (1984). Application to the global material balances. The material balances (9) are first expanded in a Taylor series in the vicinity of a point x0, which has to be a solution of eqs (9). We choose x0 to be the concentration x+ at the upper boundary of the column

MARQUARDT

section. The concentrations of the section are treated order lation

which close

x_ at the lower boundary as parameters. Only firstterms have to be regarded first, since the calcuof the Newton number depends on the points lie on the face of the corresponding polyhedron to the origin. Expansion of eqs (9) yields

If none of the elements Jik of the Jacobian J vanish, the Newton number is equal to 1, as demonstrated in the last section (see also Fig. 4). Hence multiplicity of the solution is only possible for the degenerate case. If only one of the elements Jil, vanishes and the secondorder derivatives in expression (23) are non-zero the Newton number remains equal to 1 for all non-degenerate cases (see Fig. 5). Analysis of every possible combination of vanishing elements J,, yields that only 10 out of 16 cases may reveal multiplicity. In complete analogy to Application 1 of Lyberatos et al. (1984) the conditions in Table 1 are obtained. The first case exhibits bifurcation only for degenerate parameter conditions, the others also for non-degenerate conditions. The degeneracy conditions arc compiled in column 4 of Table 1. Multiplicities of more than 1, more than 2 or more than 4 can occur for the conditions in rows 1, 2-5 and 10, respectively.

;,lgjg_ )ijg* 0

aI

k-i

I

0

O21

k-i

Fig. 5. Newton diagrams for a system of two equations with J,z = 0.

Table 1. Classification of the bifurcation points of the global material balances Case

Newton number v

1

1

2

4 5

10

4

Conditions for the elements -lit

Jir, # 0

i, k = 1,2

Degeneracy condition det J = 0

Bifurcation

analysis

and steady-state

The classification in Table 1 is not complete, since we disregarded the cases where additionally to the first-order terms the corresponding terms up to an order of equal to or greater than 2 vanish simultaneously in eq. (23). This behaviour is expected to be rather unlikely for physically relevant phase equilibrium relations and is therefore not further considered in the current state of our investigations. The conditions in Table 1 are formal in the sense that they are valid for all systems of two equations, if the derivatives of the equilibrium relation in the degeneracy conditions are replaced by the derivatives of the corresponding functions of the general eq. (19). It remains to analyze whether these criteria can be satisfied for real mixtures and feasible operation conditions. The main emphasis will be on the physical properties of the equilibrium relation and its derivatives as they determine the elements of the Jacobian directly. General statements on the existence of the bifurcation points listed in Table 1 are possible by calculating the eigenvalues li and ui from eqs (14) and (15) with respect to the conditions in Table 1. The results are listed in Table 2. The first case refers to the degeneracy condition of case 1 in Table 1. The others refer to the conditions in column 3 of Table 1. The following conclusions can be drawn. Degenerate bifurcations with a multiplicity greater than 1 in row I occur for any homogeneous mixture with positive eigenvalues for a suitable mole flux ratio A. We will discuss this case in detail later. The existence of the other (non-degenerate) bifurcation points depends on some special thermodynamical properties of the equilibrium relation. For cases 2-7 the identity be.tween the eigenvalues of the equilibrium relation Ai and the partial derivatives af/axi is required. The rest (g-10) even require multiple eigenvalues of the equilibrium

relation.

The significance of the properties of the equilibrium relation is pointed out by these results. The next section presents the results for two frequently used vapour-liquid equilibrium models.

RESULTS

FOR

Ideal liquid phase For ideal mixtures

TERNARY

MIXTURES

the equilibrium

Table

Case

Newton number v

relation

2.

Eigenvalues

relation

1763

,_p+

,

a, > Ez

>

.

. . > aN = 1.

35 # 0.

With the definitions of the elements J in eq. (23) the conditions Jll

(27) J12rJ2L

ill

2}

Eigenvalues of the linearized system IT;

)_A!4

’ I

0

points

2}, i #j

ax,

#

follow for the JitS. Thus comparison with Table 1 shows that the maximum possible Newton number is 1. For ideal mixtures no local multiplicities are therefore possible except in the degenerate case of row 1. Further results can not be attained with the Newton polyhedron method. Since the global material balances for ideal mixtures form polynomial equations, powerful tools are available from algebraic geometry for further investigation. A simple analysis of the curves joi = 0 in the phase plane shows that there exist only three different types of phase diagrams, which characterize the steady-state behaviour of the ternary system com-

[rl = uz = 0

’

JiliLof the Jacobian

+ Jzz

i,je{l,

&i=*

(25)

ax,’ ax,

&=&=A

2

(24)

(26)

a

i,jE{l,2},i#j

2-7

l(1)N

The relative volatilities are functions of the boiling temperature, even for ideal mixtures. In most cases, however, this dependency is negligible. The relative volatilities are assumed to be constant in the subsequent discussion. With this assumption the equilibrium relation (24) is formally equivalent to the adsorption equilibrium of Langmuir isotherms and the model equations become equivalent to the two-solute adsorption discussed by Rhee and Amundson (1974). Differentiation of eq. (24) and the precondition (25) result in the following properties for the equilibrium relation of an ideal, ternary mixture:

it{l,2)

ic(l,

i =

explicit function of the liquid compositions, ai is the volatility of component i relative to the boiling component. The components are numas follows:

(Ti = 0

& = A

10

is an where lowest bered

Eigenvalues of the

1

2

A=

li and oi at the bikxcation

equilibrium

1

8 and 9

multiplicity

axj

1764

A.

KIENLE

and

pletely (Kienle, 1989). The first two cases, which have been reported Rhee and (1974), exhibit two singular within the tion triangle. corresponding phase are shown Figs 6 7. The cases differ respect to behaviour of 2. In first case the phases are enriched with component 2 from the bottom to the top, whereas in the second case component 2 is enriched from the top to the bottom of the section. Following the arguments of Rhee and Amundson (1974), it can be proved that in both cases a unique solution exists, as the two singular points x_ and x+ are connected by exactly one trajectory. This result holds also for any ratio of transfer coefficients (Rhee and Amundson, 1974; Kienle, 1989).

W. MARQUARDT

A third case which was not reported Amundson exists for 0 =fi(x-)

Ax,_.

(281

The corresponding phase diagram is shown in Fig. 8. For this case the system exhibits always three singular points within the concentration triangle. Two of them are always on the boundary x2 = 0 of the concentration triangle and the third is inside. If we choose the two binary singular points as the boundary points of the column section by specifying the corresponding input flow rates at the bottom and the top of the section, we obtain infinitely many trajectories between the two singular points as shown in Fig. 8. This means that there exist infinitely many steady states for the column section.

I II

1. I.~, ja2r0 II: I., ’ 0. io*
-

by Rhee and

‘n

joIT 1,2 >O j,, > 0, ja2 <

0

In: 10l<0,102’0

m: j,,,j,,c

i.2<0

0

x

m

0

I XI

Fig. 6. Phase

diagram

for

an

ideal

mixture

with

Fig. 8. Phase

diagram for an ideal 0 =h(x_) - Ax,_.

mixture

with

II : jo,>O, jaz CO lil: jo,-=O, jo2 2-O

~:j,l,L~cO

Fig. 7. Phase

diagram A,(x)

for < A

an


ideal

mixture

with

Fig. 9. Computed

concentration O=f2(x_)-kc-.

profiles

for

Bifurcation Table

3. Parameter

values

analysis

and steady-state

multiplicity

1765

for the

A(x> T, = y,i(T, I.1537 4 2 0.067825 0.65946 scm

A

a1 12 XIx2A, = R,

This result is verified by steady-state simulation with the simulation environment DIVA (KrBner et al., 1990). The parameters of the example are presented in Table 3. For an approximate simulation of an infinitely long column section the dimensionless transfer coefficient B

B,E

Non-ideal

liquid phase

The equilibrium relation for a mixture with an ideal vapour phase and a non-ideal liquid phase is given by

i = l(l)N.

(30)

In contrast to the relation discussed in the previous section, the equilibrium composition depends now on the boiling temperature. Usually the boiling temperature is determined by the closure condition

i&!x,n -

1 = 0.

(31)

This is in general a highly non-linear equation. Thus the boiling temperature can not be eliminated explicitly in eqs (30). But for the bifurcation analysis it is more convenient to regard the equilibrium composition of the vapour phase only as a function of the liquid composition: A =A[%

V

which depends on the column length 1, the transfer coefficient of the vapour phase k, the specific interfacial area a, and the vapour flow rate Y, is enlarged until equilibrium is established at both ends of the column section. As an equilibrium criterion the gradients of the concentration profiles at the boundaries are observed. They have to vanish according to eqs (6). The result obtained this way are considered to be a reasonable approximation of the behaviour of a computationally infeasible infinitely extended column section. Figure 9 shows some of the solutions represented by the trajectories in Fig. 8. As the trajectories cross the curve j,, = 0, the corresponding concentration profiles of component 2 exhibit an extreme value according to eq. (6). The trajectory ABC runs inside the column section through another singular point, in which phase equilibrium is established; the concentration profiles therefore reach a pinch state within the section. Dynamic simulations of the transient behaviour show that all steady states are dynamically stable. This is evident because the initial amount of component 2 has to remain in the section, as the outlet as well as the inlet streams are represented by the two binary singular points A and C in the concentration triangle. No component 2 is transported in or out the column section. These results are not very surprising. They correspond very well to our physical understanding of the process. Different initial amounts of the intermediate boiling component in the column section will lead to dimerent solution structures for the same operating conditions if the leaving and entering streams only “block” the column with respect to this component. The investigation of non-ideal mixtures in the next section will reveal, however, that a similar effect is even possible with non-zero concentrations of all components in the leaving and entering streams.

X)X;,

W01

(32)

because the temperature in eq. (3 1) is also a function of the liquid composition. The vapour pressure of the pure components is calculated by the Antoine equation in the form ]n

!I!% = c. ,I Pit

ciz

_

T

+

cij

’

i = l(l)N.

(33)

The non-ideality of the liquid phase is modeled with the activity coefficients yi. As a first step we use here one of the simplest correlations for the calculation of the activity coefficients, which is the so-called Regular Solution or Two-Suffix Margules Model: RTln

yk = 0.5 g i: (A,, + Aij - Ajk)xixj. i=]

(34)

j=l

The binary interaction by

parameters

Aij are constrained

A, = A,,, Aii = 0, i, j = l(l)N.

(35)

Our own investigations as well as those of Malesinski (1965) and many others show that the predictions of this model are qualitatively and quantitatively reliable for mixtures of non-polar molecules of similar size, when the Aij are fitted to experimental data. Again, the first step of the bifurcation analysis is the calculation of the partial derivatives of the equilibrium relation (30) (see Appendix A). The complexity of the resulting expressions does not allow precise statements on the properties of the derivatives. The only general statements from thermodynamics about their properties are that each of them has to be finite (Doherty and Perkins, 1978a, b) and that the matrix of the derivatives is positive definite (Doherty, 1985) if further constant latent heats are assumed. Hence, it is impossible to determine whether the criteria for the bifurcation points given in Table 1 can be satisfied in general for physically feasible parameters. Instead, we show the existence of bifurcation points exemplarily by numerical calculation for some different mixtures. The only free system parameters for a given mixture are the flow rate ratio A and the state x_ at the bottom of the column section, if the pressure p affecting the equilibrium relation (30) is assumed to

1766

A.

KIENLE

and W.

MARQUARDT

0

I

XI A-O.85808

XI A-O.83000

A=078500

Fig. 10. Behaviour of the curves joi = 0 in the neighbourhood of a bifurcation point with v = 4 for the homogeneous azeotropic mixture Z-methyl-Z-buten+thanethiol-methylformate.

be fixed at 1.013 x lo5 Pa. Potential bifurcation points x0 are fixed by the global material balances (9). Altogether, the bifurcation problem has the three independent parameters A, xi and xZ _ According to classical singularity and catastrophe theory (Golubitsky and Schgfer, 1985) the maximum number of soIutions is expected to be 4, which corresponds to a Newton number v equal to 4 for non-degenerate parameter conditions. Thus we are searching for a mixture which satisfies the conditions in row 10 of Table 1 for certain parameter values. These criteria represent four equations to determine the three unknowns A, x1 _ and x2_. The overdetermined problem is solved in the least-square sense. This is sufficient, since it is not necessary to calculate the bifurcation point exactly. A parameter set in some neighbourhood of this point will-according to Y = Gresuit in four solutions to the global material balances. For this purpose, some of the azeotropic mixtures presented in the book of Malesinski (1965) have been investigated. The result for the mixture 2-methyl-2buten_ethanethiotLmethylformate is shown in Table 4. Besides the calculated values of A, x1 _ and x2 _ , Table 4 also contains the binary interaction parameters of the Two-Suffix Margules equation.’ Thermodynamic stability has been checked by numerical computation of the Hessian of the molar Gibbs free energy at fixed temperature and pressure, which is given in the paper of van Dongen et al. (1983) for yi from the Two-Suffix Margules Model. Since thermal equilibrium is assumed, the temperature is taken to be the boiting-point temperature for the given pressure and concentration point under consideration.

+The parameters Aij we use here are different from those presented in Table IV, row 25, by Malesinski (1965) for the same mixture. According to Malesinski (1965) the best predictions of equilibrium data can be obtained if the A, are calculated from the temperatures of the binary azeotropes. This depends on the calculation of the vapour pressure of the pure components. Because Malesinski (1965) used another correlation for this purpose, this has been repeated by using the Antoine eq. (33). The Antoine parameters are taken from Henley and Seader (1981) and transformed to SI units.

Table 4. Calculated parameters for a bifurcation point with Y = 4 and the interaction parameters of the Two-Suffix Margules Model for the homogeneous azeotropic mixture 2-methyl-2-butenethanethiol-methylformate A

X10 x20 A u OcJ/mol) A,, (kJ/mol) A L3(kJ/moU

0.85808 0.15299 0.62682 1.30331 2.38239 3.93663

The behaviour of the curves jai = 0 in the neighbourhood of the bifurcation point is depicted in Fig. 10. Remember that the points of intersection between the curves j,, = 0 and j,, = 0 mark the solutions of the global material balances. At the bifurcation value of the flow rate ratio A four solutions collapse together into point x0. The reverse phenomenon can be visualized by decreasing the parameter A, which separates the four solutions. As shown in the last diagram of Fig. 10, further solutions may occur for suitable A. which clearly demonstrates the limitation of the local theory. Analogous behaviour can be observed for many other mixtures. Table 5 gives an overview on the mixtures we examined so far together with a classification of the ternary azeotropes according to Malesinski (1965). Some of the mixtures presented in Table 5 exhibit a miscibility gap within the concentration triangle. They were only considered in the domain of a stable Iiquid phase. Our investigations do not reveal any qualitative difference between homogeneous and potentially heterogeneous mixtures as found by Sridhar and Lucia (1989) and Sridhar et al. (1989). All these mixtures exhibit bifurcation points with a Newton number v equal to 4 within the domain of a stable liquid phase. The eigenvalues of the equilibrium relation at these bifurcation points have to be double as demonstrated above. Thus our investigations are in direct contradiction to the hypothesis of

1767

Bifurcation analysis and steady-state multiplicity

Table 5. List of mixtures examined Mixture Cyclohexane-benzene-methanol 2-Methyl-2-buten-thanethiol-methylformate Hexane-chloroform%thanol Methyl acetate-acetone-methanol Cyclohexane-acetone-methanol Methanol-acetone-chloroform n-Heptane-pyridine-acetic acid

Kvaalen et al. (1985) that multiple eigenvalues of the equilibrium relation only occur in “extremely special cases”. Figure 11 shows the phase diagram of 2-methyl2-buten-ethanethiol-methylformate corresponding to Fig. 10. If we choose the singular points A and C to be the boundary values an infinite numbers of trajectories may exist, which all connect the two singular points in the same direction. Some of them intersect the curve j,,, = 0, others the curve jbz = 0. The concentration profiles of component 1 corresponding to the former must reveal therefore a maximum along the spatial coordinate. A maximum is displayed in the spatial concentration profiles of component 2 with the latter group. One single trajectory does not intersect any of the curves jai = 0. Hence the corresponding concentration profile is monotone. The family of trajectories is bounded by the trajectories ABC and ADC. As they run through another singular point inside the column

section,

the corresponding

Ternary azeotropes

Miscibility gap

No ternary azeotrope Positive azeotrope Positive azeotrope Positive azeotrope Positive azeotrope Positive-negative azeotrope Positive-negative azeotrope

YeS No No NO YeS No Yes

I

0

I

XI

Fig.

11. Phase diagram corresponding to Fig. 10.

concen-

tration profiles exhibit a pinch region. Again these results are confirmed by steady-state simulation with DIVA. Figure 12 shows some of the calculated concentration profiles represented by the trajectories in Fig. 11 discussed above. The additional parameters are A = 0.785, B1 = Bz = 5000. Dynamic simulation reveals that steady-states corresponding to trajectories inside the family drift very slowly to one of the states corresponding to a trajectory at the boundaries of the family or to the trajectory, which does not intersect

any of the curves j,; = 0. Numerical

ties or the finite approximation tended spatial

domain

instabili-

of an infinitely

are expected

p__

;,

ex-

to be a potential

0

z

reason

for this behaviour. A clear statement on the stability of the different steady states based on dynamic simulation is therefore impossible. Instead, a detailed analysis is required, which will be the subject of further investigations.

DISCUSSION

A new theoretical reasons

approach

for steady-state

umns has been presented

to clarify

multiplicity

the physical

in distillation

col-

in this paper. The column is divided into its subsystems, which are a condenser, a reboiler, a feedstage and the different column sections, because interaction between the different subsystems is expected to be a potential reason for steady-state multiplicity.

Fig. 12. Computed concentration profiles corresponding to Fig. 11.

The analysis presented in this paper is concerned with a single column section. In our opinion this has to be the first step in building a theoretical under-

1768

A.

KIENLE and

standing of the steady-state behaviour of multicomponent distillation processes. The application of the methods developped in this paper to other countercurrent separation processes is straightforward, as exhibit structural properties they the same (Marquardt, 1988). To facilitate the analysis we restrict our attention first to the theoretical limit of an infinite spatial domain. It is shown that multiple steady states as well as extrema and pinches in the concentration profiles may occur as a consequence of multiple solutions of the global material balances. This behaviour depends on the thermodynamic properties of the mixtures. A theoretical analysis reveals that non-ideality of the liquid phase is a necessary condition for steady-state multiplicity of practical relevance. This result is in agreement with the results of other investigators (Sridhar and Lucia, 1989; Sridhar et al., 1989). In contrast, a numerical analysis of ternary mixtures with a non-ideal liquid phase modeled with the Two-Suffix Margules Model for the activity coefficients of the liquid phase reveals several new and, in our opinion, important properties. Multiple steady states have been calculated for mixtures which are homogeneous in the whole concentration triangle, as well as for mixtures which are potentially heterogeneous outside the concentration domain considered. This result relaxes the role of homogeneity in more realistic non-equilibrium distillation processes, which is claimed to be of central significance in equilibrium distillation (Sridhar and Lucia, 1989; Sridhar et al., 1989). Further, it is shown that multiple eigenvalues occur inside the concentration triangle for globally homogeneous and potentially heterogeneous mixtures. This result is also supposed to have strong consequences for the equilibrium approach, which is based on the mathematical theory of hyperbolic conservation equations (Kvaalen et al., 1985). All properties predicted by the Two-Suffix Margules Model are confirmed in a qualitative sense with more common and probably more realistic models for the activity coefficients of the liquid phase such as Wilson and WNIQUAC (Kienle and Marquardt, 1990). Steady-state multiplicity in the infinitely extended column section is marked by an infinite number of steady states, which is-though exotic-a well-known phenomenon in distributed reaction systems (Eigenberger, 1972; Aris, 1975; Witmer ef al., 1986). It is expected that this infinite number of steady states reduces to a few solutions by the influence of the boundaries in the finite extended system or some other effects that we neglected such as axial dispersion or the influence of the energy balance. But as these effects are shown to be secondary (Marquardt, 198S), we do not expect a completely different behaviour, i.e. the vanishing of steady-state multiplicity as a whole. The relevance of these results for steady-state multiplicity in distiltation columns is obvious. As it turned out that even a single column section is able to exhibit multiple steady states, this may be a potential reason for multiplicity in the whole column. Thus, the thermodynamic properties causing this phenomenon as

W.

MARQUARDI

well as the interaction phenomena in the column require further investigation to contribute to a physical explanation of steady-state multiplicity in distillation columns. The significance of these results on the dynamic behaviour of distillation processes, which is marked by non-linear travelling waves (Marquardt 1988, 1989), is discussed elsewhere (Kienle and Marquardt, 1990). Acknowledgements-The authors would like to thank Prof. Dr.-Ing. E.D. Gilles for uutiation of these investigations and

his encouragement. The financial support of DFG (Deutsche Forschungsgemeinschaft) and NATO (for W. Marquardt) is apprcciatcd.

NOTATION

interfacial area per unit length, m’/m ratio of the molar flow rates in both phases (A = L/V)

Aij 4 cij

A J hi

ki L 1 N P Pis P R T V X

Y Z

interaction parameter in the Two-Suffix Margules equation, J/mol dimensionless mass transfer coefficient parameters of the Antoine equation of component i equilibrium concentration of the vapour phase Jacobian of the global material balances dimensionless mass transfer rate of component i mass transfer coefficient, mol/s/m2 molar flow rate of the liquid phase, mol/s length, m number of components common parameter vapour pressure of the pure component i, Pa pressure, Pa universal gas constant, Jjmol/K temperature, K molar flow rate of the vapour phase, mol/s vector of the mole fractions in the liquid phase vector of the mole fractions in the vapour phase dimensionless spatial coordinate

Greek letters relative volatility sections of the axes under the Newton %2 Pi %

gam sij

3;: Y

oi

Kronecker

symbol

activity coefficient eigenvalue of the equilibrium relation Newton number eigenvalue of the linearized model

Subscripts lower end of the column upper end of the column + s singular point species i, i, k

section section

dia-

Bifurcation

analysis

and steady-state

REFERENCES

Arls, R., 1975, The Mathematical Theory of D%fusion and Reaction in Permeable Catalysts, Vols 1 and 2. Clarendon Press, Oxford. Arnold, V. I., 1981, Singularity Theory. Cambridge

Univerisity Press, Cambridge. Chavez, R. C., Seader, J. D. and Wayburn, T. L., 1986, Multiple steady-state solutions for interlinked separation systems. Ind. Enyny Chem. Fundam. 25, 566576. Doherty, M. F., 1985, Properties of liquid-vapour composition surfaces for multicomponent mixtures with constant latent heat. Chem. Enqng Sci. 40, 1979-1980. Doherty, M. F. and Perkins, J. D., 1978a, On the dynamics of distillation processes. I. The simple distillation of multicomponent, nonreacting, homogeneous liquid mixtures. Chem. Engng Sci. 33,281-301. Doherty, M. F. and Perkins, .I. D., 1978b, On the dynamics of distillation processes. II. The simple distillation of model solutions. Chem. Engng Sci. 33, 569-578. Doherty, M. F. and Perkins, J. D., 1979, On the dynamics of distillation processes. III. The topological structure of ternary residue curve maps. Chem. Engng Sci. 34, 1401-1414. Doherty, M. F. and Perkins, I. D., 1982, On the dynamics of distillation processes. IV. Uniqueness and stability of the steady-state in homogeneous continuous distillations. Chem. Engng Sci. 37, 381-392. Doherty, M. F., van Dongen, D. B. and Haight, J. R., 1984. Response to comments on “Material stability of multicomponent mixtures and the multiplicity of solutions to phase-equilibrium equations. 1. Nonreacting mixtures.” Ind. Enyng Chem. Fundam. 23, 374-377. Eigenberger, G., 1972, On the dynamic behaviour of the calalytic fixed-bed reactor in the region ofmultiple steadystates-I. The influcncc of heat conduction in two phase models. Chem. Engng Sci. 27, 1909-1915. Golubitsky, M. and Schafer, D. G., 1985, Singltlarities and Groups in Bifurcation Theory, Vol. 1. Springer, New York. Henley, E. J. and Seader, J. D., 1981, Equilibrium Stage Separation Operations in Chemical Engineering. John Wiley, New York. Kienle, A., 1989, Mehrdeutigkeiten bei Stofftrennprozessen. Diploma&it, Institut fiir Systemdynamik und Regelungstechnik, Universit?it Stuttgart (unpublished, in German). Kienle, A. and Marquardt, W., 1990, On travelling waves in multicomponent, counter-current separation processes. Paper 239d Annual Meeting of AIChE. Chicago, 11-16 November 1990. Kouchnirenko, A. G., 1976, PolyCdres de Newton et nombres de Milnor. Inuentiones Math. 32, l-31 (in French). Krishna, R. and Standart, G. L., 1979, Mass and energy transfer in multicomponent systems. Chem. Engng Commun. 3, 201-275. Kriiner, A., Hall, P., Marquardt, W. and Gilles, E. D., 1990, DIVA-an open architecture for dynamic simulation. Comput. them. Engng 14, 1289-1295. Kvaalen, E.. Neel, L. and Tondeur, D., 1985, Directions of quasi-static mass and energy transfer between phases in multicomponent open systems. Chem. Engng Sci. 40, 1191-1204. Mecklenburg, J. C. and Hartland, S., 1975, The Theory of Backmixing. John Wiley, London. Prokopakis, G. J. and Seider, W. D., 1981, Azeotropic distillation towers with two liquid phases, in Foundations of Computer Aided Chemical Process Design (Edited by R. S. H. Mah and W. D. Seider). AIChE. Lucia, A., 1986, Uniqueness of solutions to single-stage isobaric flash processes involving homogeneous mixtures. A.1.Ch.E. J. 32, 1761-1770. Lyberatos, G., Kuszta, B. and Bailey, J. E., 1984, Steady-state . ._ multiplicity and bifurcation analysis via tne Newton polyhedron approach. Chem. Engng Sci. 39, 947-960. _ ._ . . . .nrm Magnussen, T., Michelsen, M. L. and hreaensluna, A., IY /y,

1769

multiplicity

Azeotropic distillation usmg UNIFAC. Int. them. Engny Symp. Ser. 56, 42/14.2/19. Malesinski, W.. 1965. Azeotropy and Other Theoretical Problems OJ” Vapour-Liquid Equilibrium. Interscience, New York. Marquardt, W., 1988, Nichtlineare Wellenausbreitung-in Weg zu reduzierten dynamischen Modellen von Stofftrennprozessen. VDI Fortschrittsberichte, Reihe 8, Nr. 161, Diisseldod (in German). Marquardt, W., 1989, Wellenausbreitung in verfahrenstechnischen Prozessen. Chemie-Ingr-Tech. 61, 362-377. English translation in Inr. them. Engng 30, 585-606. Michelsen, M. L., 1984, Comments on “Material stability of multicomponent mixtures and the multiplicity of solutions to phase-equilibrium equations. 1. Nonreacting mixtures.” Ind. Engng Chem. Fundum. 23, 373-374. Rhee, H. and Amundson, N. R., 1974, Shock layer in two solute chromatography: effect of axial dispersion and mass transfer. Ckem. Engng Sci. 29, 2049-2060. Serwinsi, M. and Gorak, A., 1983, Eine neue Berechnungsmethode fiir den axialen Konzentrationsverlauf bei der Mehrstofli-ektifikation in FiillkGrperkolonnen. Teil I, II, III. Verfuhrenstecknik 17, 44C-444, 469473, 539-545. Sridhar, L. N. and Lucia, A., 1989, Analysis and algorithms for multistage separation processes. Ind. Engng Chem. Res. 28, 793-803. Sridhar, L. N., Lucia, A. and Guo, X. I., 1989, Analysis of multicomponent, multistage separation processes: fixed temperature and pressure profiles. Paper 2Of, Annual Meeting of AIChE, San Francisco, 5-10 November 1989. van Dongen, D. B., Doherty, M. F. and Haight, J. R., 1983, Material stability of multicomponent mixtures and the multiplicity of solutions to phase-equilibrium equat$ons. 1. Nonreacting mixtures. Ind. Engnq Ckem. Fundam. 22, 472485. Witmer, G. S., Balakotaiah, V. and Luss, D., 1986, Multiplicity features of distributed systems--I. LangmuirHinshelwood reaction in a porous catalyst. Chem. Engng Sci. 41, 179-186.

APPENDIX

A: PARTIAL

DERIVATIVES

OF THE

EQUILIBRIUM RELATION BASED ON THE TWO-SUFFIX MARGLJLES MODEL FOR THE ACTIVITY COEFFICIENTS OF THE LIQUID PHASE Differentiation of eq. (30) with respect fractions xi yields

wi 1 dp;, dT -_---yX,+!5 axj

pdTL$

to the liquid

mole

*

i.j = l(1)N - 1 (Al) where 6,, 1s the Kronecker symbol. The required derivatives of the activity coefficients and the vapour pressure of the pure components can be obtained directly by differentiation of cq. (34) and eq. (33), respectively:

dpi,

ci2

dT = pi”(T dyi

-

dxj

=

FT $ (Aaj + k

1

‘i z-

+ ci#

A,k - Aa - Ai~)xt

- -Ly.lny T’

(A3)

’

The partial derivatives of the boiling temperature with respect to the liquid composition, have to be calculated from eq. (31) by implicit differentiaton: dT _=_ axj

(A5)