Multiple Steady-States Analysis and Unstable ... - ACS Publications

Article pubs.acs.org/IECR

Multiple Steady-States Analysis and Unstable Operating Point Stabilization in Homogeneous Azeotropic Distillation with Intermediate Entrainer Weisong Li,†,‡,§ Lei Zhong,†,‡,§ Yongchao He,† Jihong Meng,† Fanglian Yao,† Yusheng Guo,† and Chunjian Xu*,†,‡,§ †

School of Chemical Engineering and Technology, and ‡State Key Laboratory of Chemical Engineering, Tianjin University, Tianjin 300072, China § Collaborative Innovation Center of Chemical Science and Engineering, Tianjin 300072, China S Supporting Information *

ABSTRACT: The steady-state multiplicity in azeotropic distillation processes has been theoretically predicted and verified experimentally for decades. However, it is a tough task to detect all the steady-state solutions. The purpose of this article is 2-fold: first, it provides a systematical study on the multiplicity in homogeneous azeotropic distillation with intermediate entrainer and a novel strategy is proposed to find all the steady-state solutions using “Design Specs/Vary” in Aspen Plus. Operating on the unstable steady-state branch benefits both energy saving and reducing capital investment with the product quality specified. Because of the great benefits that operating on the unstable steady-state branch offers, the second purpose is to present suitable control schemes to stabilize the operating point on the unstable steady-state branch to achieve robust control. However, finding all the steady-state solutions for distillation is an extremely tough task. Many methods have been proposed to identify the multiple solutions in distillation, including graphical approach, homotopy continuation methods, and simulation approach.4,5,9,17−20 The multiplicities in ternary azeotropic distillation separation sequences were studied in depth by Güttinger21 and Esbjerg,22 respectively. Güttinger21 pointed out that a binary azeotrope (minimum or maximum binary azeotrope) can be separated using distillation with an entrainer whose boiling point is between the boiling points of the two components and does not show any singular points on the two boundaries connecting it to the pure components. This kind of distillation is called azeotropic distillation with intermediate entrainer, and a classic example is the separating of acetone/ n-heptane azeotropic mixtures with benzene as an intermediate entrainer. According to Matsuyama,23 the acetone/benzene/ n-heptane mixture belongs to the 001 class and can be separated in three possible distillation sequences. The three possible sequences are direct split sequence, indirect split sequence, and separation in stand-alone column, respectively. Bekiaris et al.4 showed that, in the 001 class, multiple steady states exist in a single column for any ternary feed composition. Thus, if the entrainer recycle is not closed, multiple steady states are expected for the first column of the direct or the indirect split sequence. Güttinger’s work21 proved that the existence of multiple steady states in intermediate entrainer separation sequence does neither depend on the external feed composition nor on the location of the azeotrope. Aspen Plus was utilized by Güttinger21 to detect

1. INTRODUCTION It has been long assumed that there is only one steady-state solution for a distillation process.1 Later, the existence of multiple steady states was predicted through theoretical analysis, and it was verified by experiments. The multiplicity in ternary homogeneous systems was discovered by Petlyuk2 for the first time under the assumption of constant molar flow. Magnussen’s3 simulation showed similar results for the heterogeneous mixture of ethanol/water/benzene. Bekiaris and co-workers4,5 performed a comprehensive study on multiple steady states for ternary homogeneous and heterogeneous mixtures with a powerful analytic tool, which is ∞/∞ analysis. The so-called ∞/∞ analysis refers to the analysis of a column operating under the extreme case where the number of trays (or packing height) and the reflux are infinite (the ∞/∞ case hereafter). The existence of multiple steady states was experimentally verified in a pilot-plant experiment for the methanol/methyl butyrate/toluene ternary system in 1997.6 Almost in the same period, the multiplicity in heterogeneous azeotropic distillation such as ethanol dehydration process and isopropyl alcohol−cyclohexane−water system was also verified in laboratory scale experiments.7,8 Additionally, the steady-state multiplicities in nonlinear processes such as reactive distillation have also been intensively studied in recent years.9−14 Gani and Jørgensen15 concluded three types of steadystate multiplicities, namely, output, input, and internal state multiplicity. In this article, the term “multiple steady states” (MSS) refers to output multiplicity. The implications of steadystate multiplicity for column design, synthesis (e.g., selection of entrainer, and the separation sequence) and simulation has been thoroughly discussed in Bekiaris and Morari’s16 work. They claimed that detecting only one solution may result in misleading conclusions, and some even attractive solutions may be missed when steady-state multiplicity exists. © XXXX American Chemical Society

Received: February 9, 2015 Revised: June 17, 2015 Accepted: July 16, 2015

A

DOI: 10.1021/acs.iecr.5b00572 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research the multiple solutions in the acetone/benzene/n-heptane ternary system separation sequence. In Güttinger’s work, the higher and lower stable steady states (HSS and LSS, respectively) have been traced by rigorous process simulations, but the unstable steadystate solutions were still missed. However, the existence of unstable steady states (USS) has already been experimentally verified by Dorn24 in 1998. As mentioned above, detecting all the steady-state solutions is of great importance. Recently, Baharev25 proposed a very brilliant globally convergent method for finding all the steady-state solutions of distillation columns in which initial estimates are not required and the only relevant aspect is that the Jacobian can be permuted to lower block Hessenberg form. The robustness of this method was confirmed by the successful solution of a numerically challenging reactive distillation column with seven steady states. The existence of MSS also poses a challenge to distillation control. Bekiaris et al.4,16 pointed out that when MSS existed, it was possible to jump from one steady state to another even for small disturbances. The transition between steady states has already been experimentally verified by Mohl et al.13 In some cases, operating on the USS solution branch might have considerable economic advantages. Thus, how to stabilize the unstable operating point might be of great importance. Through in-depth analysis, Jacobsen and Skogestad26 found that it was impossible to achieve the certain separation corresponding to unstable operating points under manual operation. The unstable operating point might be stabilized by feedback control of product composition or tray temperature (one-point control). During the experimental verification of the unstable steady state for methanol/methyl butyrate/toluene system, Dorn24 achieved the stabilization of the unstable operating point with PI control, but a high-purity product was not obtained. For an isopropyl alcohol/cyclohexane/water heterogeneous azeotropic column, Chien et al.27 proposed an inverse double-loop temperature control strategy to stabilize the unstable operating point. Although a high-purity product with minimum energy consumption and maximum product recovery could be obtained, the operating point should be located in the higher reflux region close to critical reflux. In this article, based on the ∞/∞ analysis of homogeneous azeotropic distillation for separating acetone/n-heptane mixtures with benzene as an intermediate entrainer presented by Güttinger,21 a novel strategy is proposed to detect all the steadystate solutions using the “Design Specs/Vary” function in Aspen Plus. This strategy takes a different aspect from the method proposed by Baharev.25 Then, the economic potentials of operations at the branches of stable and unstable steady state solutions are discussed. Further, the stabilization of the unstable operating point is explored by feedback control structures.

