Qualitative Analysis for Homogeneous Azeotropic ... - ACS Publications

Ind. Eng. Chem. Res. 2002, 41, 3943-3962

3943

Qualitative Analysis for Homogeneous Azeotropic Distillation. 2. Bifurcation Analysis Cornelius Dorn† and Manfred Morari* Automatic Control Laboratory, Swiss Federal Institute of Technology, ETHZ, CH-8092 Zu¨ rich, Switzerland

Studying global mass balances and how they interact with the shape of column profiles, Dorn and Morari (Dorn, C.; Morari, M. Ind. Eng. Chem. Res. 2002, 41, XXXX) developed a graphical method for a qualitative stability analysis of steady-state profiles in homogeneous azeotropic distillation columns. Based on the qualitative stability analysis, here it is shown how bifurcation diagrams can be constructed containing information on the stability of the equilibria. Fold, Hopf, and homoclinic bifurcations as well as more complex bifurcations can be predicted. The existence of limit cycles is explained and their approximate shape determined. It is explained how the feed composition and the reflux flow rate influence the bifurcation behavior of the column. The results are compiled in two-parameter bifurcation diagrams as well as in projections of a threeparameter bifurcation diagram. The method is applied to predict several new phenomena in distillation, e.g., operating conditions for which the only steady state of the column is unstable and surrounded by a stable limit cycle. All predicted phenomena are confirmed by steady-state and dynamic simulations. 1. Introduction Because their understanding is a necessary prerequisite for proper column design and operation, the steady state and dynamic behavior of azeotropic distillation has been studied extensively over the past decades. Studying the steady-state behavior of distillation columns assuming infinite column length and infinite reflux2,3 allowed new insights from analysis (e.g., multiple steady states) leading to practical implications for column design and operation.4 In the same way, the theory presented in this series of papers can form a basis for analyzing the dynamic behavior of distillation columns. New insights so far are the explanation of phenomena observed earlier (e.g., limit cycles) as well as the prediction of new phenomena (e.g., bifurcations). In a next step leading beyond the scope of this work, the better understanding of column dynamics should again aid column design and operation. In the preceding paper1 of this series, a qualitative method to analyze the stability of steady-state profiles in homogeneous azeotropic distillation columns was presented. In this paper, the qualitative stability analysis is extended to construct complete bifurcation diagrams and locate dynamic phenomena like limit cycles. As discussed in detail in the preceding paper, dynamically evolving column profiles at infinite reflux coincide with a section of a distillation line at any time. The exchange of holdup between stages along the column occurs infinitely fast. Thus, if a column operated at infinite reflux is initially or after a perturbation filled with holdups Mkxˆ k(0), k ) 1, ..., n, that do not correspond to a steady-state profile, the initial total column holdup n Mkxˆ k(0) will be redistributed along the column ∑k)1 infinitely fast. The total column holdup composition (TCHC) will continuously change at a slow time scale * Author to whom correspondence should be addressed. Phone: +41 1 632-2271. Fax: +41 1 632-1211. E-mail: [email protected]. † Current affiliation: McKinsey & Company, Switzerland.

n Mk ) FxF - DxD - BxB as governed by (d/dt) xavg∑k)1 until the steady-state global mass balance is satisfied. These results approximately carry over to the finite case if the internal flows are high compared to the column holdups and external flows. In this work, the columns (mainly assuming infinite reflux except when explicitly stated otherwise) are supposed to be long enough such that the product compositions will approximately be restricted in a way similar to that in the ∞/∞ case. The result is that the product paths of those finite columns can be approximated with the ∞/∞ product paths. Differences must be expected where the product paths go through singular points. These differences will be highlighted where appropriate in the later sections. The key idea of the qualitative stability analysis is to study how changes of the column inventory and changes of the shape of the column profile interact. The neighborhoods of the steady-state product compositions can each be divided into six regions as illustrated in Figure 1. In the preceding paper,1 the characteristics of the six types of neighboring regions were analyzed and summarized in Table 1. Using a detailed classification of feasible column profiles (Tables 2 and 3), the stability of each type could be determined (Table 4) by constructing schematic vector fields based on the six neighboring regions for steady-state product compositions shown in Figure 1 in this paper. The overall organization of this paper is as follows: In section 2, the qualitative stability analysis will be used to construct qualitative bifurcation diagrams including bifurcation points and information on the dynamic behavior of the column. In section 3, similarities and differences of the steady-state behavior of columns of finite and infinite lengths, respectively, are discussed because this is crucial for the following discussions. The influence of the feed composition and the reflux flow rate on the bifurcation behavior is studied in detail in section 4. It will be discussed how the standard one-

10.1021/ie010726j CCC: $22.00 © 2002 American Chemical Society Published on Web 07/12/2002

3944

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

Figure 1. Steady-state column profile and regions qualitatively categorizing the impact of eventual non-steady-state product compositions caused by perturbations of the profile.

parameter (D) bifurcation diagrams change under variation of these parameters. It will be shown how some qualitative predictions, although less rigorous when compared to the rest of the results presented here, are still possible even when the assumption of infinite reflux is dropped. The results will be generalized in section 5. Different types of two-parameter bifurcation diagrams will be introduced as a useful tool. Finally, a certain projection of a three-parameter bifurcation diagram will be shown in section 6. This allows one to identify the most influential parameters. In section 7, a short step summary is given from an application point of view. Numerical case studies confirming the predictions are presented for two ternary mixtures in section 8. 1.1. Earlier Findings of Complex Bifurcations. Bifurcations of codimension onesat least the more common ones such as fold bifurcation, Hopf bifurcation, and homoclinic bifurcationsare frequently encountered in physical and technical systems, and numerous references to their occurrence exist in the literature. (Literature overviews for multiple steady states in distillation columns are given by Bekiaris et al.,5 for oscillatory phenomena including Hopf bifurcations by Lee et al.,6 and for instability in distillation in the preceding paper.1) In contrast, only more recently mathematical theory has led to a better understanding of bifurcations of codimension two. Normal forms for all generic bifurcations of codimension two were derived, and their local behaviors were analyzed. Kuznetsov7 gives a very detailed overview of all of the bifurcations mentioned in this paper; a streamlined summary can be found in the Ph.D. Thesis of Dorn.8 Cusp bifurcations often occur when the existence of multiple steady states in a system depends on a second bifurcation parameter. This behavior has frequently been observed in processes although it was rarely attributed to a cusp bifurcation. The same holds for one of the subcases of the degenerate Hopf bifurcation and the degenerate homoclinic bifurcation. The subcase

considered here relates to the Hopf and homoclinic bifurcations, respectively, in the same way as the cusp relates to the fold bifurcation. Their appearance in physical and technical systems has been mentioned in the literature but rarely given the appropriate name. The Bogdanov-Takens bifurcation9,10 is more interesting. In the last 2 decades, it has been found in the models of many natural and technical processes: shock waves,11 convection12 and boundary layer flows,13 certain lasers14 and nonlinear optical elements,15 motion of a thin panel in a flow16 as well as dynamics of a windloaded nonlinear structure,17 and solar gravity.18 Also, they have been found in autonomous electrical and mechanical systems that can be described as Rayleighvan der Pol-Duffing oscillators.19,20 Further, BogdanovTakens bifurcations appear in certain models of population dynamics21 and even in models of romantic relationships.22 In chemical engineering, the authors are aware of Bogdanov-Takens bifurcations in reactiondiffusion systems23 and in an extended Brusselator.24 1.2. Overview of New Phenomena. A couple of phenomena in homogeneous azeotropic distillation will be discussed in this work that are explained or some of them even observed for the first time. These are Hopf bifurcations and limit cycles, which were first reported by Lee et al.6 Further, homoclinic bifurcations (i.e., two unstable steady states and one stable steady state that is the only attractor) and Bogdanov-Takens bifurcations (resulting in unstable steady states on both sides of a turning point) are explained. They were first observed by Dorn et al.25 In addition, a situation where a stable limit cycle is the only attractor is predicted and confirmed in numerical simulations. 2. Qualitative Bifurcation Analysis 2.1. Step 1: Continuation of Equilibria. In a first step, an ∞/∞ analysis3 is carried out. The only necessary information are the feed composition and the characteristics of the residue curve map, i.e., the location and type of azeotropes. (For mixtures belonging to classes more complex than the 001 case studied in this paper, also the location of the boundaries may be required.26) From the product paths derived from the continuation of all feasible profiles, a bifurcation diagram can be constructed using the lever rule. For the product paths of the ∞/∞ case shown in Figure 2, a corresponding bifurcation diagram is shown in Figure 3. That diagram contains information on the stability of the equilibria that is determined in the next step. 2.2. Step 2: Determining Stability. For the bifurcation diagram shown in Figure 3, the stability of the equilibria can be determined by application of the qualitative stability analysis. The product path consists of a finite number of parts each of which corresponds to a continuum of profiles of the same type and subtype (Figure 2). When the qualitative stability analysis is applied on a prototype of each of the types of profiles present in the actual case, the bifurcation diagram can be augmented with stability information. In Figure 2, there are profiles of type II along the product path of both subtypes: stable and unstable. Using the lever rule to construct the corresponding bifurcation diagram shown in Figure 3, it can be observed that there are unstable steady states not only on the middle branch but also on the part of the high branch that corresponds to profiles of type IIu.

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3945

Figure 4. Illustration of a bifurcation diagram including fictitious branches of amplitude for limit cycles. D1: subcritical Hopf bifurcation. D2: supercritical Hopf bifurcation. D0: fold bifurcation of limit cycles.

Figure 2. Product path of the ∞/∞ column plotted into a residue curve map.

Figure 3. Illustration of a bifurcation diagram constructed from the product paths shown in Figure 2.

