Multiple-Steady-States Identification in Homogeneous Azeotropic

4386

Ind. Eng. Chem. Res. 2005, 44, 4386-4399

Multiple-Steady-States Identification in Homogeneous Azeotropic Distillation Using a Process Simulator A. Kannan,* Manish R. Joshi,† G. Ramalinga Reddy,‡ and Denish M. Shah Department of Chemical Engineering, Indian Institute of Technology Madras, Chennai 600 036, India

Homogeneous azeotropic distillation processes may exhibit multiple steady states under certain conditions. A methodology has been described for directly capturing all the multiple steady states as well as for generating the bifurcation diagram using the HYSYS process simulator in the steady-state mode. Initial guesses for generating the multiple steady states were based on the ∞/∞ approach of Bekiaris et al. (Ind. Eng. Chem. Res. 1993, 32, 2023) and the behavior of the composition-temperature profiles. The method was tested with different systems, feed locations, and compositions and was extended to nonideal stages using the overall efficiency approach. Conventional inside-out algorithms in the steady-state mode were found to capture the stable steady states only. Using stabilizing control, these algorithms could capture the unstable steady states as well. Their ability to incorporate Murphree efficiencies in the column specifications could, hence, be exploited. Preliminary results on the sustained oscillatory behavior on the higher branch are presented. 1. Introduction Homogeneous azeotropic distillation, an enhanced separation process, is adopted for azeotropic mixtures or those with relative volatilities close to unity. An entrainer is used to accomplish the separation. Homogeneous azeotropic distillation, in contrast to the heterogeneous mode, does not lead to the formation of two liquid phases anywhere in the distillation column or its accessories such as the decanter, making the operation easier. However, it should also be noted that a suitable entrainer, which will not form azeotropes with one or both components or remain homogeneous, might not be easy to identify for the system under consideration. Homogeneous azeotropic distillation processes, the focus of the present work, may exhibit multiple steady states under certain conditions, and these have significant bearing on the column design, operation, and control. Hence, it is important to identify these steady states in an unambiguous manner using commonly available tools. 2. Definition and Occurrence of Multiple Steady States 2.1. Steady-State Simulations. Kienle and Marquardt1 define multiplicity in azeotropic distillation columns as ambiguous temperature and concentration profiles for some fixed system parameters and input variables. Gani and Jorgensen2 defined three different multiplicities that are listed below. a. Output Multiplicity: More than one set of output variables (such as temperature and compositions) is found for the same specified set of input variables. b. Input Multiplicity: For a subset of output variables, several possible values of a subset of input variables exist. * To whom correspondence should be addressed. Tel.: +91 44 2257 4170. Fax: +91 44 2257 0545. E-mail: kannan@ iitm.ac.in. † Corporate R&D Centre, Bharat Petroleum Corporation Limited, Greater Noida 201 306, India. ‡ Cognizant Technology Solutions, Chennai 600 113, India.

c. Internal State Multiplicity: For a given set of specified variables, several values of a subset of variables describing the internal state of the process are obtained. All the multiple solutions obtained should, however, have identical output values. If not, the multiplicity is defined as output multiplicity. The work reported by Magnussen et al.3 relating to multiplicity in heterogeneous azeotropic distillation spurred intense research in this area. Multiple steady states have been reported in binary distillation columns assuming constant molar overflow and ideal vaporliquid equilibria when some column specifications, such as reflux, were given in mass units rather than molar units (Jacobsen and Skogestad4). Gani and Jorgensen2 observed that even though the nonlinear equations describing the separation can give multiple solutions within a defined region, only one solution might be physically meaningful. Even without taking into account the nonlinearity induced by a rigorous thermodynamic model, they showed that the possibility of multiple solutions increases with an increasing number of stages. The nonlinear term in the model equation was related to the number of stages but not to the number of components. A graphical procedure termed as the ∞/∞ analysis was developed by Bekiaris et al.5 to analyze the multiple steady states of ternary homogeneous systems. Specifications that encouraged multiple-steady-state occurrences were made in terms of molar flow rates and infinitely long columns operating at total reflux. A geometrical, necessary, and sufficient condition for multiplicity was derived. Bekiaris and Morari6 suggested using initial estimates from the ∞/∞ analysis in column simulations; two stable steady states could be located, but the third steady state was not traceable, which was ascribed by the authors to its unstable nature. Guttinger and Morari7 extended the ∞/∞ analysis for specifications in terms of mass/volumetric flow rates. Stable-steady-state multiplicities were verified experimentally by Guttinger et al.8 More recently, Vadapalli and Seader9 used the arc length continuation

10.1021/ie049443s CCC: $30.25 © 2005 American Chemical Society Published on Web 05/14/2005

Ind. Eng. Chem. Res., Vol. 44, No. 12, 2005 4387

algorithm implemented in a FORTRAN subroutine in combination with Aspen Plus10 for computing steadystate branches in a region of multiplicity. 2.2. Dynamic Simulations. Bekiaris and Morari6 outlined the potential control issues that need to be addressed in a column exhibiting multiple steady states. These include difficulty in control, areas of attraction of potential steady states, suitable startup strategy, and control-design interactions. In an earlier study, Dorn et al.11 applied proportional-integral (PI) control to stabilize an experimental distillation column operating in the unstable regime. The control scheme was based on the propagation of temperature fronts in the column based on the boiling points of the three-component azeotropic system methanol-methyl butyrate (entrainer)-toluene. Details were not available in this work on how the dynamic simulator (Diva) handled Murphree efficiencies. Moonyong Lee et al.12 reported limit cycles and traveling waves inside a homogeneous azeotropic distillation column for the first time. Dorn and Morari13 defined a vector of average composition termed as the total column holdup composition (TCHC). They showed that, at infinite reflux, the column profiles coincided with distillation lines due to fast dynamics, while the dynamics of the total column holdup composition were slow until the global material balance was satisfied. The dynamic response of a column to perturbation was analyzed in terms of the accumulation/depletion of components and the accompanying direction of movement of the distillate composition for two cases. If the direction of response of the distillate composition was opposite to that of the perturbation, the resulting state was stable; otherwise, it was unstable. These findings were verified by dynamic simulations. In a subsequent study, Dorn and Morari14 described the construction of bifurcation diagrams with the prediction of several complex bifurcation phenomena. They also have shown through open-loop dynamic simulations that the higher steady-state branch on the bifurcation diagram can show the occurrence of a supercritical Hopf bifurcation, a homoclinic bifurcation, etc. Under these conditions, a part of the higher-steady-state branch usually considered to be stable can become unstable and be surrounded by a stable limit cycle. These authors also studied the influence of the reflux flow rate and feed composition on the bifurcation diagrams. For quantifying multiple-steady-state behavior, Bonanomi and Morari15 reduced the number of variables by averaging the stage mole fractions of a three-component system. The dynamics of the averaged mole fractions are described in terms of a driving force based on the difference between the instantaneous bottom composition and its ultimate steady-state composition, causing either accumulation or depletion of the component. Vector fields were plotted containing detailed information about the dynamic behavior. 2.3. Scope of the Current Work. The focus of the present work encompasses the following: a. Description of a methodology in which multiple steady states can be determined for a homogeneous azeotropic distillation column using the inherent advantages of rigorous thermodynamic libraries and robust solvers in a process simulator without involving additional programming. b. Investigation of multiple steady states in real columns involving nonideal stages.