Figure 1. Direct split sequence referring to the sequences proposed by Güttinger21

italics (F, D, and B) refers to the corresponding flow rates. Ex refers to the external feed which is a combination of the feed and entrainer makeup (see Figure 1), and similarly Ex refers to the external feed flow rate. xpL, xpI, and xpH refer to corresponding compositions in stream P. xAZL and xAZH refer to corresponding compositions of the L−H azeotrope. Note that all the compositions are expressed in mass basis. The NRTL activity coefficient model is used to describe the nonideality of the acetone/benzene/n-heptane ternary system. Acetone and n-heptane form a minimum-boiling azeotrope at atmospheric pressure (101.3 kPa) with xAZL = 88.02 wt %, and the feasibility of direct split sequence for separating this ternary mixture is presented in Figure 2. The heuristic optimal design of

Figure 2. Feasibility analysis of the direct split sequence with residue curve map (RCM).

this direct split sequence can be determined by the following equations. At the design point, pure L and H are obtained. Thus, the corresponding product flow rates D1 and B2 (Figure 1) are given by the amounts of these components in the external feed flow (see eq 1).

2. STEADY-STATE ANALYSIS AND SIMULATION As mentioned above, there are three possible sequences for the intermediate case, namely, the separation in stand-alone column, the direct split, and the indirect split. In this article, the direct split sequence proposed by Güttinger21 (shown in Figure 1) is chosen to separate the acetone/n-heptane mixture with benzene as intermediate entrainer for case study and an entrainer makeup stream is added to balance the tiny entrainer loss in the two product steams. The first column and the second column are named as C1 and C2 hereafter. Unless otherwise specified, L (I and H) refers to the component which has the lowest (intermediate and highest) boiling point, respectively. The feed, distillate and bottoms of a column in the flowsheet are denoted by F, D, and B, respectively, and the same notation in

D1 = x EXL·E x

and

B2 = x EXH·E x →

x B2 = EXH = C (constant) D1 x EXL

(1)

Using eq 1 and the overall mass balance Ex = D1 + B1 obtains the form shown in eq 2. D + B2 D (1 + C) Ex = 1 = 1 F1 F1 F1 B

(2) DOI: 10.1021/acs.iecr.5b00572 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Using the lever rule gives the B1/D1 ratio which can be determined by the distances D1F1and B1F1.

B1 DF = 11 D1 B1F1

realistic conditions, and the two stable solutions on the LSS and HSS, respectively, were found. However, the USS solution was still missed. On the basis of the ∞/∞ analysis, a strategy for identifying multiple steady states using Aspen Plus 2006 version is proposed here. After specifying the components, property method, columns configuration (number of stages, feed stages, types of condenser and reboiler), pressure profiles, and operating parameters (D1, D2, R1 and R2), the following strategies are implemented. 2.2.1. Tracing the LSS Branch and HSS Branch. As stated by Bekiaris,15 using profiles predicted by the ∞/∞ analysis as initial estimates for simulations, there is a better chance to find all the solutions. According to the ∞/∞ analysis results presented above, we know that the LSS branch is characterized by poor separation, and the compositions of two column distillates D1 and D2 are equivalent to L−H azeotropic composition. Then, the L−H azeotropic temperature and composition at operating pressure are supplied on the upper stages of the two columns. Using the sensitivity analysis tool in Aspen Plus, starting from D1 = 0, which is located outside the region of multiple steady states, D1 increases continuously with small steps trying to reach the first turning point (TP1). When the simulation fails to converge or the converged solution differs significantly from that of the previous step, the region of multiplicity passes. Then, choosing smaller steps near TP1 can further refine TP1. On the basis of the simulated results, the changes in each stream composition with the variation of D1 can be plotted, and the LSS branch in the bifurcation diagrams is traced. For tracing the HSS branch, from the ∞/∞ analysis results, it can be found that there is a large segment of this branch on which stream B1, B2, and D2 consist of pure heavy component H. Therefore, temperature and composition estimates can be provided for the lower stages of C1 and all the stages of C2. The boiling point of n-heptane is supplied as the temperature estimates on these stages and we can assume that the liquid on these stages is pure n-heptane only. Then, a similar procedure for tracing the LSS branch is applied, starting from D1 = Ex, which is located outside the region of multiple steady states, D1 decreases continuously with small steps trying to reach the second turning point (TP2). To refine the value of TP2, smaller steps should be used. Then, the HSS branch in the bifurcation diagrams is traced. 2.2.2. Tracing the USS Branch. It is anything but easy to find out the unstable solutions. By simulations, we found that even if pretty good temperature and composition estimates have been

(3)

Then, all streams in the flowsheet can be obtained from the material balances: F1 = D1 + B2 = E x + D2 D1 D1 1 = = F1 D1 + B1 1+

(4)

( )

B1 D B 1 = 1· 1 = F1 F1 D1 1+

B1 D1

( )

B2 E x − D1 = F1 F1

B1 D1

·

B1 D1

(5)

(6)

(7)