2.3. Step 3: Determining Bifurcations. Having constructed a bifurcation diagram as sketched in Figure 3, the next step is to determine what bifurcations are present. Fold bifurcations (turning points) are immediately obvious in bifurcation diagrams but can also be located geometrically in the product paths (cf. appendix A). Other bifurcations, for example, Hopf bifurcations, can occur when there is a local exchange of stability of the equilibria. There are other (global) bifurcations (e.g., homoclinic bifurcations) that can only be detected by studying the global behavior of the system. Hopf Bifurcation. If two branches of equilibria meet in a turning point as shown in Figure 3, one or both can be unstable, but it is not possible that both are stable. In contrast, the interval of unstable equilibria between D1 and D2 on the otherwise stable high branch is caused by more complex effects. A Hopf bifurcation must be present at D1 and D2. It is necessary to determine whether these Hopf bifurcations are subcritical or supercritical. The supercritical type of the Hopf bifurcation results in stable limit cycles surrounding the unstable equilibria in the neighborhood of the bifurcation point, while unstable limit cycles surround the stable equilibria in the case of a subcritical Hopf bifurcation.

This question can be answered from studying the residue curve map and applying the tools of the qualitative analysis of column profiles. The approach is the following: first, the steady-state profile of an equilibrium on the high branch for D1 < D < D2 is studied. It is unstable and, as will be shown, can be surrounded by a stable limit cycle. Then it will be illustrated how the amplitude of the limit cycle goes to zero as D approaches D1 or D2 from inside the interval (D1, D2). This allows the conclusion that two supercritical Hopf bifurcation points are present in D1 and D2. Note that it is necessary that the amplitude of the stable limit cycle goes to zero at both D1 and D2. Otherwise, one or both of the Hopf bifurcations points could also be subcritical, each occurring together with a fold bifurcation of limit cycles. In Figure 4, this is illustrated for the left Hopf bifurcation point being subcritical. This scenario results in a hard generation of limit cycles. Increasing D past D1 results in a sudden existence of stable limit cycles of finite amplitude. If the system is undergoing sustained oscillation, reducing the distillate flow rate below D0 causes an out-of-the-blue disappearance of the limit cycles. Consistent with the analysis to follow in the next paragraphs, this has (to the authors’ knowledge) never been observed in any published simulation or experiment of homogeneous azeotropic distillation columns. In the case of a homoclinic bifurcation, limit cycles of finite amplitude also suddenly appear or disappear. However, the period of the limit cycles becomes infinitely long as the homoclinic bifurcation point is approached. In contrast, the period of the limit cycles remains finite when a fold bifurcation of limit cycles is approached. Therefore, homoclinic bifurcations are not denoted “out-of-the-blue”. To determine the stability of a steady-state column profile, only a small neighborhood around it has to be studied. For an analysis of limit cycles, a larger neighborhood around the steady-state profile has to be examined. In Figure 5a, the neighborhood of the distillate composition of a steady-state profile of type IIu is shown in the composition space together with residue curves. The profile itself is not shown for clarity. The regions 1-6 (cf. section 1) are overlayed with a region Ω (as defined below) drawn in gray to emphasize an important observation: In all regions except region 2, the perturbed column profiles can be not only of type IIu but also of type IIs. This depends on how far away the perturbed profile is from the steady state. Profiles beginning in the direct neighborhood of the steady-state distillate composition are of type IIu.

3946


Figure 5. Residue curve diagram illustrating the neighborhood of a steady-state distillate composition of a profile of type IIu: (a) classification of regions; (b) approximate vector field for the distillate composition.

Profiles beginning further away from the steady-state distillate composition correspond to type IIs. Therefore, regions 1 and 3-6 each consist of a “stable” and an “unstable” part, while region 2 is entirely unstable. In the following, the unstable parts will be denoted as regions 1u-6u while the stable part will be denoted regions 1s and 3s-6s. Note that the unstable part of region 3 is further subdivided. For distillate compositions in region 3u getting close to pure L, the finite number of trays in the column will ultimately limit the separation power in such a way that the distillate composition moves in the direction tangent to the residue curve (cf. the discussion of profiles of type IIu for details1). That subregion will be denoted as 3v. For the union of the unstable parts of all six regions, we define the following: Definition 1 (Region Ω). The subset of the composition space that corresponds to the distillate compositions of all feasible profiles of type IIu is called region Ω. It is bounded by the curve δΩ. In the case of infinite reflux, the nonmonotonicity condition (cf. section 1) for profiles of type II to be unstable can be directly related to the geometry of the residue curves. In this case, the boundary δΩ of region Ω is defined as the set of all nontrivial maximum points of the residue curves in terms of L: δxL(H)/δH ) 0 where xL(H) are local parametrizations of the solutions to the equations defining residue curves. The regions 1 and 3-6 are separated into their respective stable and unstable parts by δΩ. As discussed before, the qualitative stability analysis allows one to determine qualitatively for each region the set of directions in which the distillate composition of a perturbed column profile will move during its transient. In Figure 5b, these directions are illustrated by arrows. The arrows in close neighborhood of the steady state (regions 1u-6u) correspond to the situation of type IIu profiles that were studied in the preceding paper1 and found to be unstable. In the same way, the arrows in the stable subregions correspond to the situations analyzed for profiles of type IIs. The arrows illustrated in Figure 5b can be interpreted as a vector field for xD, which was constructed using the arguments of the qualitative stability analysis. Next, it will be studied if the vector field allows the existence of a stable limit cycle. Assume (without loss

of generality) that the steady-state profile is subject to a perturbation that causes the distillate composition to move away slightly into region 1u. Hence (cf. Table 1 in the preceding paper1), I will subsequently accumulate in the profile. This shifts the profile onto an outer residue curve (i.e., closer to the boundary between AzLH-L-I-H), causing the distillate composition to move into region 2u and then region 3u, while moving further away from the steady-state composition and coming nearer the boundary between pure L and the L-H azeotrope. In region 3u, the profile continues its movement onto outer residue curves until finally the distillate composition enters region 3v as the separation power of the column becomes a limiting factor. Hence, the distillate composition will begin to move away from the L-H boundary and approach the boundary between regions 3 and 4. When entering region 4 coming from region 3, the distillate composition could be either in region 4u or in region 4s (by way of region 3s). In both cases, xD is going to come closer to the boundary δΩ between regions 4s and 4u. Also, xD will move toward the line separating regions 4 and 5, finally entering either region 5s or region 5u. If the distillate composition is in region 5u, the column profile will change its shape such that xD will move away from the steady-state composition, toward the line separating regions 5 and 6, and toward the curve δΩ separating regions 5u and 5s. To understand the way the distillate composition is going to move while in region 5s, it is illustrative to study two limiting cases. If the column profile is of type IIs with xD close to the line separating regions 4s and 5s, then the TCHC is mainly going to decrease in I and to increase in H while not changing much in L. In contrast, if the distillate composition is in region 5s close to the line separating regions 5s and 6s, the TCHC is going to increase in L and to decrease in I while not changing much in H. In both cases, the distillate composition will therefore approach the boundary δΩ between regions 5s and 5u and finally leave region 5 toward region 6. At this point, the question might arise whether xD can, for example, be driven to increase in L while approaching the boundary between regions 5s and 5u, then entering region 5u, there decreasing in L while being driven, entering region 5s again, and thus possibly forming a limit cycle inside region 5, perhaps including


Figure 6. Approximate vector field for the distillate composition in the neighborhood of a steady state of type IIu and the surrounding limit cycle.

also region 4? A quick glance at Figure 5b might suggest this scenario. However, this is impossible because there are no equilibria around which a limit cycle could form. Instead, once on δΩ, xD can be assumed to slide along the boundary until entering region 6. The same arguments that were used in region 5 can also be applied to discuss the movement of the distillate composition in region 6. Finally, the distillate composition will enter region 1 again. The distillate composition will now start its movement in region 1, but in contrast to the starting point of this discussion, it will start further away from the steady-state composition. Therefore, the distillate composition will initially spirally diverge from the steady state. However, that divergence is bounded. In regions 1u, 2u, and 3u, the distillate composition will approach the boundary between pure L and the L-H azeotrope but obviously cannot cross it. In regions 1u and 3u-6u, the movement away from the steady-state composition is bounded by the neighboring regions of the stable type s in which the distillate composition will drift back in the direction toward the steady state. Therefore, the existence of a stable attracting cyclic orbit (limit cycle) can be concluded. The limit cycle has the approximate shape illustrated in Figure 6. As was discussed before, the existence of stable limit cycles is not sufficient for the Hopf bifurcations to be supercritical. Therefore, the amplitude of the limit cycles will be studied next. The boundary δΩ of region Ω played an important role for the location of the limit cycles of xD in the composition space. In regions 4, 5, and 6, the limit cycle and δΩ were found to be closely related, approximately coinciding (cf. Figure 6). Approaching the two Hopf bifurcation points on the product path (under variation of the bifurcation parameter D′), the steady-state distillate composition of the unstable equilibria approaches δΩ by definition. The amplitude of the limit cycle will therefore get smaller. When the situation in Figure 6 is studied for xD on δΩ, the amplitude of the limit cycles is zero. Hence, the Hopf bifurcations are supercritical. Focusing on the distillate composition was useful in the above discussion because it is the main source of changes in the TCHC. However, because the distillate composition is basically just one end of the column profile, it is obvious that not only xD is going to oscillate