Figure 1. Representation of ∞/∞ analysis (after Bekiaris et al.5).

c. Illustration of how conventional and popular algorithms such as inside-out can be made to capture all the multiple steady states using stabilizing control. d. Finally, in a preliminary study, demonstration with open-loop dynamic simulations of how sustained oscillations on the higher branches reported by Morari and co-workers12-15 can be obtained, thereby extending the reach of process simulators for design, simulation, and analysis. Due to the relevance of the ∞/∞ analysis to the current study, it is described below. 3. The ∞/∞ Analysis (Bekiaris et al.5) This analysis is essentially based on graphical representation of the residue curves for a three-component system in triangular coordinates. A ternary 001 system16 that represents a minimum-boiling azeotrope (Az) between the heavy-boiling component (H) and the lightboiling component (L) is considered. Multiple steady states are shown to occur for this system when tracing the bifurcation path using the distillate flow rate as the parameter. In Figure 1 initially, the distillate composition is fixed at the azeotropic composition (Az) and the residue composition at the feed composition (F ) B), implying zero distillate output by the material balance and tie line concepts. The distillate flow rate is increased while maintaining the azeotropic distillate composition by moving the location of residue (B) from the feed location (F) toward the intermediate-component vertex (I) as shown in Figures 1a and 2. Subsequently, moving the material balance line in an anticlockwise fashion along the edges of the residue curve diagram (Figure 1), the location of the distillate is traced in succession from Az, along the L-H edge, to the light-component vertex (L), along the L-I edge, and, subsequently, into the triangle along B′H (Figure 1d) until the distillate flow rate equals the feed flow rate. This is sketched in Figure 1a-d in the form of the loci of movement of the material balance lines. These lines have the distillate

4388

Ind. Eng. Chem. Res., Vol. 44, No. 12, 2005 Table 1. Specifications for Case Study 1

Figure 2. Demonstration of MSS using ∞/∞ analysis.

Figure 3. Qualitative bifurcation diagram plotted in terms of the heavy component in the residue (xw,H):14 s stable; - - - unstable.

composition fixed at the upper end and the residue composition fixed at the lower end. The zone of multiple steady states falls within the distillate flow rates of D ) L and D ) M and is marked by turning points TP1 and TP2 at the ends (Figure 2). Though lines (b) and (c) overlap, implying the same values of compositions, they are shown to be separate for the sake of clarity. As discussed in section 2.2, Dorn and Morari14 have shown through dynamic simulations that, in addition to the unstable middle branch, even a section on the higher branch can become unstable with complex bifurcation behavior. The qualitative bifurcation diagram, based on dynamic simulations, is depicted in Figure 3. 4. Procedure Adopted for Finding Multiple Steady States Using HYSYS17 With steady-state simulations, when the distillate flow rate is chosen as the bifurcation parameter, its lower bound is nearly zero and its upper bound is nearly the feed flow rate. Outside the zone of multiplicity, distinct steady states exist. The bifurcation parameter close to the lower limit will possess a distinct type of steady state (type (a) profile in Figure 2), and that close to the upper limit will possess a different type of steady state (type (d) profile in Figure 2). 4.1. Case Study 1. The procedure is demonstrated through the case study based on Dorn et al.11 (Table 1). The system has a binary azeotrope for the light-heavy pair of methanol (L)-toluene (H) located at 77 mass % (or 88 mol %) methanol. In this system, the intermediate-boiling entrainer, methyl butyrate (I), separates the binary azeotropic constituents L-H. 4.1.a. Switching between Types of Multiple Steady States. The bifurcation plot involving the composition of the heavy component in the residue is

feed conditions

column specifications

state: saturated liquid pressure: 101.32 kPa feed flow rate: 3.58 kg/h composition (mass fractions): methanol (L): 0.66 methyl butyrate (I): 0.06 toluene (H): 0.28

pressure: 101.32 kPa number of stages: 40 (excluding condenser and reboiler) feed stage location: 21 from top type of condenser: total reflux rate: 13.50 kg/h thermodynamic models: liquid phase: Wilson vapor phase: ideal gas solver: sparse continuation

plotted against the distillate flow rate (Figure 3) and will be used as the basis for further discussion. 4.1.a.1. Switching to the Lower Concentration Steady State (LSS). This steady state is characterized by poor separation and low purity of products in the outlet streams. Temperature estimates alone are sufficient to switch to this stable steady state. Both the reboiler and condenser temperature estimates are specified close to the temperature corresponding to the L-H mixture’s azeotropic condition, implying negligible separation. The same temperature may also be specified at the feed stage. 4.1.a.2. Switching to the Upper Stable Steady State. This corresponds to the “higher steady state” (HSS) with reference to Figure 3. Again, the temperature estimates alone are sufficient to attain this type of steady state. The temperature estimate of the reboiler is specified to be higher than the temperature of the condenser. A reasonable guess would be the boiling point of the heavy component. The increased temperature takes into account the higher concentration of the heavy component in the bottom product, and consequently, the solution shifts to the high concentration steady-state branch. 4.1.a.3. Switching to the Unstable Steady State (USS). Column profiles of this type of steady state correspond to type (b) in Figure 2, where the top of the material balance relation moves along the L-Az edge of the triangle while the bottom of the material balance line correspondingly moves along the I-H edge. Specifying the temperature estimates alone was found to be insufficient for switching to the unstable steady state. Hence, composition profile estimates based on the ∞/∞ analysis may be used to switch to the unstable steady state. First, the location of a material balance line is assumed in the zone of unstable solutions. The composition profile in the column may be approximated by the residue curve passing through both the distillate (D) and residue (W) compositions located at the ends of the material balance line. However, locating such a curve may be difficult, involving composition estimates of the interior of the triangle. Hence, the composition profile estimates may be based along the edges of the triangle, with the condenser and reboiler compositions determined by the ends of the material balance line as illustrated by Bekiaris and Morari.6 The composition estimates for the condenser may be specified at the azeotropic composition or, preferably, close to the lightcomponent vertex. The feed tray is set to a binary mixture containing light and intermediate components. The reasons for this are given in section 5. The bottom product estimate is fixed at a binary mixture of heavy and intermediate species defined by the material balance line drawn from the azeotropic point through the feed composition and terminating on the I-H edge. Composition estimates for the remaining stages are