Now the compositions and flow rates of all streams are specified for the ∞/∞ case and all the material balances are fulfilled. The overall material balance relationship is the basis of ∞/∞ predictions. 2.1. ∞/∞ Analysis. ∞/∞ analysis is very useful in determining the existence of multiple steady states for a column sequence,21 and the composition profiles predicted by this method can be used as the initial estimates for process simulations. In the direct split sequence, if the external feed flow rate Ex, composition, and the flow rate of recycle stream D2 are specified, only one degree of freedom is left. Then, B2 (or D1) can be easily determined with the specified D1 (or B2) according to the material balance equations. Thus, both D1 and B2 can be chosen as the bifurcation parameter. In this paper, D1 is chosen as the bifurcation parameter for the ∞/∞ bifurcation analysis. Then, the ∞/∞ bifurcation analysis of the direct split sequence is made and the composition bifurcation diagrams of D1, B2, B1 and D2 are constructed. The ∞/∞ bifurcation diagrams are shown in the Supporting Information (Figure S1). In the bifurcation diagrams, the continuation paths are divided into three branches: (1) D1 ∈ [0, Lp1], as D1 increases gradually, D1 and D2 are located at the L−H azeotrope, the fraction of heavy component H in B2 increases; (2) D1 ∈ [Lp1,Lp2], as D1 decreases gradually, the fraction of light component L in D1 increases to 1, and the fraction of heavy component H in B2 also increases. It is noteworthy that the compositions of the two desired product (D1 and B2) are unique when D1 is set to be Lp1 or Le, while the compositions of intermediate streams B1 and D2 can be varied, and so do the corresponding column temperature and composition profiles. This branch is unstable and there exists an unstable steady-state solution; (3) D1 ∈ [Lp2, Ex], as D1 increases gradually, the fraction of light component L in D1 decreases from 1 to xExL, while the fraction of heavy component H in B1 increases to 1 and then keeps constant. In this paper, considering the light component fraction in D1 (xD1L) within the multiplicity region D1 ∈ [Lp1,Lp2], the three branches mentioned above are called lower steady state (LSS), unstable steady state (USS), and higher steady state (HSS), respectively. 2.2. Strategies for Identifying Multiple Steady States Using Aspen Plus. In Güttinger and Morari’s work,21 Aspen Plus was used to detect the multiple solutions in the intermediate entrainer separation sequences. Their results confirmed that the steady-state multiplicities may exist in the column operating at

Figure 3. Case study for identifying multiple steady states using Aspen Plus. C


Article

Industrial & Engineering Chemistry Research

Figure 4. Simulated composition bifurcation diagrams for the two products D1 and B2 (in the figure, panels a, b, and c correspond to D1 and panels d, e, and f correspond to B2).

D1, and this might be the reason for the failure of the sensitivity analysis. In this paper, a novel strategy using the “Design Specs/ Vary” function is proposed to detect the unstable steady-state solutions, and the effectiveness of this strategy is proven by the simulation results. This four-step strategy lays out as follows: first, the temperature and composition profiles predicted by the

provided for the two columns, and sensitivity analysis is performed with small steps, the simulations may fail to converge. That is to say, sensitivity analysis is not always very effective in finding the unstable solutions. Basing on the ∞/∞ bifurcation diagrams, we find that process variables such as xD1H, xD1L, xB2L, and xB2H, etc. are extremely sensitive to the bifurcation parameter D


Article


Figure 5. Simulated composition bifurcation diagrams for the two internal streams B1 and D2 (in the figure, panels a, b, and c correspond to B1 and panels d, e, and f correspond to D2).

∞/∞ analysis is supplied to facilitate the convergence of the simulation; second, one of the process variables such as the fraction of heavy component H in B2 (xB2H) is assigned as the “Design Specs”, and the bifurcation parameter D1 is adjusted as the “Vary” to meet this specification; third, the specified value of the process variable is adjusted continuously with very small

steps, and the corresponding D1 is obtained as it was done in the second step. Note that, in all these cases, the temperature and composition profiles simulated in the last one are supplied to accelerate the convergence of the next; finally, the unstable steady-state solution can be traced, and the USS branch is plotted using the simulated results. E


Article


Figure 6. Temperature profiles for three steady states at D1 = 2551.02 kg/h.

Figure 7. Liquid composition profiles for three steady states at D1 = 2551.02 kg/h.

2.2.3. Case Study. The case study with detailed column operating parameters is shown in Figure 3. Both columns have 50 theoretical trays and operate under atmospheric pressure. We use the notation of numbering stages from the top, with stage 0 and 51 being the reflux drum and reboiler, respectively. The fresh feed is a 50/50 wt % of acetone/n-heptane mixture with a flow rate of 5000 kg/h and 8000 working hours per year. To balance tiny entrainer losses in these two product streams, a small entrainer makeup with flow rate of 55.2 kg/h is added to the system. The flow

rate of the recycle stream is specified at REC = 5000 kg/h, and the two product specifications are made: xD1L ≥ 98 wt %, xD1I ≤ 0.2 wt %, xB2H ≥ 98 wt %. Using the above-mentioned strategy, all the branches can be traced and the bifurcation diagrams are obtained; the simulated results of the D1, B2, B1, and D2 compositions along the continuation paths are shown in Figures 4 and 5. For a better representation of local continuation paths in the vicinity of two turning points, local magnified images are also placed on the corresponding bifurcation diagrams. F


Article