on the limit cycle but all compositions on all trays. However, the amplitude of the oscillation will differ significantly. The bottoms composition, for example, was assumed to coincide with pure H as a pinch point at all times at the very beginning of this discussion. This was justified as a useful assumption for the qualitative stability analysis, but, of course, xB will never stay exactly on pure H. Studying dynamic simulations of distillation columns exhibiting sustained oscillations reveals that the bottoms composition is in fact oscillating, but with an extremely small amplitude. If that amplitude grows above negligible values, a new bifurcation can occur as studied in the following. Homoclinic Bifurcation. As the behavior of profiles of type II is studied in this section 2, the bottoms composition was assumed to remain at approximately pure H at all times. Situations in which the bottoms composition moves away from pure H more than just slightly behave very differently. They can occur when the distillate composition enters region 1 (coming from region 6) relatively far away from the steady state. As a result, the total change of the TCHC will be bigger because it lasts longer until the distillate composition leaves region 1 again toward region 2. The impact of this change in the TCHC (accumulation of L and I and depletion of H) on a profile of type IIu was discussed in the course of the stability analysis:1 The accumulation of I causes the profile to shift onto an outer residue curve, i.e., closer to the boundary from the L-H azeotrope over pure L and pure I to pure H. It was found that the depletion of H in the column causes trays between pure H and pure I to move toward pure I, with the trays closest to the pinch point I passing it and already changing direction toward pure L, all along the residue curve, a section of which the column profile must coincide with. As a result of the accumulation of L in the TCHC, trays between pure I and pure L move toward pure L. For that reason and because the column cannot get “longer”, i.e., the total separation power in the column remains constant, the trays between pure L and the L-H azeotrope also move toward L. If those changes are too big, there will be no more trays at the bottom of the column having a composition of approximately H because more and more have already started to move away from that pinch point before. Thus, the pinch point in H does no longer exist and the bottoms composition will start moving toward pure I. The perturbation of the bottoms composition causes a depletion of I and an accumulation of H in the TCH, partially compensating the opposite changes resulting from the perturbation of the distillate composition. From here on, two different scenarios can govern the ongoing transient. Scenario 1. If the change of I and H in the TCHC through the distillate continues to dominate over the change through the bottoms, the transient can continue with xD moving into region 2 and quickly into region 3. Then, at the latest, the accumulation of H in the column will start to bring the bottoms composition back to pure H. In the following transient, xD moves into regions 4, 5, and so on. In conclusion, this case is a limit cycle as described before, but this time the bottoms composition is oscillating more significantly. Scenario 2. If, on the other hand, the change of I and H in the TCHC is dominated by the perturbation of the bottoms composition, there will be an accumula-

3948


tion of L and H and a depletion of I. The depletion of I would cause the profile to shift onto an inner residue curve, and as a result the distillate composition will eventually start to move in the opposite direction, i.e., toward the L-H azeotrope. The increasing accumulation of L in the column will make it most likely that the distillate composition will approach the L-H azeotrope. In that case, the column has converged toward a type I profile. Having found a situation where the unstable steadystate profile of type IIu can or cannot be surrounded by a stable limit cycle, it remains to determine what mechanism causes the switch between these two situations. As mentioned at the beginning, the latter situation can occur when the distillate composition enters region 1 relatively far away from the steady state. In other words, if the size of the limit cycle grows above a certain size, it ceases to exist because it violates the condition of the bottoms composition remaining at pure H. More precisely, the profile violates the condition of xB remaining at pure H too much. This indicates the occurrence of a homoclinic bifurcation, in which a limit cycle increases in amplitude until touching a saddle equilibrium and then disappearing. Without referring to any rigorous norm here, the amplitude of the limit cycle can be seen as the size of the corresponding closed orbit in the composition space. From the arguments above describing the possible existence of limit cycles, it follows that the size of the limit cycle depends on the distance of the steady-state distillate composition from the boundary of region Ω. Note that the arguments used above cannot give a quantitative prediction of the operating conditions for which the homoclinic bifurcation will occur. This is in contrast to the Hopf bifurcation, which can be predicted quantitatively as the intersection of the distillate product path with the boundary of region Ω. What is known is that the homoclinic bifurcation points must lie between Hopf bifurcation points. Oscillation Region. A phenomenon that is of interest in many applications is the occurrence of underdamped oscillations in the open-loop response of a system. Because supercritical Hopf bifurcations were found to be present in distillation columns for certain conditions, there will also be steady states where the eigenvalues with the largest real parts are a pair of complex conjugates. Thus, the column will exhibit underdamped oscillations when approaching these steady states after perturbations or changes in the operating conditions. This will hold in the neighborhood on both sides of a Hopf bifurcation point where a stable focus exchanges stability. In the following, a region on a steady-state branch corresponding to foci will be called an oscillation region. Note that an equilibrium changing its type from node to focus does not undergo a bifurcation in the mathematical sense.7 3. Profiles Revisited: Finite Columns As mentioned in section 1, also columns of sufficient finite length will contain at least one pinch point. However, the composition profile will only come near the respective singular points of the residue curve map, but it will in no point coincide exactly with one of them as happens in the ∞/∞ case. Profiles of type I and II contain only one pinch point, and that pinch point is at the end of the profile. In contrast, profiles of type IIIa and IIIc contain one pinch point in their middle, and profiles of type IIIb contain two pinch points. Keeping

Figure 7. Illustration of the product paths for a column of finite and infinite lengths, respectively.

the separation power of the column constant, it is obvious that the pinch composition can come nearer to a singular point in the case of profiles of type I and II. In Figure 7, this is illustrated for the distillate composition. The distillate composition in the finite column cannot reach pure L. Hence, the distillate product path turns near pure L. Depending on the actual conditions (number of trays, position of feed tray, feed composition, and reflux), the distillate composition might not even reach the boundary between L and I. The same holds for the bottoms product path near pure I. These two cases correspond to profiles of type IIIb. Knowing the shape of the product path in the finite case, it is still necessary to classify all of the profiles of type III into the subclasses a-c. One approach could be to differentiate the types by use of the geometric condition for multiple steady states3 (appendix A). After location of the turning points, they could be used to segment the product paths. In this way, the subtype of the profiles of type III can be related to the branch on which the corresponding equilibria are. However, the consequence of this approach would be that the stability properties of the subtypes IIIa-IIIc would not carry over to the finite case under all conditions.8 In the following, a different approach is pursued here. The subclasses for the profiles of type III will be defined through the monotonicity properties of their corresponding column composition profiles. Profiles of type IIIa are nonmonotonic in L and monotonic in I and H. Profiles of type IIIc are monotonic in L. Profiles of type IIIb are nonmonotonic in all three components. 4. Influence of Parameters In section 2 and in most of the literature, the main bifurcation parameter in the diagrams shown is the distillate flow rate. However, other parameters too can have a qualitative impact on the dynamic behavior of the distillation column. In the following, it will be discussed how the standard one-parameter (D) bifurcation diagrams change under variation of other parameters, including reduction of the internal flow rates in the column from infinite to finite values. This is a subset of the general multiparameter bifurcation problem.


First, however, it must be determined what parameters are interesting. For a given feed composition and quality, the constant molar overflow (CMO) model of a homogeneous distillation column with a defined pressure profile has 2 degrees of freedom. Of those, the dimensionless distillate flow rate D′ ) D/F and the dimensionless reflux flow rate L′ ) L/F will be chosen here because they represent the internal and external flows. The composition of the ternary feed allows two specifications. The parametrization of the feed composition suggested in a previous paper25 is extended. For the 001 class mixtures studied here, we choose to parametrize the feed composition by xFLI ) xFL/xFI and xFLH ) xFL/xFH. Because the degenerate cases of feeds are excluded, 0 < xFLI < ∞ and 0 < xFLH < ∞ hold. The Hopf bifurcation (and indirectly the homoclinic bifurcation) discussed in section 2 corresponded to an exchange of the local stability of profiles of type II. Therefore, the parametrization of the feed composition was chosen such that its influence on that part (i.e., type II) of the distillate product path becomes as obvious as possible. In the following subsection, it will be studied how this choice of parametrization results in a certain decoupling of the impacts of xFLI and xFLH on the bifurcation behavior of the column. The influence of the feed composition on the location of the turning points can be determined algebraically. With the ∞/∞ analysis,3 it is predicted that for mixtures of class 001 there will be multiple steady states (i.e., two turning points) for all feed compositions. Using the relations

xL )

xLIxLH xLI + xLH + xLIxLH

xI )

xLH xLI + xLH + xLIxLH

xH )

xLI xLI + xLH + xLIxLH

the location of the left and right turning points D′L1 and D′L2 can be calculated to be at

[

D′L1 )

xLIxLH xLI + xLH + xLIxLH

1 D′L1 xAz L D′L2 ) xLI(1 + xLH) xLI + xLH + xLIxLH

xFLH e xAz LH xFLH g xAz LH

]

(1)

(2)

In the above equations, xAz LH denotes the ratio of light to heavy components in the binary azeotrope between L and H. For all positive xLI and xLH, D′L1(xLI,xLH) is monotonically increasing in both xLI and xLH. Thus, also D′L2(xLI,xLH) is monotonically increasing in both xLI and F Az xLH if xFLH e xAz LH. For xLH g xLH, D′L2(xLI,xLH) is monotonically increasing in xLI and monotonically decreasing in xLH. Yet, D′L2(xLI,xLH) is continuous in both xLI and xLH. The following analysis for xFLI and xFLH is in analogy to the concepts presented in section 2, based on the