Figure 4. Initial-guess generation for the unstable profile and all three converged steady states (FS ) feed stage). Converged steady-state profiles: -O- higher; -]- unstable; -4- lower; s ∞/∞ initial guesses for the unstable profiles (interpolated by HYSYS); - - - initial guesses based on top and bottom plate specifications without specifying feed stage composition.

Figure 7. Liquid-stream mole-fraction profiles for the unstable steady state D ) 3.1 kg/h: -[- methanol; -9- methyl butyrate; -b- toluene.

Figure 8. Liquid-stream mole-fraction profiles for the low concentration steady state D ) 3.1 kg/h: -[- methanol; -9methyl butyrate; -b- toluene.

Figure 5. Temperature profiles for three steady states at D ) 3.1 kg/h: -[- higher steady state; -9- unstable steady state; -b- lower steady state.

Figure 6. Liquid-stream mole-fraction profiles for the high concentration steady state D ) 3.1 kg/h: -[- methanol; -9methyl butyrate; -b- toluene.

interpolated automatically by HYSYS. The initial guess generated is shown in Figure 4. The converged multiple steady states using the sparse continuation solver in HYSYS for the conditions given in Table 1 are shown in Figures 5-8. 4.1.b. Locating Turning Point 1 in the Bifurcation Diagram. The first turning point (TP1) as shown in Figure 3 is obtained by tracking solutions starting from a very low value of the distillate flow rate.

Specification of the lower limit should be the same as the currently converged value at some low specified value of the distillate flow rate (D). The upper limit should be set close to the value of the feed flow rate (F). When the distillate flow rate increases, the simulator tries to continue to stay at the same lower stable steady state. This is true even in the zone of multiplicity. As the zone of multiplicity passes, the present solution branch no longer exists. The program fails to converge and suddenly jumps to the higher steady state. This jump can be observed by plotting the distillate flow rate versus the chosen dependent variable, which is usually the composition of the light or heavy species in the distillate or residue. This highest value of the distillate flow rate at which the lower stable steady state continues to exist is designated as turning point 1 (TP1). Choosing smaller steps near the calculated value can further refine this value of D. Alternatively, the method based on norms can also be used as discussed below. 4.1.c. Locating Turning Point 2. The second turning point is obtained by tracking solutions starting from a very high value of the distillate flow rate (D). The value of D is decreased to a value below which it is not possible to continue with the same type of steady state. This lowest value of D is designated as TP2 (Figure 3). The distillate rates that lie within the bounds of TP1 and TP2 represent the zone of multiple solutions. To have a quick change over to different types of steady states, the earlier converged solutions should not be used as initial guesses. The norm of the solution based on all the state variables of the column, defined according to eq 1, was calculated for successive solutions (“n+1” and “n”) along a particular branch. It is noticed that the norm of the

4390

Ind. Eng. Chem. Res., Vol. 44, No. 12, 2005

Figure 9. Identification of turning points at the higher-steadystate (HSS) and the lower-steady-state (LSS) branches.

solution jumps to a very high value at the turning points (Figure 9).

norm ) ||Xn+1 - Xn||

Figure 10. Initial temperature profiles generated by different methods and the unstable solution: s final converged unstablesteady-state temperature profile; -O- initial temperature profile based on ∞/∞-based initial guess procedure; -0- initial temperature profile based on specifying only top and bottom compositions without specifying feed stage temperature (does not converge to any steady state). f Indicates direction of convergence.

(1)

HYSYS has different options for solvers that include inside-out algorithms,18 simultaneous correction procedures,19 and the sparse continuation method.17 While the inside-out methods could converge only to the lower or upper steady states, it is notable that the sparse continuation solver could converge to the unstable solution as well. Vadapalli and Seader,9 in an earlier study, illustrated the use of the Aspen Plus process simulator in tracking the multiple steady states. This required interfacing a user-programmed FORTRAN routine serving as a corrector to the Aspen Plus predicted values. 5. Analysis of the Proposed Methodology Convergence to the USS through steady-state simulations is facilitated through initial guesses that are close to the actual solution. However, the unstable solutions are not known apriori for a given system, and the ∞/∞ procedure facilitates their identification. The solver, based on these composition estimates, generates the bubble point temperature profile along the length of the column. It was observed that the composition estimates (as described in section 4 on switching to unstable steady states) in turn led to the generation of the temperature profile that was very similar to the final converged unstable-steady-state temperature profile (Figure 10). This was subject to further closer scrutiny. For the HSS (Figure 6), the entrainer composition peak occurs in the enriching section. The early occurrence of the entrainer peak permits the full development of the heavy-component profile, leading to the second sharp front in the temperature profile. However, the USS profile is distinctly characterized by the persistence of the light component up to the feed stage, pushing the entrainer maximum to the stripping section (Figure 7). The unstable profile is characterized by only one sharp temperature front (Figure 5), corresponding to the intermediate-boiling entrainer, and a diffused front, corresponding to the increasing amount of heavy species in the residue. The initial-guesses procedure for USS generation reflects these by providing compositions that are rich in the light component in the condenser, a mixture of light and intermediate components in the

Figure 11. Convergence pattern for inside-out algorithm using ∞/∞ procedure: s ∞/∞ initial guess; - - - intermediate iterations; -9- final converged higher steady state. f Indicates direction of convergence.