Industrial & Engineering Chemistry Research As seen in Figures 4 and 5, the simulated results reveal that three steady states exist when D1 ∈ [2533.0, 2848.2] kg/h. However, for this case study, the two turning points predicted by ∞/∞ analysis are located at D1 = 2844.7 kg/h and D1 = 2500.0 kg/h, respectively. This small deviation between the simulated and ∞/∞ analysis predicted results can be attributed to the differences between practical operating conditions and the ∞/∞ case. The rigorously simulated bifurcation diagrams are in excellent agreement with the ∞/∞ predicted results shown in Figure S1 (Supporting Information). A closer examination reveals that nonmonotonic curves are presented for all the investigated variables, and some segments of the continuation paths show curious behaviors. For example, in the vicinity of the second turning point (TP2) in Figure 4, there is a small segment on the HSS branch, in which the fraction of heavy component H in the distillate D1 decreases when D1 increases; that is, the reflux ratio decreases and the result is counterintuitive. To illustrate the steady-state multiplicities more precisely, the bifurcation parameter D1 is specified at 2551.02 kg/h, and the unstable solution corresponds to xD1L = 0.98 and xB2H = 0.98 at this specified value. The corresponding feed stages in the case study showed in Figure 3 are exactly the feed stages that minimize the total reboiler heat duty QRtot (QRtot = QR1 + QR2) at this specification. For the three steady states at D1 = 2551.02 kg/h, the temperature profiles and liquid composition profiles of the two columns are shown in Figures 6 and 7, respectively. For the acetone/benzene/n-heptane system, especially for the results locating on the USS branch, the composition profiles depicted in Figure 7 are in agreement with those reported in Kannan’s19 work to some extent. The consistency of the results verifies the effectiveness of the method proposed in this paper for multiple steady states tracing. For the LSS branch (shown in Figures 6 and 7), with the exception of sharp changes in the column condensers and reboilers, the temperatures for both columns rise slowly from the top to the bottom. The composition of D1 is very close to that of the acetone/n-heptane azeotrope. The C2 bottom product B2 is primarily composed of acetone and n-heptane. The liquid compositions on the upper stages for both columns are almost unchanged. Drastic changes in the liquid composition have been observed near the bottom of both columns; for the HSS branch, the temperatures change significantly for both columns. Two temperature fronts are separately located at the upper part and the bottom of C1, and there is a temperature front at the bottom of C2. The product streams D1 and B2 are mainly composed of acetone (98.00 wt %) and n-heptane (99.39 wt %), respectively, with benzene being the major impurity in both streams. The liquid compositions vary greatly in the vicinity of corresponding temperature fronts. Throughout the second column, an extremely low level of acetone is observed; for the USS branch, the two temperature fronts are both located at the lower part of C1, while the temperature profile for C2 is similar to that of HSS. D1 consists mainly of acetone (98.00 wt %), with n-heptane being the major impurity. B2 is composed primarily of n-heptane (98.00 wt %), with benzene being the major impurity. In the vicinity of temperature fronts, the liquid compositions for both columns also show remarkable changes. 2.3. Operating at Stable Solution or Unstable Solution Branch. After all the steady-state solution branches being traced, one may ask the following question: For given design specifications, on which branch is the operating point more reasonable? The case study in section 2.2.3 is carried out for the economical analysis of locating the operating point on different

Table 1. Economical Comparison for Operations on the USS and HSS Branches parameters

USS

HSS

ID1 (m) HC1 (m) QC1 (kW) QR1 (kW) AC1 (m2) AR1 (m2) vessel and tray costs (103 $) HX costs (103 $)

1.227 30 1593.8 1713.7 113.9 25.85 322.058 242.070

1.612 30 2867.5 2984.9 204.1 44.58 435.374 350.959

ID2 HC2 (m) QC2 (kW) QR2 (kW) AC2(m2) AR2 (m2) vessel and tray costs (103 $) HX costs (103 $)

1.228 30 1405.3 1388.9 41.58 26.46 322.481 162.955

1.920 30 3455.4 3443.8 102.37 66.46 528.646 294.483

283.812 10.866 1049.565 644.533

588.069 22.908 1609.462 1147.465

C1

C2

steam cost (103 $/year) cold water cost (103 $/year) FCI (103 $) TAC (103 $/year)

Figure 8. Detailed flowsheet information with the operating point located on the USS branch.

steady-state branches. The product specifications (xD1L ≥ 98 wt %, xD1I ≤ 0.2 wt %, xB2H ≥ 98 wt %) can be achieved by operating either at a stable or at an unstable steady state. As seen in Figure 10, the composition of D1 is very close to that of the acetone/n-heptane azeotrope when the operating point is located on the LSS branch, the product quality requirements cannot be met. Thus, only the economical reasonability of operating at the USS and HSS branch is going to be discussed. For the case study, total annual cost (TAC) is used as the objective function. TAC mainly consists of two parts: the costs related to fixed capital investments (FCI) and process variable costs (steam, cooling water, and electricity consumption, etc.). On the basis of the TAC analysis, a comprehensive comparison for operating on the USS and HSS branches is made, and the detailed results are listed in Table 1. For more information on the TAC estimation based on Douglas’ work28 and the corresponding flowsheet details for operating on the USS or HSS branch, please refer to the Supporting Information. G


Article


Figure 9. Temperature profiles for both columns.

Figure 10. Liquid composition profiles for both columns.

Figure 11. Basic control scheme.

branch can be cut down to 56.17% of that being operated on the HSS branch. The results of this case study remind us to be aware of the possibility of multiple steady states in process design. For the separation process existing multiple steady states, all the solutions should be found out using a proper strategy. Otherwise, eligible and even more economical design alternatives may be omitted, and the process will be inevitably operated at the expense of excessive energy and fixed capital cost.

According to Table 1, with the operating point located on the USS branch, the steam consumption (QR1 + QR2) and the cooling water consumption (QC1 + QC2) are only 48.26% and 47.43% of that being operated on the HSS branch, respectively. Lower utility consumption means smaller heat exchanging areas (ARn and ACn) and smaller column sizes (IDn), and lower fixed capital investments (FCI) consequently. Taking both FCI and utility costs into account, the TAC of operating on the USS H


Article


Figure 12. Dynamic responses of the basic control structure: 20% feed flow rate disturbances.

3. UNSTABLE OPERATING POINT STABILIZATION

steady-state design with the operating point located on the USS branch, the control strategy investigations are carried out with Aspen Dynamics. Figure 8 gives the detailed flowsheet information when the operating point is located on the USS branch. The corresponding column temperature and liquid composition profiles are shown in Figures 9 and 10, respectively. To convert the steady-state design to the dynamic model, the sizes of the equipment (column shells, reflux drums, and column bases, etc.) and the plumbing system such as valves’ pressure drops need to be specified. The “Tray Sizing” in the Aspen Plus is utilized to determine the column diameters and both of the column diameters in this case are rounded to 1.3 m with the consideration of standardized equipment production. With the guide of Luyben’s book,31 the commonly used heuristic for reflux drums and column bases sizing is to provide ∼5 min of liquid holdup when half full and all the control valves’ pressure drops are ∼3 bar with the valve half open at the design flow rate. After all the above-mentioned items being specified, the flowsheet is

As stated above, it is much more profitable to operate the separation on the USS branch. However, the existence of multiple steady states poses a great challenge on stabilizing the operating point on the USS branch. In the last two decades, the implications of steady-state multiplicity on the control system design have been systematically analyzed by several scholars.4,10,16,26,29,30 They pointed out that transitions between parallel steady states may occur for small manipulated variable disturbances. The steady states transition has also been experimentally verified by Mohl.13 According to Jacobsen,26 the separations corresponding to the unstable operating point cannot be achieved with manual operation. However, the unstable operating point might be stabilized by a robust feedback control of a product composition or tray temperature. Bearing all these implications in mind, several temperature control structures are proposed to stabilize the operating point on the USS branch. On the basis of the I