methodology for the analysis of infinite reflux in columns of finite length. In contrast, the study of the influence of L′ will require new ideas. 4.1. Feed Composition: xFLI. The impact of changes in xFLI on the location of the turning points is seen from eqs 1 and 2. As xFLI increases, both turning points move toward larger values of D′. Also, the size of the MSS region (D′L2 - D′L1) increases. The part of the distillate product path that coincides with a section of the straight line through pure H and the feed composition corresponds to profiles of type II (cf. Figure 7). If xFLI is increased, that part of the distillate product path will move closer to the boundary between pure L and pure H. In Figure 8a, the part of the distillate product paths corresponding to profiles of type II is illustrated for different values of xFLI. Because this is important to capture some complex bifurcations, the figure shows the type II part of the product paths of a column of finite length. As discussed in section 3, the product paths of a column of finite length cannot go through saddle singular points such as pure L. Hence, for high values of xFLI, the distillate product paths will not reach the boundary between pure L and pure I. In case a1 in Figure 8, all steady states on the part of the high branch corresponding to profiles of type II will be stable because there are no profiles of type IIu. If xFLI is increased, the product path moves closer to the boundary between pure L and pure H and thus also approaches region Ω. As is illustrated in Figure 8a in case a2, for some value of xFLI, when the product path just touches the boundary of region Ω (cf. definition 1), there will be one single steady state on the high branch which is on the boundary between stability and instability. In this point, the two Hopf bifurcation points are born. If xFLI is further increased (case a3), there are two distinct intersections of the product path, with the boundary of Ω each corresponding to a Hopf bifurcation point as discussed in section 2.3. For higher values of xFLI, the two points move away from each other on the high branch. If xFLI is increased even further, the amplitude of the limit cycle can grow enough to finally cause two homoclinic bifurcations to occur. After a further increase of xFLI, the situation shown in case a4 of Figure 8a can arise: there is only one Hopf bifurcation, and all steady states between that Hopf bifurcation point and the turning point are unstable. This was first observed in numerical simulations.25 The mechanisms for going from case a3 to case a4 will be described and classified in section 5. 4.2. Feed Composition: xFLH. The impact of changes in xFLH on the location of the turning points is seen from eqs 1 and 2. Both turning points move toward larger values of D′, and the size of the MSS region (D′L2 - D′L1) F increases as xFLH is increased until xFLH ) xAz LH. As xLH is increased further, the right turning point D′L2 moves toward lower values of D′ while D′L1 is still increasing. Hence, the MSS region becomes narrower again. The influence of the changes in the feed composition when xFLI is held constant can be studied in a way analogous to studying the influence of xFLI. However, because the influence of xFLH partially depends on xFLI, they must be studied together. In Figure 8, the part of the distillate product paths corresponding to profiles of

3950


Figure 8. Illustration of the impact of xFLI and xFLH on the bifurcation behavior.

Figure 9. Illustration of the impact of xFLH on the distillate product path in a column of finite length.

type II is illustrated for different values of xFLI and xFLH. Figure 8b shows that the influence of xFLH is basically to determine where the distillate product path ends, while Figure 8a shows that changing xFLI has an impact on the shape of the product paths and thus on what bifurcations can occur. Again, it is important to note that in Figure 8 the part of the distillate product path corresponding to profiles of type II does not reach the boundary between pure L and pure I if xFLI is high. When cases a1 and b1 in Figure 8 are compared, it is obvious that there will be no qualitative difference between the corresponding bifurcation diagrams. However, for higher values of xFLI, increasing xFLH can have the effect illustrated in cases a3 and b3. As the feed composition moves into region Ω, one of the Hopf bifurcation points disappears. (Technically speaking, it moves out of the physically meaningful interval of distillate flow rates 0 e D′ e 1.) This allows an interesting prediction: if the feed composition is inside region Ω, the steady states on the high branch are stable between the left turning point and the Hopf bifurcation point and are unstable between the Hopf bifurcation point and the upper limit D′ ) 1 of the parameter interval. In this case, the steady states on the high branch are unstable and surrounded by a stable limit cycle for D′L1 e D′ e 1. This means that operating conditions exist for which a homogeneous azeotropic

distillation column has no stable steady state, with a stable limit cycle being the only attractor. To the authors’ knowledge, this phenomenon has not been observed or predicted before for homogeneous azeotropic distillation columns. It is illustrated with dynamic simulations in section 8.2. From the geometrical properties of region Ω and the distillate product path, it is obvious that this phenomenon can only occur for xFLH > xAz LH. The intersection of the distillate product path with δΩ allows a quantitative prediction of the occurrence of this phenomenon. In contrast, the phenomena described in the following are predicted in a more qualitative way. Increasing xFLH beyond certain values will result in additional effects caused by the finite number of stages. Parts a and b of Figure 9 illustrate that the maximum achievable purity in the distillate (max xD L ) depends D only a little on the feed composition. This only holds as long as the feed composition is relatively far away from the distillate composition corresponding to max xD L. D F However, the closer x comes to that distillate composition (i.e., xFLH increases), the more also max xD L has to D increase. This is illustrated in Figure 9c where increasing xFLH finally caused the resulting distillate product

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3951 Table 1. Legend for One- and Two-Parameter Bifurcation Diagrams branches bifurcation points

bifurcation lines

s -f O + CP DH DP BT T H P

stable unstable oscillation (nonbif) Hopf homoclinic cusp degenerate Hopf degenerate homoclinic Bogdanov-Takens fold Hopf homoclinic

path to be “pushed” out of region Ω again. Therefore, steady states on the high branch next to the left turning point are stable again. A rigorous analysis of this case would be very complex because the assumptions on the bottoms product path do not hold anymore: the bottoms composition for the profiles corresponding to the distillate composition close to the upper turning point will not be on the boundary between pure I and H anymore. This effect will therefore not be studied in more detail in the following. As a conclusion for the influence of xFLH on the bifurcation behavior in this case, it can be stated that increasing xFLH above the point in which the maximum achievable purity in L remains relatively unchanged causes additional bifurcations. 4.3. Reflux. Studying the influence of the reflux flow rate on the bifurcation behavior is difficult because it implies dropping of the assumption of infinite reflux. Two aspects have to be considered: How do the product paths and the shape of profiles change as the reflux flow rate is decreased? As the reflux is decreased, the product paths (Figure 10a) move away from the boundaries of the composition space (cf. section 3). In addition, the shape of the profiles (Figure 10b) can change such that a profile of type IIu (unstable) may turn to type IIs (stable). This results in a reduction of the size of region Ω. Hence, Hopf and other bifurcations disappear once reducing L′ caused region Ω to shrink so much that its boundary δΩ does not intersect with the distillation product path anymore. The mechanisms are discussed in more detail in the Ph.D Thesis.8

5. Bifurcation Diagrams In the previous section, it was studied how the parameters xFLI, xFLH, and L′ influence the bifurcation behavior of a distillation column separating a 001 class mixture. In this section, these results will be generalized and presented in two-parameter bifurcation diagrams. Using the theory presented in section 4, the location of bifurcation points can be predicted in dependence of two parameters. The two-parameter bifurcation diagrams presented here are qualitatively correct for all 001 class mixtures. Two-parameter diagrams can also be determined numerically for given mixtures. The advantage of two-parameter bifurcation diagrams is their compactness in which they illustrate regions of different types of dynamic behavior in a two-dimensional parameter space. If it is determined which type of bifurcation behavior is present for given operating conditions, studying two-parameter bifurcation diagrams allows an instantaneous prediction of how the bifurcation behavior is going to change under variation of the parameters. Another advantage of two-parameter bifurcation diagrams is that they highlight the role of certain bifurcation points of codimension two as organizing centers for the boundaries of these regions in the parameter space. A legend explaining the symbols used in the following illustrations is provided in Table 1. 5.1. Two-Parameter (xFLI, D′) Bifurcation Diagram. For a given xFLH and infinite (or sufficiently high) reflux flow rate, the influence of xFLI on the bifurcation behavior is illustrated in the two-parameter (xFLI, D′) bifurcation diagram in Figure 11. As discussed in section 4.2, xFLH was chosen to be sufficiently small (xFLH < xAz LH) such that the feed composition is outside region Ω (cf. Figure 8b). The two-parameter (xFLI, D′) bifurcation diagram is illustrated with a series of standard one-parameter (D′) bifurcation diagrams. Each of the one-parameter (D′) bifurcation diagrams corresponds to fixing xFLI at a selected value. For case g, the detailed correspondence between the two-parameter (xFLI, D′) bifurcation diagram and a selected one-parameter (D′) bifurcation diagram is illustrated as an example in Figure 12. Because the low and middle branches are less interest-

Figure 10. Illustration of the impact of the reflux flow rate on the bifurcation behavior.

3952


Figure 11. Two-parameter (xFLI, D′) bifurcation diagram with bifurcations of codimension one and two (legend in Table 1). Further, the regions of oscillation are illustrated (gray shaded). Note that the low branch and a part of the middle branch are omitted for clarity.

steady states on the low and middle branches are always stable and unstable, respectively. Thus, the phenomena illustrated in Figure 11 refer to the high branch of the corresponding one-parameter (D′) bifurcation diagrams. The only exception is the right curve T of fold bifurcation points because these correspond to the fold bifurcations connecting the low and middle branches. The diagram in Figure 11 reflects all predictions from section 4.3. Assume, for example, that a column is operating on conditions corresponding to the 0 in Figure

Figure 12. Illustration of the qualitative correspondence between a two-parameter (xFLI, D′) bifurcation diagram and a one-parameter (D′) bifurcation diagram when xFLI is constant (legend in Table 1).

ing in the cases studied here, they are omitted in the one-parameter (D′) bifurcation diagrams in Figure 11 for clarity. In general, additional information would be needed to construct one-parameter bifurcation diagrams from a two-parameter bifurcation diagram in the presence of MSS. Without it, it could not be decided on which of the branches the other bifurcations occur. In the two-parameter (xFLI, D′) bifurcation diagram in Figure 11, there are bifurcations of codimension one and two. Further, a region where the system will exhibit underdamped oscillations near the equilibria on the high branch is illustrated. Because it is not a bifurcation when two real eigenvalues meet and turn into a pair of complex conjugate eigenvalues, the oscillation region is not bounded by a bifurcation line. As discussed in the previous sections, for the case of 001 class mixtures, the