feed stage, and a mixture of heavy and intermediate species in the residue. The enriching section temperatures do not change much during subsequent iterations, while those in the stripping section move toward the final converged unstable steady state. On the other hand, if a linear variation based on only top and bottom product compositions had been provided as the initial composition estimates, the resulting temperature profile would not be similar to the final unstable-steady-state temperature profile and may not converge to the unstable steady state or even to the stable steady state. These composition estimates were also shown in Figure 4, while Figure 10 shows the initial and final converged temperature profiles. Further, the choice of the solver is also important. It is necessary for the solver to stay close to a given solution path if it is in the proximity of the final solution. A continuation algorithm (sparse continuation17), unlike the other conventional insideout algorithm,18 would tend to stay close to the estimated solution. Even when ∞/∞ guesses were provided to the inside-out algorithms, the enriching section temperature profile started to increase gradually with successive iterations, forming the first temperature front, and the second temperature front was formed in the stripping section. These led to the stable highersteady-state solution (Figure 11). Further, the simultaneous correction procedure was also tested and tended to capture only the higher stable steady state. The


inside-out algorithm was found to provide the lower and higher stable steady states depending upon the initial guesses. 6. Tracing of the Bifurcation Diagram 6.1. Tracing the Lower Stable-Steady-State Branch. The lower limit of the distillate flow rate is taken to be less than that obtained at TP2 (see Figure 3). The lower branch is traced until the first turning point TP1 as described in section 4.1.b. The distillate flow rate is defined as the independent variable, and the output variables such as the liquid-phase mole fraction of the heavy component in the bottoms are chosen as the dependent variables. 6.2. Tracing the Higher Concentration Branch. The simulation is converged at TP2. The same set of dependent variables as those selected for obtaining the lower-steady-state branch is chosen. The lower limit of the distillate flow rate in this case study is set equal to the value of TP2. The upper limit is set to be equal to a distillate rate greater than the value of TP1 and close to that of the feed rate to ensure that the complete higher-steady-state branch is generated. 6.3. Tracing the Unstable-Steady-State Branch. An unstable-steady-state solution, preferably at the midpoint of the two turning points, is obtained following the method presented under the section on switching to unstable steady states. Obtaining this midpoint is not difficult since the turning points were obtained accurately by the methods described above. This helps in reducing the influence of other steady states on rapidly converging to the unstable steady states on the path. The distillate rate is gradually reduced and convergence to the unstable-steady-state profile is attained at each value. To track the unstable profiles along the bifurcation path, the previous converged unstable solution may be used as the initial guess for the next distillate flow rate. This is continued with similar types of column profiles up to the value of D that is very close to TP2. Simulations in this unstable-steady-state case study are performed with the same set of dependent variables selected earlier. The lower limit of the distillate flow rate is set to be equal to TP2. The upper limit is selected to be slightly higher than the value of the distillate flow rate at TP1 in order to ensure that the unstable branch does not exist beyond TP1. The step size is specified, and the case study with positive increments is run. Step sizes were kept at 0.01 kg/hr. Executing the three steps described above generates the complete bifurcation diagram. The bifurcation diagram is shown in Figure 12. It may be observed that the turning points identified are in good agreement with those based on the values indicated by the norms. The method described above was implemented for a range of feed stage locations and feed compositions. The results are shown in Figures 13 and 14. 7. Validation of the Proposed Methodology (Case Study 2) The published case study of Vadapally and Seader9 was used for validating the proposed methodology. The specifications for this case study are given in Table 2. The typical results of the simulation obtained using the above methodology, shown in Figures 15 and 16, are in excellent agreement with the published results of Vadapally and Seader.9

Figure 12. Bifurcation diagram for the toluene (H) mole fraction in the bottoms stream (case study 1): ] LSS; O USS; 0 HSS.

Figure 13. Bifurcation diagrams for different feed stage locations (case study 1): -0- 5th stage; -]- 13th stage; -O- 21st stage; -4- 31st stage.

Figure 14. Bifurcation diagrams for feed streams with different methanol mass fractions (case study 1).

8. Detection of Multiple Steady States in an Azeotropic Distillation Column with Nonideal Stages The sparse continuation solver of HYSYS is suitable for generating the bifurcation diagrams for unity Murphree efficiency only. The inside-out solvers that can incorporate Murphree stage efficiencies cannot converge to the unstable steady states. To overcome this problem through steady-state simulations, the following procedure was adopted with reference to the conditions in Table 1. The AIChE method is used to estimate the efficiency of the sieve tray whose configurations are specified in Table 3. Perfect mixing and negligible entrainment on the trays was assumed. The profiles of the compositions

4392


Table 2: Specifications for Case Study 2 feed conditions

column specifications

state: saturated liquid pressure: 101.32 kPa feed flow rate: 2.0 kg/h composition (mass fractions): methanol (L): 0.6566 methyl-Butyrate (I): 0.0628 toluene (H): 0.2806

pressure: 101.32 kPa number of stages: 40 (excluding condenser and reboiler) feed stage location: 21 from top type of condenser: total reflux rate: 3.72 kg/h thermodynamic models: liquid phase: Wilson vapor phase: ideal gas solver: sparse continuation

in this case study were such that they divided the column into different sections. In each section, separation is found to take place effectively between the light and heavy key components, with the keys being different for each section. The pseudobinary approach as suggested by Lockett20 was adopted to estimate the Murphree tray efficiencies. After defining the pseudobinary components, the AIChE method outlined in King21 is used to predict the vapor phase Murphree efficiencies for binary systems. The AIChE method was coded separately. Options were provided for choosing the lengths of the section for applying the pseudobinary analysis as well as for selecting the light and heavy key components in the respective sections. The values of the Murphree efficiency were estimated to be approximately constant on most trays. An average value of Murphree efficiency of 70% was found to be suitable for the entire column over a wide range of distillate flow rates. The entire high concentration branch is generated with a gradually reducing number of trays until the hypothetical column closely matches the actual column (i.e., with EMV ) 70%). Equivalent performance when

Table 3. Specifications for Tray Sizing Based on the Column in Case Study 1 type of tray fractional approach to flooding foaming factor maximum tray diameter (cm) Weir height (cm) Weir length (cm) total area (m2) active area (m2) (down comer/total area) ratio (m2/m2) (sieve hole area/active area) ratio (m2/m2) tray spacing (cm) tray number with this maximum diameter