Article


Figure 13. Dynamic responses of the basic control structure: feed composition disturbances.

pressure checked and then it can be converted to the dynamic model. On the basis of this dynamic model, several control structures are proposed to stabilize the operating point on the USS branch, and the effectiveness of these control structures is tested by feed flow rate and composition disturbances. 3.1. Basic Control Structure. Temperature is widely used to provide inferential control of compositions in industrial practice because its measurement is fast, reliable, and inexpensive. Thus, the control strategy investigation in this paper is confined to the temperature control, and the traditional PI control-based strategy is overwhelmingly preferred for the simplicity. Since this intermediate entrainer separation process is similar to the extractive distillation process, a control structure that is similar to that of extractive distillation may perform well. The temperature control of the extractive distillation process has been studied by Luyben,32−34 Arifin,35 and Wang,36 and these control structures are the same in nature, but used for different systems.

Figure 14. Pressure effect on the liquid phase composition of benzene/ n-heptane mixture at 370.2 K. J


Article


Figure 15. Differential temperature control structure.

Figure 16. Dynamic responses of the differential temperature control structure: 20% feed flow rate disturbances. K


Article


Figure 17. Dynamic responses of the differential temperature control structure: feed composition disturbances.

by the flow rates of the heating steam into C1 and C2 reboilers, respectively (reverse acting). Most of the control loops discussed here are commonly used distillation control strategies. In this paper, all the flow controllers are PI with normal settings: Kc =0.5 and τI = 0.3 min. Maintaining a liquid level at set-point value is often not necessary, so all the level loops are P-only to provide maximum flow smoothing. A noteworthy feature revealed here is the inventory control suggested by Luyben,31−33 Arifin,34 and Wang,35 the reflux drum level in C2 is held by the makeup benzene flow rate. In the C2 reflux drum level loop, Kc =10 is used to simulate the very fast control for the entrainer makeup stream. Kc = 2 as suggested by Luyben30 is used in the other level loops. The pressure controllers are PI with default values: Kc = 20 and τI =12 min. The tray pressure profiles are dynamically calculated in the simulator. In considering small measurement and actuator lags, a 1 min dead-time element is inserted into each temperature control loop to offer a realistic simulation of dynamic performance. Taking the interaction of the two control loops into account, the sequential tuning method recommended by Luyben30 is used to

On the basis of the control strategy investigations of extractive distillation processes, the basic control structure is proposed. Figure 11 gives the basic control structure for this process and the SVD analysis is used to select temperature control stages. In this basic control structure, the temperature control stages are 38 and 46 for C1 and C2, respectively. The details of selecting temperature control stages can be referred to the Supporting Information. As shown in Figure 11, the details of the basic control scheme are described as follows: (1) Feed is flow-controlled (reverse acting); (2) base levels in C1 and C2 are held by manipulating the corresponding flow rate of bottom stream (direct acting); (3) reflux drum level in C1 is held by manipulating distillate flow rate (direct acting); (4) reflux drum level in C2 is held by manipulating the makeup benzene flow rate (reverse acting); (5) flow rate of stream REC is ratioed to the feed flow rate; (6) top pressures in both columns are controlled by manipulating the corresponding cooling water flow rate in the condenser (reverse acting); (7) reflux ratios in both columns are fixed; (8) temperatures on stage 38 in C1 and stage 46 in C2 are controlled L


Article


Figure 18. Pressure-compensated temperature control structure.

TC2 is fixed, when the feed flow rate increases, the pressure on the temperature control tray increases, remarkable decrease in the n-heptane concentration on this stage is expected due to its high sensitivity to pressure, thus considerable benzene escapes from the bottom, and large offset occurs, and the vice versa. In one word, due to the pressure effect, the basic control structure cannot handle feed flow rate disturbances well. Several techniques have been proposed for overcoming the adverse effect of pressure fluctuation problems: pressure-compensated temperature control, differential temperature control, and doubledifferential temperature control.37 In this paper, the differential temperature control and pressure-compensated temperature control will be discussed later. 3.2. Differential Temperature Control Structure. Differential temperature control has been widely used in control of multicomponent distillation processes. Kister38 introduced the differential temperature control in details on the basis of precedent researches. For the differential temperature control loop in a single column, the temperature Tn on the temperature control tray n is measured as usual, and the temperature Tm on a second tray m where the temperature is relatively insensitive to composition is also measured. Then, the temperature Tm is subtracted from the temperature Tn, giving a differential temperature ΔT. When the column pressure changes, both Tn and Tm change equally, and the differential temperature ΔT will be

tune the two temperature controllers TC1 and TC2. For the two temperature control loops, relay-feedback tests are run to achieve the ultimate gain Ku and period Pu, and then the Tyreus− Luyben tunings are used to obtain the gain Kc and integral time constant τI. The effectiveness of this basic control structure is tested by fresh feed flow rate and composition disturbances. For each case, a 20% step change is introduced to the fresh feed stream at t = 0.2 h, and the dynamic responses for the positive and negative disturbances are represented by the solid and dashed lines, respectively. Note that, the positive (negative) feed composition disturbance means the composition increase (decrease) of acetone in the feed. The dynamic responses of this basic control structure for the feed flow rate and composition disturbances are shown in Figures 12 and 13. As can be seen in Figure 12, it takes about 7 h to come to a new steady state for ±20% feed flow rate disturbances. The acetone product purities are held fairly close to their desired values at the new steady states, but this is not the case for the n-heptane product. There is large offset in the purity of n-heptane product; the n-heptane product purity drops to 96.52 wt % for the +20% increase in feed flow rate, and rises to 98.70 wt % for the −20% increase in feed flow rate. In Figure 13, the system can be brought back to a new steady state within 10 h, and both of the two product streams can be maintained closely at the corresponding specifications. These results indicate that this basic control structure cannot handle feed flow rate disturbances well. The reasons for the ineffectiveness of this basic control structure are discussed here. In this separation sequence, the separation implemented in C2 can be regarded as binary separation between benzene and n-heptane. Vapor−liquid calculations are carried out for the benzene−n-heptane binary system at the set point of TC2 (370.2 K), and the effect of pressure on the liquid phase composition of benzene−n-heptane binary system at 370.2 K is plotted in Figure 14. It reveals that the liquid phase composition of the benzene−n-heptane binary system is very sensitive to pressure when the system temperature is fixed. When the feed flow rate increases, the holdup in each stage increases, even if the top pressure of C2 is controlled pretty well, the pressures on all stages increase accordingly, the pressure changes on the stages closed to the column base is more striking. If the set point of the