11. Then, changes in xF such that xFLI decreases would result in the subsequent disappearance of the homoclinic bifurcations, the Hopf bifurcations, and then the oscillation region on the high branch. The lines of the left fold, Hopf, and homoclinic bifurcations end in a Bogdanov-Takens bifurcation point BT for increasing values of xFLI (see Figure 11e). There the left Hopf bifurcation point and the left homoclinic bifurcation point together reach the left turning point and disappear. This corresponds to the transient from case a3 to case a4 in Figure 8 when the upper turning point enters region Ω. Also, the left boundary of the oscillation region touches T in BT. In appendix B, it will be discussed why there must be a Bogdanov-Takens bifurcation between cases d and f. The dynamic behavior of the column in the neighborhood of BT is discussed in appendix B. In some cases, a region where the dominant eigenvalue is positive real will appear on the high branch between the left turning point and the homoclinic point, thus disconnecting the region of oscillation. The left region of oscillation subsequently shrinks and disappears, while the right oscillation region becomes gradually narrower. The value of xFLI for which the oscillation region starts being disconnected depends on the actual mixture. It can be lower (as illustrated in Figure 11) or higher (not illustrated) than the value of xFLI for which


Figure 13. Two-parameter (xFLH, D′) bifurcation diagram with bifurcations of codimension one and two (legend in Table 1). Further, the regions of oscillation are illustrated (gray shaded).

the Bogdanov-Takens bifurcation occurs. In other cases, no disconnection of the region of oscillation is observed. 5.2. Two-Parameter (xFLH, D′) Bifurcation Diagram. The influence of xFLH is twofold, as was discussed in section 4.2. As long as xFLH e xAz LH, the main influence on the bifurcation behavior is to shorten the last part of the distillate product path which corresponds to profiles of type II. Hence, the characteristic features (e.g., region of oscillation or occurrence of Hopf bifurcation points) are nearly unaffected except being shifted correspondingly to higher values of D′. This is illustrated in cases a-c in the two-parameter (xFLH, D′) bifurcation diagram shown in Figure 13. Note that D > F (i.e., D′ > 1) is physically irrelevant. However, for mathematical models (e.g., CMO), D′ > 1 is feasible if B′ < 0 is allowed. Therefore, some of the region D′ > 1 is included in the illustrations because it eases the presentation. Assuming infinite reflux L′ and xFLI high enough such that bifurcations are present, increasing xFLH beyond xAz LH can cause the feed composition to move into region Ω. When xF crosses the boundary δΩ from the outside, the Hopf bifurcation point corresponding to the intersection of the distillate product path and δΩ occurs for D′ ) F′. Increasing xFLH further causes the Hopf bifurcation point to move toward values of D′ > 1, hence moving out of the region of physically meaningful operation. If xFLH is increased so far that the feed composition approaches the boundary δΩ from the inside, the maximum achievable distillate purity in L will increase, finally moving steady states near the left turning point out of region Ω again. Hence, at some value of xFLH, a Bogdanov-Takens bifurcation must occur at the left fold bifurcation. For values of xFLH > (xFLH)BT, there is now a new Hopf bifurcation as well as a new homoclinic bifurcation on the high branch (Figure 13, case d). The steady states between the left fold bifurcation and the Hopf bifurcation are stable.

Further increasing xFLH causes the two homoclinic bifurcation points to meet in a degenerate homoclinic bifurcation and disappear (Figure 13, case e) and the Hopf bifurcation point to move to higher values of D′ and to disappear finally. This happens when xFLH is increased so far that xF is outside region Ω again. In the series of events described above, steps c and d could be exchanged. In that case, the Bogdanov-Takens bifurcation on the left fold bifurcation would occur before xF enters region Ω. The corresponding sequence of bifurcations is not illustrated here. 5.3. Two-Parameter (L′, D′) Bifurcation Diagram. It was discussed in section 4.3 that, for the 001 class mixtures studied here, multiplicities occur for all ternary feed compositions if L′ is sufficiently large and the column has a sufficiently large number of trays. As L′ decreases, the MSS region becomes smaller and finally disappears because of insufficient separation. This is illustrated in Figure 14 for a given feed composition xFLI and xFLH. An appropriate norm ||x|| of the state vector suitable for illustration is plotted over D′ and L′ for all feasible steady states. For small values of L′, there is one and only one steady state for every value of D′. For L′ > L′CP, multiple steady states exist for an interval of D′. The curve T of fold bifurcation points (codimension one) is projected onto T ˆ in the (D′, L′) parameter plane. Both branches of T meet in a cusp bifurcation point CP (codimension two). The curve T is smooth everywhere. The geometric singularity of T ˆ in the cusp bifurcation point is caused by the projection (cf. appendix B). For a given feed composition xFLI and xFLH, the influence of L′ on the bifurcation behavior is illustrated in the two-parameter (L′, D′) bifurcation diagram in Figure 15. There xFLH was chosen such that the feed composition is not located within region Ω (cf. Figure 8b). In addition, xFLI was set at a value sufficiently high such that at high reflux all of the bifurcations discussed in section 4.3 occur as D′ and L′ are varied. In the same

3954


Figure 16. Classification of the dynamic behavior in parameter space L′ - xFLI (cf. Table 3). Table 2. Types of Bifurcation Behavior Corresponding to the Regions Shown in Figure 16

Figure 14. Two-parameter (L′, D′) bifurcation diagram illustrating a line of fold bifurcations T and a cusp bifurcation point CP.

region

description

Figure

I IIa IIb IIIa IIIb IV

unique stable SS MSS MSS, damped oscillations MSS, two Hopf bifurcations MSS, two Hopf and two homoclinic bifurcations MSS, one Hopf and one homoclinic bifurcation

11a 11b 11c 11d 11g

Table 3. Correspondences between Boundaries Separating the Regions in Figure 16 and Bifurcations of Codimension Two boundary between

bifurcation

I and IIa IIa and IIb IIb and IIIa IIIa and IIIb IIIb and IV

cusp none (oscillation region) degenerate Hopf degenerate homoclinic Bogdanov-Takens

other cases, the disconnection of the region of oscillation was not observed. 6. Classes of Dynamic Behavior

Figure 15. Two-parameter (L, D′) bifurcation diagram with bifurcations of codimension one and two (legend in Table 1). Further, the regions of oscillation are illustrated (gray shaded).

way as in the two-parameter (xFLI, D′) bifurcation diagram in Figure 11, all interesting phenomena except for one of the fold bifurcations occur on the high branch of the corresponding one-parameter (D′) bifurcation diagram. Appendix B gives details on why there must be a Bogdanov-Takens bifurcation between cases d and f where lines H and P meet and end on T in BT. In some cases, a region where the dominant eigenvalue is positive real will appear on the high branch between the left turning point and the homoclinic point thus disconnecting the region of oscillation. The left region of oscillation subsequently shrinks and disappears, while the right oscillation region becomes gradually narrower. The value of L′ for which the oscillation region starts being disconnected depends on the actual case of mixture. It can be lower (as illustrated in Figure 15) or higher (not illustrated) than the value of L′ for which the Bogdanov-Takens bifurcation occurs. For

In the preceding section, it was illustrated that the qualitative behavior of the one-parameter (D′) bifurcations is mainly influenced by xFLI and L′, whereas the main effect of xFLH is to shift both the MSS region and the other bifurcations without changing their main characteristics. Increasing xFLH above certain values causes some bifurcations to disappear from the feasible interval of D′ as they would occur for D > F. This allows one to focus on L′ and xFLI as an example to characterize the types of dynamic behavior. In this section, it will be studied how the (L′, xFLI) parameter space is divided into regions corresponding to qualitatively different steady state and dynamic behavior. For the example studied in this paper, the regions are illustrated in Figure 16. The properties of the regions are listed in Table 2. The diagram shown in Figure 16 is a qualitative illustration, constructed using the predictions presented in section 4. Because one of the parameters is L′, a quantitative prediction of this diagram is not possible (cf. sections 3 and 4.3). Figure 16 can also be interpreted more formalized. Spanning a three-dimensional parameter space with xFLI, L′, and D′, the boundaries separating the regions in Figure 16 are the projections of the lines of bifurcation points of codimension two onto the (xFLI, L′) space. The correspondences are given in Table 3. Similar to Figure 16, the lines of bifurcation points of codimension two in


the three-dimensional parameter spaces (xFLI, xFLH, D′) and (xFLH, L′, D′) can be projected onto the (xFLI, xFLH) and (xFLH, L′) spaces, respectively. 7. Summary From an application point of view, a qualitative bifurcation analysis for a given mixture starts with determining the corresponding residue curve map and class27 of the mixture to be studied. Given a feed composition, the corresponding product paths and bifurcation diagram for the infinite reflux case are determined with the ∞/∞ analysis (section 2.1). Classifying the different types of profiles along the product paths (corresponding to different sections of the branches in the bifurcation diagram), the stability of each type can be determined with the qualitative stability analysis (section 2.2). Note that the qualitative stability analysis is applied to classes of mixtures. Hence, the stability of profiles does not need to be determined again if it has been examined earlier for other mixtures belonging to the same class. Bifurcation points in the bifurcation diagram can be determined geometrically (e.g., fold bifurcations) or by studying qualitative changes of schematic vector fields near equilibria exchanging stability (section 2.3). In a next step, the effect of dropping the assumption of infinite reflux can be evaluated. After construction of the two-parameter (L, D′) bifurcation diagram as shown in Figure 15, one or more numerical simulations are needed to locate the operating reflux in the diagram. As discussed before, the reason for this is that the influence of the reflux flow rate can only be studied qualitatively. The theory does not predict, however, at which reflux flow rate, for example, multiple steady states of homoclinic bifurcations begin or cease to exist. Finally, the influence of changes in the feed composition can be studied by constructing the two-parameter (xFLH, D′) and (xFLI, D′) bifurcation diagrams as shown in Figures 11 and 13.