sieve 0.80 0.85 11.0 5 8.0 0.0095 0.0076 0.10 0.12 20.0 41

comparing the two columns was defined in terms of matching exit concentrations as well as intermediate profiles. The number of ideal stages inside the column which gave equivalent performance for this case was estimated by trial and error to be 28. The feed stage was maintained at the same relative position. The higher-concentration branch was generated for a 40-tray column with a Murphree efficiency of 70% as well as its equivalent column with 28 ideal stages. The results are very close, as can be seen in Figure 17. Due to the small difference in the compositions despite a tray difference of 12, the HSS branch with 100% column efficiency also appears to be close to the other two profiles, especially at high distillate flow rates. But there is a noticeable difference for lower values of the specified distillate flow rates for the USS and HSS branches. Now, assuming that the column operating in the unstable region would also behave with approximately 70% overall efficiency, the unstable profiles were generated using the sparse continuation algorithm with 28 stages as well. The unstable-steady-state branch and lower-steady-state branch were also generated using procedures outlined earlier. The entire bifurcation diagram is illustrated in Figure 17. It is noted that the first turning point characterizing the lower-steady-state branch is located at nearly the same distillate flow rate for both the ideal and real columns while a significant deviation occurs in the location of the second turning point. The extent of branching of the unstable steady state and the range of distillate flow rate covered by the zone of multiplicity diminishes with a decreasing number of stages. For comparison with the results of Dorn et al.11, the turning points for a case involving a 40-tray bubble cap column (with a Murphree efficiency

Figure 15. Bifurcation diagram for the methanol mole fraction in the distillate stream (case study 2): 4 LSS; ] USS; O HSS.

Figure 16. Bifurcation diagram for the methyl butyrate mole fraction in the bottoms stream (case study 2): 4 LSS; ] USS; O HSS.

Figure 17. Bifurcation diagram showing the effects of Murphree efficiency. Efficiency ) 1: O HSS; 4 USS; ] LSS. Murphree efficiency: 0 0.70 HSS (HYSIM inside-out). Equivalent 28 stages: s LSS; 2 USS; b HSS.

Ind. Eng. Chem. Res., Vol. 44, No. 12, 2005 4393 Table 4. Specifications for Case Study 3 feed characteristics and column detail condition pressure (kPa) flow rate (kmol/h) composition (mole fractions): acetone (L) benzene (I) n-heptane (H) number of stages type of condenser pressure (kPa) feed stage location (from top) reflux rate (kmol/h) thermodynamic models: liquid phase vapor phase column solver

feed 1

feed 2 (entrainer)

saturated liquid 101.32 100

saturated liquid 101.32 1

0.9123 0.0 0.0877

0 1 0

42 (exclusive of condenser and reboiler) total 101.32 21 180 UNIQUAC Ideal gas sparse continuation

of 0.5) were 3.05 and 3.25 kg/h distillate flow rates. Here, for a tray column with an efficiency of 0.7, the turning points were found to be ∼2.86 and 3.26 kg/h. The higher the efficiency, implying more stages, the broader is the range of multiplicity. Hence, the second turning point in this work occurred at a lower distillate flow rate than those of Dorn et al.11 For comparison with the ∞/∞ plots of Dorn et al.,11 the turning points occur at 2.58 and 3.27 kg/h, while in this work for the ideal 40-stage column, they occur at 2.66 and 3.26 kg/h. 9. Testing the Proposed Methodology for Different Systems 9.1. Case Study 3. The proposed methodology involving the ∞/∞ procedure of generating initial guesses and the sparse continuation solver for detecting multiple steady states was tested for the system acetone (L)benzene (I)-n-heptane (H). The system studied considers the separation of the acetone-n-heptane mixture using an intermediate entrainer benzene. The specifications for this case study are given in Table 4. The initialguess generation for this case study was also based on the ∞/∞ methodology and the bifurcation diagram is given in Figure 18. Further, the methanol (L)-n-butyraldehyde (I)toluene (H) system6 was also analyzed, and all the multiple steady states could be captured using the methodology given above. The results are given as part of the Supporting Information. The thermodynamic models and their binary interaction parameters from HYSYS are provided as Supporting Information for these systems.

Figure 18. Bifurcation diagram based on the n-heptane (H) mole fraction in the bottoms stream for case study 3: ] LSS; 4 USS; O HSS.

efficiency concept. The methodology for rigorously incorporating Murphree efficiencies and converging to the unstable steady state using the inside-out solver would require dynamic simulations with stabilizing control as described in the following section. 11. Procedure Adopted for Stabilizing the Unstable Steady State in Ideal Staged Columns The concept of stabilizing control using the PI controller as outlined by Dorn et al.11 was adopted to capture the unstable-steady-state solution. The multiple-steadystate regions are identified using the procedures described in section 4, and temperature maps as a function of distillate flow rate are plotted. The column is first simulated corresponding to the higher-steady-state (HSS) region in the dynamic mode. The temperature on a particular tray that corresponds to the HSS region is specified for applying stabilizing control. On this basis, the dynamic simulation will eventually settle to the distillate flow rate corresponding to the higher-steadystate region. Next, the set point temperature on the chosen tray is changed to a new value from the temperature map such that the simulator will still stay on the HSS branch but at a lower distillate flow rate. This is continued until the temperature for stabilizing control is changed to that corresponding to the unstable steady state. The unstable-steady-state temperature value was determined from the map produced by the steady-state sparse continuation solver. Now the distillate flow rate that was decreasing with successive dynamic simulations as long as the higher-steady-state branch was traced will start to increase. In other words, the converged solution is in the unstable-steady-state mode for a distillate flow rate for which a higher stable-steadystate solution was detected earlier. 12. Case Study 4

10. Limitations of Solvers The approach described above may be recommended for solely steady-state process simulators. The limitations of the commonly used and robust inside-out algorithms for detecting multiplicities were demonstrated and the use of the sparse continuation solver in HYSYS17 to capture the unstable steady states through suitable initial guesses was illustrated. The inside-out algorithms were able to incorporate Murphree efficiencies but could not detect unstable solutions, whereas the sparse continuation solver could converge to unstable solutions but could not incorporate efficiencies. The latter limitation was overcome through the overall