Figure 19. Effect of system pressure on the bubble point of liquid mixture on stage 46 in C2. M


Article


Figure 20. Flowsheet equation for pressure-compensated temperature control.

Figure 21. Dynamic responses of the pressure-compensated temperature control structure: 20% feed flow rate disturbances.

selected. The details of selecting the temperature control trays can also be referred to the Supporting Information. For the case study in this article, the temperatures on tray 46 and tray 16 in C2 are selected as reference temperature and the second temperature, respectively. Then, the differential temperature ΔT = T2,46−T2,16 can be controlled by manipulating the flow rate of the heating steam. The differential control structure is presented in Figure 15. In this differential temperature control structure, the same tuning method is used and the differential temperature controller (DTC) is tuned with other controllers with the same tuning parameters as those in the basic control structure.

held constant consequently. Thus, using the different temperature ΔT as a controlled variable might be effective in overcoming the adverse effect of pressure fluctuation. As mentioned above, the poor dynamic responses of the basic control structure for feed flow rate disturbances are caused by the pressure effect in C2. Thus, the differential control loop is only implemented in C2 and the other control loops are the same with those in the basic control structure. For the differential control, selecting the second temperature is critical to the dynamic performance of this control scheme. According to the SVD analysis method proposed by Ling,39 the temperature control trays for this differential temperature control structure are N


Article


Figure 22. Dynamic responses of the pressure-compensated temperature control structure: feed composition disturbances.

temperature control, resulting in more heating steam being fed to the C2 reboiler; thus, the amount of benzene acting as an impurity in the C2 bottom stream will be remarkably reduced, and vice versa. 3.3. Pressure-Compensated Temperature Control Structure. A brief description of the pressure-compensated temperature control can be found in Buckley’s book.40 This control strategy has been introduced by Luyben33,41 and Yu42 to the control of heat-integrated pressure swing distillation, and robust control can be achieved. To overcome the adverse effect of pressure fluctuation, the pressure-compensated temperature control loop is implemented on stage 46 in C2 and the other control loops are the same with those in the basic control structure. This pressurecompensated temperature control scheme is shown in Figure 18. Before setting up the pressure-compensated temperature control structure, the relationship between the bubble point of the liquid phase on the stage 46 in C2 (0.00/72.61/27.39 wt % wt % acetone/n-heptane/benzene) and the pressure is required. For the liquid phase on stage 46 in C2, the bubble point temperatures within a specified pressure range (1.15−1.45 bar) are calculated

After the control structure is constructed, the effectiveness of this differential temperature control scheme is tested with the same disturbances posed in the basic control structure. The dynamic responses of this differential control structure for the feed flow rate and composition disturbances are presented in Figures 16 and 17, respectively. As seen in Figures 16 and 17, this differential temperature control structure can handle both feed flow rate and composition disturbances fairly well. The system takes about 10 h to come to a new steady state. Both products are maintained close to their specification levels at the new steady state. Compared with dynamic responses shown in Figure 12, the n-heptane product purity increases from 96.52 to 97.65 wt % for the +20% increase in feed flow rate, and drops from 98.70 to 98.22 wt % for the −20% increase in feed flow rate. The explanation for the improvement in dynamic performance lays out as follows: when the feed flow rate increases, the pressures on all stages and the temperature on stage 16 in C2 will increase accordingly; then, the temperature on stage 46 in C2 also increases rather than be held constant under this differential O


Article


Figure 23. Comparison of dynamic performances between differential temperature control and pressure-compensated temperature control: feed flow rate disturbances.

and the results are plotted in Figure 19. As seen in Figure 19, there is a linear relationship between the bubble point temperature and pressure, with a slope of 26.8966. Since the temperature controller PCTC is reverse-acting, the pressurecompensated temperature TPC can be calculated by the following equation:

TPC/K = T2,46/K − (P2R /bar − 1.3235) × 26.8966

where the P2R is the pressure of C2 reboiler and P2R is 1.3235 bar at the design conditions. In the pressure-compensated temperature control structure, the pressure-compensated temperature TPC is calculated by the P


Article


Figure 24. Comparison of dynamic performances between differential temperature control and pressure-compensated temperature control: feed composition disturbances.

calculation unit CAL, and TPC is fed to the pressure-compensated temperature controller (PCTC). TPC is held constant by manipulating the heating steam flow rate. In the Aspen Dynamics, the calculation of TPC is facilitated via the “Flowsheet Equations”, as shown in Figure 20. Note that the Aspen Dynamics numbers the stages from the top, with stage 1 being the reflux drum and the

last stage being the reboiler (stage 52 for a column with 50 stages), which is different from our notation. The effectiveness of this pressure-compensated temperature control structure is tested and the dynamic responses are shown in Figures 21 and 22. The results shown in Figures 21 and 22 reveal that both of the feed flow rate and composition disturbances can Q


Article


simulation. However, the unstable steady-state solution branch is still omitted. A novel four-step strategy combining the ∞/∞ analysis and the “design specs/vary” function is proposed in this article to trace the unstable steady-state solution branch, and the case study for separating acetone/n-heptane with benzene as intermediate entrainer confirms that this four-step strategy is effective. Qualitative and quantitative analysis is made for the case study. For the same product specifications, the steam consumption and TAC of operating on the unstable steady-state branch are only 48.26% and 56.17% of those operating on the stable steady-state branch. The results of the case study verify the conjecture that operating on the unstable steady state may be much more profitable, and remind us to be cautious to the occurrence of multiple steady states in process design. Because of the great advantages offered by operating on the unstable steady-state branch, several control structures are proposed to stabilize the operating point on this branch. The basic control structure proposed in this article provides an effective control for the feed composition disturbances. However, for the feed flow rate disturbances, large offsets occur in the n-heptane product composition as a result of pressure effect. To alleviate the adverse effect of pressure fluctuation, differential temperature control and pressure-compensated temperature control are proposed, and the effectiveness of these tow control structures are tested. The testing results reveal that both feed flow rate and composition disturbances can be fairly well handled by these two modified temperature control structures, and a somewhat similarity can be observed for these two alternative control schemes. Thus, both of the modified control schemes can be used for stabilizing the design, which provides a useful guide for the application of an unstable operating point in chemical industry.