is chosen as a scalar index for the state variables of the column: n

Tavg(t) )

Tj(t) Mj ∑ j)1 n

Mj ∑ j)1

8. Case Studies

This quantity was chosen to provide some engineering insight. Obviously, any other state variable could have been used. An advantage of choosing Tavg is that it is monotonically increasing along the continuation path for all examples studied in this work. Therefore, the resulting bifurcation diagrams have the typical S-curve shape which facilitates illustrations. The TCHC xavg is used for illustrations of the trajectories in phase diagrams. Because a rather large number of one- and twoparameter bifurcation diagrams will be shown, a legend explaining the symbols used in the illustrations is provided centrally in Table 1. 8.1. MMBT. Until today, the mixture MMBT is the only homogeneous azeotropic mixture for which the existence of multiple steady states has been studied experimentally.28,29 Further, it was the first mixture for which simulations were published showing the existence of limit cycles in open-loop models of homogeneous azeotropic distillation columns.6 The binary minimumboiling azeotrope between methanol and toluene is located at 88 mol % methanol (xAz LH ) 7.33). To study the influence of xFLI, a bifurcation diagram (shown in Figure 17-I) is computed for xFLI ) 10, xFLH ) 6, and L′ ) 48 (corresponding to region IIa in Figure 16). Increasing xFLI results in the successive appearance of an oscillation region on the high branch (Figure 17II), two Hopf bifurcation points (Figure 17-III), two homoclinic bifurcation points between the Hopf bifurcation points (not illustrated), and a Bogdanov-Takens bifurcation resulting in the disappearance of the left Hopf and homoclinic bifurcation points (Figure 17-IV). For this sequence of bifurcation diagrams, the correspondences are summarized in Table 4. In Figure 18,

In the following, two mixtures of the class 001 are studied with steady-state and dynamic simulations to confirm and illustrate the bifurcations predicted in the previous sections. The mixtures used are methanolmethyl butyrate-toluene (MMBT) and methanol-butyraldehyde-benzene (MBB). While most simulations are performed using the CMO model, some simulations performed with RADFRAC, the distillation column model included in Aspen PLUS V.10 and Aspen Dynamics V.10, are also included. All simulations with the CMO model are performed open loop, where open loop refers to perfect level control (i.e., constant molar holdup) in the condenser and the reboiler. Hence, there are no controllers closing any loop around the column and hence no controller dynamics that could interact with the column dynamics. Simulation details are given in appendix C. In addition to the results presented here, the mixture acetone-benzene-heptane (ABH) was studied in detail in an earlier publication.25 All findings presented in that paper agree with the theory presented in this work. For the bifurcation diagrams, the weighted average temperature of all trays, the reboiler, and the condenser

Figure 17. MMBT: Series of bifurcation diagrams (showing Tavg over D′) for different values of xFLI. xFLH ) 6 and L′ ) 48 (legend in Table 1).

3956


Figure 18. MMBT: Distillate product paths corresponding to the bifurcation diagrams in Figure 17. (I) xFLI ) 10; (II) xFLI ) 20; (III) xFLI ) 29.9; (IV) xFLI ) 80. Figure 20. MMBT: Series of bifurcation diagrams for different values of L′. xFLI ) 80 and xFLH ) 6 (legend in Table 1).

Figure 19. MMBT: Zoom-in of the four cases shown in Figure 18 with steady-state column profiles for fold and Hopf bifurcation points. Compare Figure 18 for legend and units of the axes.

Figure 21. MMBT: Response of average liquid composition in the column starting at an unstable steady state on the high branch for D′ ) 0.97.6 xFLI ) 173.4, xFLH ) 6.8, and L′ ) 12.

Table 4. Correspondences for the Series of Bifurcation Diagrams Shown in Figure 17

The effect of increasing L′ on the bifurcation behavior is studied in a similar series of bifurcation diagrams illustrated in Figure 20. The first (shown in Figure 20I) is computed for the operating conditions xFLI ) 80, xFLH ) 6, and L′ ) 1 (corresponding to region IIa in Figure 16). Increasing L′ results in the successive appearance of an oscillation region on the high branch (L′ ) 8; Figure 20-II), two Hopf bifurcation points (L′ ) 16; Figure 20-III), two homoclinic bifurcation points between the Hopf bifurcation points (L′ ) 20; Figure 20-IV), and a Bogdanov-Takens bifurcation resulting in the disappearance of the left Hopf and homoclinic bifurcation points (L′ ) 48; Figure 20-V, which is the same case as that in Figure 17-IV). The correspondences are summarized in Table 5. In Figures 21-24, dynamic simulations are presented to examine in detail the behavior of a distillation column going through a limit cycle. The operating conditions xFLI ) 173.4, xFLH ) 6.8, D′ ) 0.97, and L′ ) 12

bifurcation diagram in Figure 17 case in Figure 11 region in Figure 16

I a IIa

II b IIb

III c IIIa

IV g IV

the corresponding distillate product paths are presented together with residue curves. The locations of the Hopf bifurcation points along the distillate product path agree with the predictions that were illustrated in Figure 8a. For each of the four cases, the distillate product path is illustrated separately in Figure 19. The column profiles for the equilibria corresponding to the left limit point and the Hopf bifurcation points (if applicable) are also shown. Especially for the case xFLI ) 80, the profiles will be nonmonotonic in all three components (type IIu) on both sides of the left limit point. This means that the predictions agree with Figure 17-IV, where the steady states on both sides of the left fold bifurcation point were found to be unstable.

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3957 Table 5. Correspondences for the Series of Bifurcation Diagrams Shown in Figure 20 bifurcation diagram in Figure 20 case in Figure 15 region in Figure 16

I a IIa

II b IIb

III c IIIa

IV d IIIb

V g IV

correspond to the first published case of sustained oscillations in open-loop homogeneous azeotropic distillation.6 Figure 21 shows the time domain responses of the TCHC in the column to a small perturbation in D′ starting from the unstable steady state on the high branch. The period of oscillation is ≈110 h. In the upper left subfigure of Figure 22, the boundary of region Ω is illustrated in the composition space together with the limit cycle of the distillate composition. Note that region Ω was determined by assuming infinite reflux for simplicity. In regions 4-6 (cf. Figure 5), the limit cycle approximately coincides with the boundary δΩ. The difference (cf. Figure 6) is due to the fact discussed in section 4.3 that region Ω shrinks as the reflux flow rate is reduced. In the other subfigures in Figure 22, column profiles are plotted at selected times. As the distillate composition moves from region 1 to 6 and finally back into region 1, the shapes of the profiles are in excellent

agreement with the theoretical predictions for Hopf bifurcations and limit cycles given in section 2.3, if the influence of the finite reflux flow rate on the shape of column profiles is considered (section 4.3). Figure 23 illustrates the time domain response of the liquid composition on tray 12 in the composition space. In Figure 24, the movement of the temperature front inside the column during the oscillation is given at intervals of ∆t ) 10 h. Note that the velocity and shape of the front moving up and down is different. Simulations confirming the findings with the more complex column model RADFRAC were also published by Lee et al.6 8.2. MBB. The homogeneous mixture MBB was studied by Bekiaris and Morari.30 The minimum-boiling binary azeotrope between methanol and benzene is at 61.31 mol % methanol (xAz LH ) 1.58). Note that there is uncertainty whether methanol and butyraldehyde form an azeotrope or not. Our decision to neglect that ambiguous azeotrope is discussed in appendix C. All results for the mixture MMBT (e.g., shown in Figures 17 and 20) were reproduced for the mixture MBB.8 In the following, the influence of xFLH on the bifurcation behavior of the column is studied. A starting point is the bifurcation diagram illustrated in Figure 25-I.

Figure 22. MMBT: Column profiles during the limit cycle. The distillate composition moves from region 1 to 6 and back into region 1 (see also Figure 21).

3958


Figure 23. MMBT: Phase plot of liquid composition on tray 12 starting at an unstable steady state on the high branch6 (see also Figure 21).

Figure 25. MMBT: Series of bifurcation diagrams for different values of xFLH. xFLI F ) 605, and L′ ) 5 (legend in Table 1).

Figure 24. MMBT: Temperature profiles (∆t ) 10 h) during the sustained oscillation6 (see also Figure 21).

This case corresponds to Figure 13a. Increasing xFLH lets the feed composition, and hence the end of the distillate product path, move closer and finally into region Ω. Thus, the right Hopf bifurcation point moves to higher values of D′ and finally disappears as the feed composition crossed δΩ. This is shown in parts II and III in Figure 25 corresponding to parts b and c of Figure 13. At the same time, the maximum achievable purity in methanol increases as a result of the feed being more pure in methanol. Hence, the left limit point approaches δΩ from the inside, at first causing an oscillation region to appear. Crossing δΩ corresponds to a BogdanovTakens bifurcation. This is illustrated in parts IV and V of Figure 25, corresponding to parts c and d of Figure 13. Further increasing xFLH causes the feed composition to approach δΩ from the inside and finally crossing it, resulting in the Hopf bifurcation point (that was born in the Bogdanov-Takens bifurcation) to move along the high branch toward higher values of D′ and finally to disappear as it moves out of the interval of feasible D′. This is illustrated in parts VI and VII of Figure 25, corresponding to parts e and f of Figure 13. Further, I