The published case study from Dorn et al.,11 which has the added advantage of experimental data, is referred to in this work. However, the column efficiencies are set to unity and a different control scheme will be used as discussed in section 13. The feed conditions and column specifications are according to Table 1. For the present case study, the multiple-steady-state temperature maps are shown in Figure 19. This diagram would provide guidelines to the temperature estimates that can be fixed as the set point on a suitably chosen tray for implementing the dynamic simulations. The dynamic simulation will eventually settle to the corresponding steady state, which may be either stable

4394


Figure 19. Temperature map as a function of stage number and distillate flow rates.

Figure 21. Control configuration of the distillation column.

13. Control Scheme for Dynamic Simulations

Figure 20. Comparison of steady-state and dynamic simulation temperature profiles for the distillate flow rate of 3.1 kg/h (HSS ) higher steady state and USS ) unstable steady state): [ dynamic simulation; s steady-state simulation (HSS); 2 dynamic simulation; - - - steady-state simulation (USS).

or unstable depending on the temperature set point value. The corresponding distillate flow rate at steady state is noted. For validation purposes, a distillate flow rate of 3.1 kg/h is selected both in the higher- and unstable-steady-state regions. From Figure 19, the corresponding higher-steady-state (HSS) temperature on the 10th tray is noted (viz. 101.35 °C). This temperature is maintained on the 10th tray by stabilizing PI control in the dynamic simulations. After reaching the steady state, the temperature profiles from the dynamic simulations are compared to those from the steady-state simulations. The temperature profiles from the steadystate and dynamic simulations are in excellent comparison as seen from Figure 20. Similarly, the temperature value corresponding to the unstable steady state (103.3 °C) is fixed on the 31st tray. The solutions converged to the unstable profile and are again in good comparison with the predictions of direct steady-state simulations (Figure 20). The steady-state simulation procedures were described in sections 4 and 6. It is interesting to note that, in both dynamic simulations, the HYSIM inside-out solver could be used without any problems. As noted earlier, if direct steady-state simulations had been carried out for the distillate specifications of 3.1 kg/h, the HYSIM inside-out solver would have converged to the higher steady state only while the sparse continuation solver would have required suitable initial guesses to converge to the unstable-steady-state profiles.

A different scheme to that of Dorn et al.11 is used in the current study (Figure 21). In their work, the distillate flow rate was used to control the temperature of the chosen tray, while the condenser and reboiler liquid levels were controlled by the reflux flow rate and bottoms flow rate, respectively. The column operated at atmospheric pressure. In the present study, the reboiler heat input was varied to control the temperature on a particular location, as it was found to give a comparatively faster response. The liquid levels in the reboiler and condenser were controlled by manipulating the residue and distillate flow rates, respectively, while the cooling water flow rate was used to control the pressure inside the column. PI controllers were used. 14. Effect of Murphree Efficiencies Incorporating Murphree efficiencies in columns exhibiting multiple steady states introduces additional complications. These include the reduction in the region of multiplicities and the unavailability of the temperature map that was generated for ideal columns as described in section 11. In other words, the temperatures that should be specified for stabilizing control on the trays are unknown. Further, as discussed in section 10, the solvers have limitations for incorporating efficiencies and tracking unstable steady states. The following methodology was adopted to get around these limitations. Using the inside-out solver, the entire higher-steadystate branch (HSS) was generated. The concept of overall efficiency was applied for tracing the unstable profile. If the number of theoretical stages that gives equivalent separation performance as the actual column can be identified, then the temperature map over the range of distillate flow rates can be easily generated. This will include the unstable-steady-state solutions as well by using the sparse continuation solver. With this map in hand, the procedure described in section 11 for dynamic simulations may now be executed for the actual column with Murphree efficiencies. A distillate flow rate in the unstable-steady-state region is chosen for the equivalent, ideal column, and for this specification, the


Figure 22. Temperature map for a 28-stage column (inside-out algorithm for stable steady states and sparse continuation solver for unstable steady states).

Figure 23. Temperature map for a 40-stage column with 70% Murphree efficiency (inside-out algorithms used for all steady states).

unstable-steady-state temperature profile is known. One temperature in the stripping section is selected, as it is known that the steep temperature front is occurring in this section in the unstable-steady-state region. This temperature was assigned as the set point to one of the appropriate stages in the stripping section of the actual column with Murphree efficiencies. The dynamic simulations were carried out with stabilizing control and the inside-out algorithm that permits efficiencies to be incorporated in the trays. The final distillate flow rate attained was compared with the distillate flow rate specification of the ideal-staged, equivalent column. If the distillate flow rates from the dynamic simulations did not match with the equivalent, ideal-stage column, then the temperature location was changed to a nearby tray. By trial and error involving only a few iterations, the correct stage for the temperature specification that matched very closely with the distillate flow rate specification was determined. On the basis of the specifications in case study 1, a column with a Murphree efficiency of 70% was calculated. For this column, it was found that an equivalent ideal-staged column with 28 stages gave the closest separation performance. The temperature maps for the 28-ideal-staged column are shown in Figure 22. Using this and the procedure described above, all the steady states, including those in the region of multiplicity, could be identified using the inside-out solver with a Murphree efficiency of 0.7. The resulting temperature

Figure 24. Comparison of ideal-staged column with real column with EMV ) 70%: s 40 stages, EMV ) 70% (steady state); 2 40 stages, EMV ) 70% (dynamic); - - - 28 ideal stages (steady state).

Figure 25. Comparison of simulation results with the experimental results of Dorn et al.:11 - - - simulation results; b experimental results.