be fairly well handled by this pressure-compensated temperature control structure. The purity of each product stream can be maintained closely at the desired value. Compared with the basic control structure, at the new steady states, the n-heptane product purity increases from 96.52 to 98.06 wt % for the +20% increase in feed flow rate, and drops from 98.70 to 97.98 wt % for the −20% increase in feed flow rate. 3.4. Comparisons between the Differential Temperature Control and Pressure-Compensated Temperature Control. As presented above, both differential temperature control and pressure-compensated temperature control can cure the problems caused by pressure fluctuation. In this section, a direct comparison of these two control structures is made. This direct comparison is presented in Figures 23 and 24, and the dynamic responses of differential temperature control and pressure-compensated temperature control are represented by the solid and dashed lines, respectively. As seen in Figure 23, for the feed flow rate disturbances, the acetone product purities at the corresponding new steady states are almost the same. However, for the differential temperature control structure, the steady-state offset in the purity of n-heptane product is larger than that of the pressure-compensated temperature control scheme. Thus, the pressure-compensated temperature control structure performs somewhat better than the differential temperature control structure for the feed flow rate disturbances. This tiny difference between these two control structures can be explained as follows: considering an increase of the feed flow rate here, all the stage pressures and temperatures in C2 increase accordingly; since the temperature on stage 16 in C2 is less sensitive to the heat input of C2 reboiler, the increase in T2,16 is smaller than that in T2,46; then, a smaller increase in the temperature on stage 46 is expected in the differential temperature control structure; thus, compared with the pressurecompensated temperature control, the C2 heat input is smaller, resulting in more benzene escaping from the C2 bottom, and vice versa. This explanation has been verified by the rigorous simulations. For the +20% feed flow rate disturbance, the temperature on stage 46 in C2 rises from 370.17 to 371.71 and 372.52 K for the differential temperature and pressurecompensated temperature control structures, respectively. However, the comparisons shown in Figure 24 reveal that the differential temperature control scheme performs somewhat better than the pressure-compensated temperature control structure. Smaller offsets occur under the differential control scheme for a similar reason stated above. When the feed is switched to 60/40 wt % of acetone/n-heptane, the feed to C2 is going to decrease. Then, a smaller decrease in T2,46 is expected under the differential temperature control structure. Therefore, a smaller amount of benzene escaping from C2 bottom is achieved, and vice versa. In brief, both differential temperature and pressure-compensated temperature control structures can handle the feed flow rate and composition disturbances fairly well, and the dynamic performances of these two alternatives are somewhat similar. Thus, both control schemes can be used for stabilizing the unstable operating point.

■

ASSOCIATED CONTENT

S Supporting Information *

(1)∞/∞ bifurcation diagrams for the acetone/benzene/ n-heptane mixture; (2) cost estimation; (3) selecting the temperature control stages for the basic control structure; (4) selecting the temperature control stages for the differential temperature control structure. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b00572.

■

AUTHOR INFORMATION

Corresponding Author

*Tel: +86 022-27892145. Fax: +86 022-27404440. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Doherty, M. F.; Perkins, J. D. On the Dynamics of Distillation Processes IV: Uniqueness and Stability of the Steady State in Homogeneous Continuous Distillations. Chem. Eng. Sci. 1982, 37, 381−392. (2) Petlyuk, F. B.; Avetyan, V. S. Investigation of Three Component Distillation at Infinite Reflux. Theor. Found. Chem. Eng. 1971, 5, 499− 506. (3) Magnussen, T.; Michelsen, M. L.; Fredenslund, A. Azeotropic Distillation Using UNIFAC. Inst. Chem. Eng. Symp. Ser. 1979, 56, 1−9. (4) Bekiaris, N.; Meski, G. A.; Radu, C. M.; Morari, M. Multiple Steady States in Homogenous Azeotropic Distillation. Ind. Eng. Chem. Res. 1993, 32, 2023−2038.

4. CONCLUSIONS In this article, the steady-state multiplicities of the intermediate entrainer separation sequence are thoroughly studied with the ∞/∞ analysis. With the initial estimates provided by the ∞/∞ analysis, the two stable steady-state solution branches can be easily traced by the sensitivity analysis in a rigorous Aspen R


Article