refers to case a in Figure 13, II to case b, III and IV to case c, and V to case d. While the starting point of the series of bifurcation diagrams illustrated in Figure 25 corresponded to region I in Figure 16, increasing xFLH further led to situations not corresponding to any of the regions in Figure 16. A classification of those cases into regions is possible but will inevitably lead to complications because in the end the three-dimensional parameter space would have to be subdivided into three-dimensional subregions. This approach is not pursued here because it is hard to illustrate without giving additional insight. The bifurcation diagrams illustrated in Figure 25 indicated the existence of a situation newly observed in open-loop models of distillation columns: For certain values of xFLH, operating conditions exist where the column has only one steady state. However, different to all published cases known to the authors, that single steady state is unstable and surrounded by a stable limit cycle. For the feed composition used to compute the bifurcation diagram shown in Figure 25-VI, nonlinear dynamic simulations were performed at D′ ) 0.999 99 for various initial conditions. Figure 26 shows selected phase plane trajectories of the average liquid composition in the column during each simulation. A pronounced stable limit cycle surrounds the unstable steady state. The simulation packages Aspen PLUS V.10 and Aspen Dynamics V.10 by Aspen Technology, Inc., were used to reproduce some steady-state bifurcation dia-


between the two homoclinic bifurcation points. For the first hour of the dynamic simulations, D′ is decreased or increased by 10-5D′. As expected, all transients end on the corresponding steady state on the low branch, which is the only attractor for values of D′ within the interval between the two homoclinic bifurcation points. 9. Conclusion

Figure 26. MBB: Phase diagram of trajectories to different initial conditions in an operating region where the only steady state is unstable and surrounded by a stable limit cycle. D′ ) 0.999 99, xFLI ) 605, xFLH ) 1.7, and L′ ) 5.

A novel approach for the qualitative analysis of the dynamic behavior of homogeneous azeotropic distillation columns separating mixtures of the class 001 was presented. Based on the methodology (published in the preceding paper1) to determine the local stability of column profiles, the global behavior of the column was examined by construction of bifurcation diagrams. Phenomena like Hopf bifurcations, limit cycles, and homoclinic bifurcations were predicted and explained. The influence of the operating conditions (feed composition and reflux flow rate) was analyzed. An extension of the methodology for a qualitative analysis of mixtures whose residue curve maps have more complex topologies is in preparation.26 It is important to point out that the insight is gained on physical grounds without numerical simulations of the distillation column. Instead, all results can be derived by interpretation of the residue curve map of the corresponding mixture. Note that the residue curve map illustrates some physical properties of a mixture and therefore does not depend on the thermodynamic model used because it could also be measured. Acknowledgment The authors thank Jan Ulrich for his support in the preparation of the final version of this paper and Prof. Moonyong Lee for his contributions to the case studies presented in section 8. Further, the authors thank Thomas E. Gu¨ttinger, who programmed the distillation column model used for all CMO simulations. Appendix A: Geometrical Detection of Fold Bifurcations

Figure 27. MBB: Steady-state bifurcation diagram and dynamic simulations showing the average temperature in the column for L′ ) 2 (Aspen PLUS, Wilson).

grams and the dynamic simulations, respectively. RADFRAC, the distillation column model included in these packages, is a model that includes mass and energy balances, volumetric holdups, rigorous pressure drop calculations, level control for reboiler and condenser, pressure control, etc. For many purposes, this model can be called “rigorous” in comparison to the “simple” CMO model. Compared to the CMO model, most of the design parameters such as the number of trays, feed composition, etc., were left unchanged. Other parameters that do not exist in the simple model like volumetric holdups were chosen in a way to give a meaningful approximation of the CMO case. Simulations were performed for the feed composition xFLI ) 600 and xFLH ) 1.5, i.e., xF ) (0.5995, 0.001, 0.3995). For the case L′ ) 2 shown in Figure 27, the stability of the steady states was checked with dynamic simulations using Aspen Dynamics. As illustrated, the column is initialized on steady states on the middle and high branches for D′ ) 0.967, which is in the interval

In one-parameter bifurcation diagrams (e.g., Figure 3), fold bifurcation points can easily be detected as turning points. They can also be detected directly in the product paths in the composition space by reformulation of the geometric condition for multiple steady states.3 In Figure 28, a part of the distillate and bottoms product paths are shown. When D is chosen as the bifurcation parameter, the distillate product path starts (D ) 0) at

Figure 28. Geometric condition for the detection of fold bifurcations in the product paths.

3960


the unstable node of the residue curve map while the bottoms product path ends (D ) F) at the stable node. Consider an equilibrium with a corresponding profile and product compositions xD(2) and xB(2). A condition for a generic fold bifurcation to occur at that equilibrium is that the tangent vectors to the product paths in xD(2) and xB(2) are in an angle of 180°. For nondegeneracy of the fold bifurcation, the angle between the two tangent vectors must not have a maximum or minimum in the bifurcation point. Appendix B: On Bifurcations of Codimension Two in Homogeneous Azeotropic Distillation In this appendix, some aspects of the bifurcation points of codimension two reported in sections 4 and 5 are discussed in more detail. The cusp bifurcation served as an example in section 5.3 and is studied in detail in Figure 14. It is seldom recognized as a bifurcation under its own name, although it is present in most systems where the number of steady states depends on the value of continuous parameters. In analogy to the cusp bifurcation that can be interpreted as the birth of twofold bifurcations, the degenerate Hopf bifurcation is the birth of two Hopf bifurcation points. There are different ways for a Hopf bifurcation point to degenerate, but no clear terminology exists in the literature. The physical understanding of the process allows one to determine in which way that the Hopf bifurcation is degenerate in the cases studied here: the transversality condition is violated. The different role of D′ and xFLI as bifurcation parameters is obvious. While D′ is used to parametrize the product path and thus localize the Hopf bifurcation points, xFLI influences whether there are Hopf bifurcations at all (cf. Figure 8a). The transient from zero to two Hopf bifurcations occurs in case a2, where the high branch is stable except for the point in which the distillate product path touches δΩ. The corresponding equilibrium is semistable. The only way in which this can happen is illustrated in Figures 29a29c. Starting with an equilibrium on the middle branch, the eigenvalue (1) with the largest real part is in the right half-plane. At the fold bifurcation, it crosses through the origin into the left half-plane. Going further, it meets with the eigenvalue (2) with the second biggest real part and splits into a pair of complex conjugate eigenvalues. The pair of eigenvalues approaches the imaginary axis from the left. In Figure 29a, corresponding to case a1 in Figure 8, the eigenvalues do not cross the imaginary axis. All steady states on the high branch are stable. Figure 29c (cf. case a3 in Figure 8) corresponds to two Hopf bifurcations on the high branch. Between these two cases, the degenerate Hopf bifurcation occurs (Figure 29b; cf. case a2 in Figure 8). The argumentation for the degenerate homoclinic bifurcation is similar. However, because the homoclinic bifurcation is not defined by specific locations of eigenvalues, it is difficult to illustrate and details are omitted here. Increasing xFLI further led to case a4 in Figure 8a. The qualitative stability analysis of that case found that the equilibria on the high branch between the left fold bifurcation and the only Hopf bifurcation point are unstable. The eigenvalues or that case are illustrated in Figure 29e. At the fold bifurcation, the second eigenvalue (2) crosses through the origin into the right

Figure 29. Illustration of the movement of the two most dominant eigenvalues of the steady states on the high branch as the distillate flow rate is varied: (b) degenerate Hopf bifurcation; (d) BogdanovTakens bifurcation.

half-plane. Going further, it meets with eigenvalue (1), splits into a pair of complex conjugate eigenvalues, and crosses the imaginary axis, corresponding to the Hopf bifurcation point. The situation illustrated in Figure 29e emerged from case c under variation of xFLI. The point in the complex plane in which the two real eigenvalues meet after the fold bifurcation and split into a pair of complex conjugate eigenvalues moved from the left half-plane (Figure 29c) into the right half-plane (Figure 29e). Assuming smoothness of a system, the eigenvalues of the local linearizations move continuously under variation of the system parameters. Hence, there must be a xFLI ) (xFLI)BT in which for D′ ) (D′)BT there are two eigenvalues in the origin (Figure 29d). Thus, there is a Bogdanov-Takens bifurcation. This type of bifurcation can be interpreted as a fold and a Hopf bifurcation occurring simultaneously. In addition, there is also a homoclinic bifurcation contained in the BogdanovTakens bifurcation. The Bogdanov-Takens bifurcation point BT is illustrated in more detail in Figure 30. The codimensionone bifurcation lines meeting in BT separate the (xFLI, D′) parameter space locally in four regions. Regions 2 and 4 are further separated by the boundary (sometimes called a nonbifurcation line) of the oscillation region. The illustrations of the phase portraits corresponding to these regions show the main characteristics of the local neighborhood of the equilibria for xFLI and D′ near the Bogdanov-Takens point. Note that region 1 in Figure 30 corresponds to operating conditions outside the region of multiplicities; i.e., there is only the one steady state corresponding to the lower branch. Locally, there is no equilibrium because there is no middle or high branch. Appendix C: Simulation Details All simulations were performed using the CMO model as described in the preceding paper.1 For the dynamic simulations, the CMO model is integrated with


Figure 30. Illustration of the two-parameter (xF, D′) bifurcation from Figure 11 close to the Bogdanov-Takens bifurcation point (legend in Table 1). Table 6. Column Design Used for the Simulations no. of trays (including condenser and reboiler) feed tray (counting from condenser) tray liquid holdup [kmol] condenser liquid holdup [kmol] reboiler liquid holdup [kmol] column pressure [atm] feed flow rate [kmol/h]

46 41 3 3 3 1.0 100