map for the 40-stage column is shown for comparison in Figure 23. In Figure 24, the two bifurcation diagrams are compared and observed to be in reasonable agreement. As seen from Figure 17, the nonideal, 40-stage column has a slightly different turning point (TP2) from that of its equivalent column with 28 ideal stages. As explained above, the temperature maps are based on the nonideal column on the higher branch and the equivalent, ideal column on the middle branch. This accounts for the bump which occurs in the transition region between the unstable and higher-stable-steadystate regions near the turning point in Figure 24. The same procedures were repeated for the experimental results of Dorn et al.,11 where the Murphree efficiency was reported to be 0.5 and the comparison is given in Figure 25. 15. Stability of the Solutions Dynamic simulations using HYSYS also enable the demonstration of the stability of the different branches. Among the three steady states, the higher and lower steady states are open-loop stable and the unstable steady state is open-loop unstable. The response given in terms of the averaged temperature of the distillation column in the HSS region for feed disturbance ((10% for a duration of 20 min) is shown in Figure 26. For the positive feed flow rate disturbance for a fixed period of time, it was observed that the average temperature of the column initially changed and subsequently regained

4396


Figure 26. Dynamics of the temperature profile in the higher steady state for a (10% disturbance in the feed. The distillate flow rate is 3.1 kg/h.

Figure 27. Dynamics of methanol composition in the distillate for a (10% disturbance in the feed in the HSS region for a distillate flow rate of 3.1 kg/h.

Figure 29. Dynamics of temperature profile corresponding to the USS indicating instability for a +10% disturbance in the feed at the distillate flow rate of 2.9 kg/h.

Figure 30. Dynamics of methanol composition in the distillate for a +10% disturbance in the feed for the USS at a distillate flow rate of 2.9 kg/h.

Figure 28. Dynamics of toluene composition in the bottoms for a (10% disturbance in the feed corresponding to the HSS for a distillate flow rate of 3.1 kg/h.

Figure 31. Dynamics of toluene composition in the bottoms for a +10% disturbance in feed corresponding to the USS at a distillate flow rate of 2.9 kg/h.

its original value. The responses to the feed disturbances in the product compositions (methanol in the distillate and toluene in the bottom) are shown in Figures 27 and 28. From these observations, it may be concluded that the distillation column is open-loop stable. The system open-loop dynamics are such that nearly 12 h are required to bring the column back to its original steady state, resulting in the production of large quantities of off-specification product. The temperature response of the distillation column in the USS region for a feed disturbance (+10% for a duration of 20 min) is shown

in Figure 29. For a feed disturbance, the system exhibits oscillatory behavior, indicating the unstable nature of the solution. Figures 30 and 31 shows the oscillatory behavior of methanol and toluene composition in the distillate and bottoms streams in response to feed disturbances. 16. Closed-Loop Behavior in the Region of Multiplicity for Feed Disturbances Figure 32 shows the temperature-profile behavior in the HSS region for feed disturbances. In both cases, the


Figure 32. Closed-loop temperature-profile behavior in the HSS region for feed disturbances at a distillate flow rate of 3.1 kg/hr.

Figure 33. Closed-loop temperature-profile behavior in the USS region for feed disturbances at a distillate flow rate of 2.9 kg/h. Table 5. Specifications for Case Study 5 feed conditions temperature: 76.84 °C pressure: 163.8 kPa feed flow rate: 100 kmol/h composition (mole fractions): methanol (L): 0.8671 methyl butyrate (I): 0.005 toluene (H): 0.1279

column specifications pressure: 151 kPa (condenser) 165.8 kPa (reboiler) number of stages: 44 (exclusive of condenser and reboiler) internals: sieve tray feed stage location: 40 (from top) type of condenser: total reflux rate: 2000 kmol/h distillate rate: 97.5 kmol/h tray holdup: 9 kmol (on all trays) thermodynamic models: liquid phase: Wilson vapor phase: ideal gas

temperature returns to its original state after the feed rate is brought back to its original value. In the closedloop mode, the system returns to its original state much faster when compared with the open-loop case. Figure 33 shows the closed-loop temperature-profile behavior in the USS region for feed disturbances. For a positive feed disturbance, initially the temperature is found to decrease and, after some time, it returns to its original value. For a negative feed disturbance, initially the temperature is found to increase. In both cases, the temperature is showing oscillatory behavior before returning to its original state. The reason for this is that, because of the instability of the solution, it will try to go to either the lower steady state or higher steady state for the feed disturbance. For positive feed disturbances, it will try to go to the lower stable steady state, since during this dynamic the temperature along the length

Figure 34. Time history of methyl butyrate composition on tray 11 with subsequent step changes in D′.

Figure 35. Transition of the methyl butyrate composition profile on tray 11 from damped oscillations to sustained oscillations.

of the column decreases. Due to the control action, it returns to the original steady state. A similar phenomenon was observed for a negative feed disturbance. In this case, the column tries to go to the HSS region, where the temperature increases along the length of the column. 17. Open-Loop Dynamic Simulations In a preliminary study reported for the first time with HYSYS, open-loop dynamic simulations were also carried out to obtain the limit cycle behavior on the higher stable branch. Details of the case study are given in Table 5. Perturbations to the distillate-to-feed ratio (D′) were given in small decrements typically of 0.001 starting from 0.975. Typical results are shown in Figures 34-36. In these figures, damped oscillations are observed on the higher branch (D′ > 0.9325). At D′ ) 0.9325, sustained and bounded oscillations begin to occur, increasing in amplitude with subsequent step changes. Continuing the perturbations beyond D′ ) 0.927 results in the profiles suddenly moving to the lower steady state with a consequent disappearance of oscillations (Figure 36). It is interesting to note that Dorn et al.,22 in a different study involving the acetone (L)-benzene (I)n-heptane (H) system, obtained similar behavior in a narrow region where there was only one solely attracting stable steady state and two unstable steady states. Additional plots on the oscillations in the compositions are provided in the Supporting Information. Lee et al.12 observed the limit cycle behavior for a light-to-heavy mole fraction ratio in the feed (xFLH) ) 6.78, a reflux-

4398


Figure 36. Methyl butyrate composition on tray 11 as the system immediately jumps to the lower steady state at D′ ) 0.9265.

to-feed ratio (L′) ) 12, and D′ between 0.9641 and 0.9735. Both the constant molar overflow (CMO) model and Aspen Plus simulator (RADFRAC option) were used by these authors with qualitatively similar but quantitatively different results. In our present simulations also, while the column configuration is similar, the holdup, pressures and pressure drop, feed conditions, column diameter, and L′ are different. Experimental investigation of the interesting effect of limit cycle behavior on the higher branch is recommended. Further, a detailed investigation on the equipment parameters, such as feed stage location, holdup, column diameter, pressure drop, and tray details, is also required. 18. Conclusions The detection of multiple steady states, especially the unstable branch in a homogeneous azeotropic distillation column, is not a trivial task. A methodology based on the ∞/∞ approach5 was developed to capture the unstable steady state. The method was validated with published case studies and was observed to work with different systems, feed stage locations, numbers of stages, and feed compositions. Dynamic simulations with stabilizing control were also carried out to overcome the shortcomings of conventional solvers in tracking the unstable steady state. The difficulty involving Murphree efficiencies was overcome using the overall efficiency approach. Both open- and closed-loop responses of the column to feed disturbances were illustrated for higher and unstable steady states. As a preliminary investigation, the occurrence of sustained oscillations on the higher branch was demonstrated with HYSYS. The inherent advantages of commercial process simulators, if combined with the methodology of detecting multiple steady states without additional programming, will be of great use to both research and industrial practice. The issues of operation and control in the multiple-steady-state regimes can be easily handled in the industries. From a theoretical viewpoint, the rigorous analysis of the highly nonideal and nonlinear nature of the homogeneous azeotropic distillation column represents a formidable challenge. Quickly identifying the multiple steady states and their regions of occurrence will be useful to test simpler models that may be developed to understand the nonlinear behavior of