Industrial & Engineering Chemistry Research (5) Bekiaris, N.; Meski, G. A.; Morari, M. Multiple Steady States in Heterogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1996, 35, 207−227. (6) Güttinger, T. E.; Dorn, C.; Morari, M. Experimental Study of Multiple Steady States in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1997, 36, 794−802. (7) Müller, D.; Marquardt, W. Experimental Verification of Multiple Steady States in Heterogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1997, 36, 5410−5418. (8) Wang, C.; Wong, D.; Chien, I. L.; Shih, R. F.; Wang, S. J.; Tsai, C. S. Experimental Investigation of Multiple Steady States and Parametric Sensitivity in Azeotropic Distillation. Comput. Chem. Eng. 1997, 21, S535−S540. (9) Singh, B. P.; Singh, R.; Kumar, M. V. P.; Kaistha, N. Steady State Analysis of Reactive Distillation Using Homotopy Continuation. Chem. Eng. Res. Des. 2005, 83, 959−968. (10) Kumar, M. V. P.; Kaistha, N. Steady-state Multiplicity and Its Implications on the Control of an Ideal Reactive Distillation Column. Ind. Eng. Chem. Res. 2008, 47, 2778−2787. (11) Purohit, J. L.; Mahajani, S. M.; Patwardhan, S. C. Analysis of Steady-State Multiplicity in Reactive Distillation Columns. Ind. Eng. Chem. Res. 2013, 52, 5191−5206. (12) Chen, F.; Huss, R. S.; Doherty, M. F.; Malone, M. F. Multiple Steady States in Reactive Distillation: Kinetic Effects. Comput. Chem. Eng. 2002, 26, 81−93. (13) Mohl, K. D.; Kienle, A.; Gilles, E. D.; Rapmund, P.; Sundmacher, K.; Hoffmann, U. Steady-State Multiplicities in Reactive Distillation Columns for the Production of Fuel Ethers MTBE and TAME: Theoretical Analysis and Experimental Verification. Chem. Eng. Sci. 1999, 54, 1029−1043. (14) Delgado-Delgado, R.; Barroso-Muñoz, F. O.; Segovia-Hernández, J. G.; Hernández-Escoto, H.; Castro-Montoya, A. J.; Rico-Ramírez, V.; Hernández, S. Multiple Steady States in Thermally Coupled Distillation Sequences: Revisiting the Design, Energy Optimization, and Control. Ind. Eng. Chem. Res. 2014, 53, 17515−17521. (15) Gani, R.; Jørgensen, S. B. Multiplicity in Numerical Solution of Non-linear Models: Separation Processes. Comput. Chem. Eng. 1994, 18, S55−S61. (16) Bekiaris, N.; Morari, M. Multiple Steady States in Distillation: ∞/∞ Predictions, Extensions, and Implications for Design, Synthesis, and Simulation. Ind. Eng. Chem. Res. 1996, 35, 4264−4280. (17) Kienle, A.; Marquardt, W. Bifurcation Analysis and Steady-State Multiplicity of Multicomponent, Non-Equilibrium Distillation Processes. Chem. Eng. Sci. 1991, 46 (7), 1757−1769. (18) Ramzan, N.; Faheem, M.; Gani, R.; Witt, W. Multiple Steady States Detection in a Packed-Bed Reactive Distillation Column Using Bifurcation Analysis. Comput. Chem. Eng. 2010, 34, 460−466. (19) Kannan, A.; Joshi, M. R.; Reddy, G. R.; Shah, D. M. MultipleSteady-States Identification in Homogeneous Azeotropic Distillation Using a Process Simulator. Ind. Eng. Chem. Res. 2005, 44, 4386−4399. (20) Vadapalli, A.; Seader, J. D. A Generalized Framework for Computing Bifurcation Diagrams Using Process Simulation Programs. Comput. Chem. Eng. 2001, 25, 445−464. (21) Güttinger, T. E.; Morari, M. Multiple Steady States in Homogeneous Separation Sequences. Ind. Eng. Chem. Res. 1996, 35, 4597−4611. (22) Esbjerg, K.; Andersen, T. R.; Müller, D.; Marquardt, W.; Jørgensen, S. B. Multiple Steady States in Heterogeneous Azeotropic Distillation Sequences. Ind. Eng. Chem. Res. 1998, 37, 4434−4452. (23) Matauyama, H.; Nishimura, H. Topological and Thermodynamic Classification of Ternary Vapor-Liquid Equilibria. J. Chem. Eng. Jpn. 1977, 10, 181−187. (24) Dorn, C.; Güttinger, 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−515. (25) Baharev, A.; Neumaier, A. A Globally Convergent Method for Finding All Steady-State Solutions of Distillation Columns. AIChE J. 2014, 60, 410−414.

(26) Jacobsen, E. W.; Skogestad, S. Multiple Steady States and Instability in Distillation: Implication for Operation and Control. Ind. Eng. Chem. Res. 1995, 34, 4395−4405. (27) Chien, I. L.; Wang, C. J.; Wong, D. S. H. Dynamics and Control of a Heterogeneous Azeotropic Distillation Column: Conventional Control Approach. Ind. Eng. Chem. Res. 1999, 38, 468−478. (28) Douglas, J. M. Conceptual Design of Chemical Processes; McGraw Hill: New York, 1998; pp 568−577. (29) Sneesby, M. G.; Tadé, M. O.; Smith, T. N. Steady-state Transitions in the Reactive Distillation of MTBE. Comput. Chem. Eng. 1997, 22, 879−892. (30) Kumar, M. V. P.; Kaistha, N. Role of Multiplicity in Reactive Distillation Control System Design. J. Process Control 2008, 18, 692− 706. (31) Luyben, W. L. Distillation Design and Control Using Aspen Simulation; Wiley: New York, 2006. (32) Luyben, W. L. Plantwide Control of an Isopropyl Alcohol Dehydration Process. AIChE J. 2006, 52, 2290−2296. (33) Luyben, W. L. Comparison of Extractive Distillation and Pressure-Swing Distillation for Acetone−Methanol Separation. Ind. Eng. Chem. Res. 2008, 47, 2696−2707. (34) Luyben, W. L. Effect of Solvent on Controllability in Extractive Distillation. Ind. Eng. Chem. Res. 2008, 47, 4425−4439. (35) Arifin, S.; Chien, I. L. Design and Control of an Isopropyl Alcohol Dehydration Process via Extractive Distillation using Dimethyl Sulfoxide as an Entrainer. Ind. Eng. Chem. Res. 2008, 47, 790−803. (36) Wang, Q.; Yu, B.; Xu, C. Design and Control of Distillation System for Methylal/Methanol Separation. Part 1: Extractive Distillation using DMF as an Entrainer. Ind. Eng. Chem. Res. 2012, 51, 1281−1292. (37) Yu, C. C.; Luyben, W. L. Use of Multiple Temperatures for the Control of Multicomponent Distillation Columns. Ind. Eng. Chem. Process Des. Dev. 1984, 23, 590−597. (38) Kister, H. Z. Distillation Operations; McGraw-Hill, Inc.: New York, 1990; pp 563−566. (39) Ling, H.; Luyben, W. L. Temperature Control of the BTX Divided-wall Column. Ind. Eng. Chem. Res. 2010, 49, 189−203. (40) Buckley, P. S.; Luyben, W. L.; Shunta, J. P. Design of Distillation Column Control Systems; Instrument society of America: Research Triangle Park, NC, 1985; pp 234−239. (41) Luyben, W. L. Design and Control of a Fully Heat-Integrated Pressure-Swing Azeotropic Distillation System. Ind. Eng. Chem. Res. 2008, 47, 2681−2695. (42) Yu, B.; Wang, Q.; Xu, C. Design and Control of Distillation System for Methylal/Methanol Separation. Part 2: Pressure Swing Distillation with Full Heat Integration. Ind. Eng. Chem. Res. 2012, 51, 1293−1310.

S


Multiple Steady-States Analysis and Unstable ... - ACS Publications

Recommend Documents