DDASAC,31 an index-1 DAE solver. AUTO32 is used for a steady-state bifurcation analysis of the CMO model. To investigate the local dynamic behavior and stability, the eigenvalues are calculated with LAPACK33 by linearizing the nonlinear CMO model at every steady state. The simulations are based on the column design given in Table 6. The vapor phase is modeled as an ideal gas. Vapor pressures and liquid activity coefficients are calculated with the Antoine equation and the Wilson model, respectively. The Wilson parameters for the mixture MMBT are taken from Gu¨ttinger et al.28 The Antoine coefficients are from Thermopack (thermodynamic subroutines provided by Prof. M. F. Doherty and J. Knapp, University of Massachusetts, Amherst, MA). For the mixture MBB, the Antoine coefficients are from Thermopack. Bekiaris and Morari30 discussed the ambiguity of the azeotrope between methanol and butyraldehyde. In a massive collection of azeotropic data,34 it is reported that methanol and butyraldehyde do not form any azeotrope. However, the reference cited35 must be wrong or mistaken because butyraldehyde is not mentioned in it at all. Another collection of information on azeotropes36 does not contain any information on the mixture of methanol and butyraldehyde. In conclusion, the author is not aware of reliable information on whether that azeotrope exists or not. Another approach, yet less reliable, is to use group contribution methods. Using Aspen PropertiesPLUS, the UNIFAC model predicts a maximum-boiling azeotrope between methanol and butyraldehyde while UNIFAC with the Dortmund corrections (UNIFAC-DMD) does not. In the first case, the ternary mixture would be of the class 301-S and in the latter of the class 001. Here, we will put more trust in the predictions from UNIFAC-DMD. For the Wilson model, the coefficients for the binaries methanol-benzene and butyraldehyde-benzene were taken from ASPEN PropertiesPLUS. For the binary methanol-butyraldehyde, the coefficients were regressed with the AspenPLUS properties estimation

option from data generated with UNIFAC with Dortmund corrections.37 Notation n ) number of trays nc ) number of components L ) light boiling component I ) intermediate boiling component H ) heavy boiling component Mk ) holdup on tray k kF ) feed tray F ) feed flow rate xF ) feed composition xFLI ) xFL/xFI xFLH ) xFL/xFH L, V ) reflux and boilup flow rate D, xD ) distillate flow rate and composition B, xB ) bottoms flow rate and composition D′ ) D/F L′ ) L/F xk ) liquid composition on tray k yk ) vapor composition in equilibrium with xk xki ) fraction of component i in xk Ω ) set of all xD of profiles of type IIu see also Table 1 Abbreviations CMO ) constant molar overflow TCH ) total column holdup TCH in i ) total column holdup in component i TCHC ) total column holdup composition ABH ) acetone-benzene-heptane MBB ) methanol-butyraldehyde-benzene

Literature Cited (1) Dorn, C.; Morari, M. Qualitative Analysis for Homogeneous Azeotropic Distillation. 1. Local Stability. Ind. Eng. Chem. Res. 2002, 41, 3930-3942. (2) Petlyuk, F.; Avet’yan, V. Investigation of the Rectification of Three-Component Mixtures with Infinite Reflux. Theor. Found. Chem. Eng. 1971, 5 (4), 499-507. (3) Bekiaris, N.; Meski, G.; Radu, C.; Morari, M. Multiple Steady States in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1993, 32 (9), 2023-2038. (4) Ulrich, J.; Morari, M. Influence of the Impurities on the Control of Heterogeneous Azeotropic Distillation Columns. Ind. Eng. Chem. Res. 2002, 41 (2), 230-250.

3962


(5) Bekiaris, N.; Gu¨ttinger, T.; Morari, M. Multiple Steady States in DistillationsVLLE Inaccuracies, Multiplicities and ∞/∞ Predictions. AIChE J. 2000, 46 (5), 955-979. (6) Lee, M.; Dorn, C.; Meski, G.; Morari, M. Limit Cycles in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1999, 38 (5), 2021-2027. (7) Kuznetsov, Y. A. Elements of applied bifurcation theory, 2nd ed.; Springer: New York, 1998. (8) Dorn, C. Nonlinear Dynamics in Homogeneous Azeotropic Distillation. Ph.D. Thesis, ETH No. 13947, ETH, Zu¨rich, Switzerland, 2000; available via http://www.control.ethz.ch/. (9) Takens, F. Forced Oscillations and Bifurcations. Applications of Global Analysis. Comm. Math. Inst. Rijksuniversitat Utrecht 1974, 3, 1-59. (10) Bogdanov, R. Versal Deformations of a Singular Point on the Plane in the Case of Zero Eigenvalues. Sel. Math. Sov. 1981, 1 (4), 389-421. (11) Keyfitz, B. Shocks Near the Sonic Line. A Comparison Between Steady and Unsteady Model for Change Type. IMA Appl. Math. 1990, 27, 88-106. (12) Ondarç uhu, T.; Mindlin, G.; Mancini, H.; Pe´rez-Garcıá, C. The Chaotic Evolution of Patterns in Beńard-Marangoni Convection with Square Symmetry. J. Phys.: Condens. Matter 1994, 6, A427-A432. (13) Ehrenstein, U.; Koch, W. Homoclinic Bifurcation in Blasius Boundary-Layer Flow. Phys. Fluids 1995, 7 (6), 1282-1291. (14) Labate, A.; Ciofini, M.; Meucci, R.; Boccaletti, S.; Arecchi, F. Pattern Dynamics in a Large Fresnel Number Laser Close to Treshold. Phys. Rev. A 1996, 56 (3), 2237-2241. (15) Degtiarev, E.; Wataghin, V. Takens-Bogdanov Bifurcation in a Two-Component Nonlinear Optical System with Diffractive Feedback. J. Mod. Opt. 1998, 45 (9), 1927-1942. (16) Holmes, P. Bifurcation in Coupled Oscillators with Applications to Flutter and Divergence. ASME Nonlinear Syst. Anal. Synth. 1980, 2, 389-416. (17) Yagasaki, K. Dynamics of a Simple Model for a WindLoaded Nonlinear Structure: Bifurcations of Codimension One and Two. J. Appl. Mech. 1998, 65 (6), 505-512. (18) Merryfield, W.; Toomre, J.; Gough, D. Nonlinear Behavior of Solar Gravity Modes driven by 3He in the Cores. I. Bifurcation Analysis. Astrophys. J. 1990, 353, 678-697. (19) Rodriguez-Luis, A.; Freire, E.; Ponce, E. On a Codimension 3 Bifurcation Arising in an Autonomous Electronic Circuit. In Bifurcation and Chaos; Seydel, R., Ed.; Birkhaüser Verlag: Basel, Switzerland, 1991; Vol. 97, pp 301-306. (20) Freire, E.; Gamero, E.; Rodriguez-Luis, A. Study of a Degenerate Bogdanov-Takens Bifurcation in a Family of Mechanical Oscillators. Mech. Res. Commun. 1998, 25 (3), 287-297. (21) Burchard, A. Substrate Degradation by a Mutualistic Association of Two Species in the Chemostat. J. Math. Biol. 1994, 32, 465-489. (22) Gragnani, A.; Rinaldi, S.; Feichtinger, G. Cyclic Dynamics

in Romantic Relationships. Int. J. Bifur. Chaos 1997, 7 (11), 26112619. (23) Schu¨tz, P.; Bode, M.; Purwins, H.-G. Bifurcations of Front Dynamics in a Reaction-Diffusion System with Spatial Inhomogeneities. Phys. D 1995, 82, 382-397. (24) Li, R.-S.; Liu, Q.-H. Takens-Bogdanov Bifurcation in TwoVariable Chemical Reaction Systems. J. Chem. Phys. 1992, 97 (5), 3871-3872. (25) Dorn, C.; Lee, M.; Morari, M. Stability and Transient Behavior of Homogeneous Azeotropic Distillation. Comput. Chem. Eng. 1999, 21 (Suppl.), S185-S188. (26) Dorn, C.; Morari, M. Qualitative Analysis for Homogeneous Azeotropic Distillation. 3. Extension to Complex Mixtures. In preparation, 2002. (27) Matsuyama, H.; Nishimura, H. Topological and Thermodynamic Classification of Ternary Vapor-Liquid Equilibria. J. Chem. Eng. Jpn. 1977, 10 (3), 181-187. (28) Gu¨ttinger, T.; Dorn, C.; Morari, M. Experimental Study of Multiple Steady States in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1997, 36 (3), 794-802. (29) Dorn, C.; Gu¨ttinger, T.; Wells, J.; Morari, M.; Kienle, A.; Klein, E.; Gilles, E. Stabilization of an Azeotropic Distillation Column. Ind. Eng. Chem. Res. 1998, 37 (2), 506-515. (30) Bekiaris, N.; Morari, M. Multiple Steady States in Distillations∞/∞ Predictions, Extensions and Implications for Design, Synthesis and Simulation. Ind. Eng. Chem. Res. 1996, 35 (11), 4264-4280. (31) Stewart, W.; Caracotsios, M.; Sorensen, J. DDASAC Software Package Documentation, 1995. (32) Doedel, E.; Champneys, A.; Fairgrieve, T.; Kuznetsov, Y.; Standstede, B.; Wang, X. AUTO97: Continuation and Bifurcation Software for Ordinary Differential Equations; Concordia University: Montreal, Canada, 1997 software documentation. (33) Anderson, E.; Bai, Z.; Bischof, C.; Demmel, J.; Dongarra, J. Lapack Users’ Guide: Release 2.0; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1995. (34) Horsley, L. Azeotropic Data; Advances in Chemistry Series 116; American Chemical Society: Washington, DC, 1973; Vol. III. (35) Carleton, L. Phase Equilibria in DimethylhydrazineWater System. Chem. Eng. Data Ser. 1956, 1, 21-24. (36) Gmehling, J.; Menke, J.; Krafczyk, J.; Fischer, K. Azeotropic Data; VCH: Weinheim, Germany, 1994. (37) Aspen Plus Release 9 Reference Manual: Physical Property Methods and Models; Aspen Technology Inc.: Cambridge, MA, 1995.

Received for review September 4, 2001 Revised manuscript received April 23, 2002 Accepted May 15, 2002 IE010726J

Qualitative Analysis for Homogeneous Azeotropic ... - ACS Publications

Recommend Documents