multicomponent, multistage chemical processes. It is also recommended that simulators provide tools for bifurcation analysis owing to their increasing importance in the simulation and control of separation processes. Acknowledgment The partial support provided by the Ministry of Human Resources Development, India, under the R&D scheme to implement this work is gratefully acknowledged. Supporting Information Available: For the systems (I) methanol (L)-methyl butyrate (I)-toluene (H), (II) acetone (L)-benzene (I)-n-heptane (H), and (III) methanol (L)-n-butyraldehyde (I)-toluene (H), the information provided includes the liquid-stream molefraction profiles for systems II and III; the feed and column conditions for system III; the binary interaction parameters for systems I, II, and III; the Wilson and UNIQUAC models from the HYSYS documentation;17 and additional plots on the oscillations of the compositions, described in section 17. This material is available free of charge via the Internet at http://pubs.acs.org. Nomenclature B ) residue D ) distillate D ) distillate flow rate (kmol/hr) D′ ) ratio of distillate-to-feed molar flow rate E ) Murphree efficiency H ) heavy component I ) intermediate component L ) light component L′ ) ratio of reflux-to-feed molar flow rate HSS ) higher steady state LSS ) lower steady state TP1 ) turning point 1 TP2 ) turning point 2 USS ) unstable steady state X ) state variable (equation 1) x ) mole fraction Superscripts F ) feed

Ind. Eng. Chem. Res., Vol. 44, No. 12, 2005 4399 Subscripts B ) residue H ) heavy species MV ) Murphree vapor LH ) light-to-heavy component Symbols ∞ ) infinity

Literature Cited (1) Kienle, A.; Marquardt, W. Bifurcation Analysis of Multicomponent, Non-equilibrium Processes. Chem. Eng. Sci. 1991, 46, 1757. (2) Gani, R.; Jorgensen, S. B. Multiplicity in Numerical Solution of Nonlinear Models: Separation Processes. Comput. Chem. Eng. 1994, 18, S55. (3) Magnussen, T.; Michelsen, M. L.; Fredunslund, A. Azeotropic Distillation Using UNIFAC. In Third International Symposium on Distillation. I. Chem. E. Symposium Series 56; The Institution of Chemical Engineers: London, 1979; p 1. (4) Jacobsen, E. W.; Skogestad, S. Multiple Steady States in Ideal Two Product Distillation. AIChE J. 1991, 37, 499. (5) Bekiaris, N.; Meski, G. A.; Radu, C. M.; Morari, M. Multiple Steady States in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1993, 32, 2023. (6) Bekiaris, N.; Morari, M. Multiple Steady States in Distillation: ∞/∞ Prediction, Extensions, and Implications for Design, Synthesis, and Simulation. Ind. Eng. Chem. Res. 1996, 35, 4264. (7) Guttinger, T. E.; Morari, M. Predicting Multiple Steady States in Equilibrium Reactive Distillation. 2. Analysis of Hybrid Systems. Ind. Eng. Chem. Res. 1999, 38, 1649. (8) Guttinger, T. E.; Dorn, C.; Morari, M. Experimental Study of Multiple Steady States in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1997, 36, 794. (9) Vadapalli, A.; Seader, J. D. A Generalized Framework for Computing Bifurcation Diagrams Using Process Simulation Programs. Comput. Chem. Eng. 2001, 25, 445.

(10) Aspen Plus v. 10 User Guide; Aspen Technology, Inc.: Cambridge, MA, 1999. (11) Dorn, C.; Guttinger, T. E.; Wells, G. J.; Morari, M.; Kienle, A.; Klein, E.; Gilles, E. D. Stabilization of an Unstable Distillation Column. Ind. Eng. Chem. Res. 1998, 37, 506. (12) Lee, M.; Dorn, C.; Meski, G. A.; Morari, M. Limit Cycles in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1999, 38, 2021. (13) Dorn, C.; Morari, M. Qualitative Analysis for Homogeneous Azeotropic Distillation. 1. Local Stability. Ind. Eng. Chem. Res. 2002, 41, 3930. (14) Dorn, C.; Morari, M. Qualitative Analysis for Homogeneous azeotropic Distillation. 2. Bifurcation Analysis. Ind. Eng. Chem. Res. 2002, 41, 3943. (15) Bonanomi, E.; Morari, M. Phase-Plane Analysis for Homogeneous Azeotropic Distillation Columns. Ind. Eng. Chem. Res. 2002, 41, 3963. (16) Matsuyama, H.; Nishimara, H. Topological and Thermodynamic Classification of Ternary Vapour-Liquid Equilibria. J. Chem. Eng. Jpn. 1977, 10, 181. (17) HYSYS v. 3.0.1 User Guide; Hyprotech, Ltd.: Calgary, Canada, 2002. (18) Boston, J. F.; Sullivan, S. L., Jr. A new class of solution methods for multicomponent, multi-stage separation processes. Can. J. Chem. Eng. 1974, 52, 52. (19) Naphthali, L. M.; Sandholm, D. P. Multicomponent separation calculations by linearization. AIChE J. 1971, 17, 148. (20) Lockett, M. J. Distillation Tray Fundamentals; Cambridge University Press: New York, 1986. (21) King, C. J. Separation Processes; McGraw-Hill Book Company: New York, 1980. (22) Dorn, C.; Lee, M.; Morari, M. Stability and transient behavior of homogeneous azeotropic distillation. Comput. Chem. Eng. 1999, 23 (Suppl.), 185.

Received for review June 25, 2004 Revised manuscript received March 29, 2005 Accepted April 8, 2005 IE049443S

Multiple-Steady-States Identification in Homogeneous Azeotropic

Recommend Documents