Ind. Eng. Chem. Res. 2009, 48, 3873–3879
3873
PROCESS DESIGN AND CONTROL Time-Optimal Control of Startup Distillation of Nonideal Mixtures Alexandru Woinaroschy* Department of Chemical Engineering, “Politehnica” UniVersity of Bucharest, Polizu Str. 1-5, Bucharest, Romania
In this work, the time-optimal control of startup distillation columns by iterative programming proposed by Woinaroschy [Ind. Eng. Chem. Res. 2008, 47, 4158.] is extended to the case of nonideal mixtures. The minimization of distillation startup time is performed by iterative dynamic programming employing randomly chosen candidates for admissible control. The control variables are the reflux ratio and/or the reboiler heat duty. The dynamic distillation model proposed by the author in the mentioned work is applied with a preliminary computation of suitable correlation relations for equilibrium constants. An illustrative application for separation of a propene-propane mixture at plant scale is presented. Because of the small relative volatility of this mixture and respective high reflux ratio, the decrease of the reboiler heat consumption corresponding to the minimization of the startup time is significant. Introduction The ability of iterative programming (IDP) to solve high dimensional time-optimal control (TOC) problems involving complex models, such as startup distillation columns, was demonstrated in a previous article.1 Optimal control policies for minimization of startup time of distillation columns were obtained on the basis of a dynamic distillation model (DDM) representing a good compromise between correctness and complexity. This result was possible mainly because of the transformation of the differential algebraic equations model into an ordinary differential equations model, which is computationally much more efficient. The proposed DDM was applied for ideal and azeotropic mixtures, with expressions of equilibrium constants that can be derived analytically with respect to temperature and composition. For nonideal mixtures (distillation at low temperatures and/or high pressures), the calculation of equilibrium constants is made with implicit methods (e.g., Soave, Peng-Robinson, etc.). In this case, the application of numerical differentiation involved in the plate temperature equation is expected to fail. A more promising way to extend the use of the proposed DDM for this case is to obtain previous correlation relations for equilibrium constants and to derive them analytically. The extension to the case of nonideal mixture is presented in connection with an important industrial application: the separation of the propene-propane mixture.
consideration was given by HYSYS authors to these components in developing the binary interaction coefficients for the Peng, Robinson and Soave, Redlich, and Kwong equations of state to ensure that these methods correctly model the system. In the present article, the equilibrium and thermodynamic data were calculated on the basis of the Soave, Redlich, and Kwong equations of state. This separation needs many stages, and therefore two separate columns are considered: a stripper and a rectifier. The stripper is modeled in HYSYS as a reboiled absorber and contains 89 theoretical stages. The rectifier is a refluxed absorber containing 61 theoretical stages. The stripper contains two feed streams: one is the external feed (on stage 10 of the stripper), and the other is the bottom from the rectifier. Propane is recovered from the stripper bottoms, and propene is taken off from the top of the rectifier. The HYSYS flowsheet of propene-propane separation is presented in Figure 1, and the main data of simulation are given in Table 1.
Preliminary Steady-State Simulation of Propene-Propane Separation The separation of the propene-propane mixture is a tutorial application of the HYSYS simulator,2 which was used in the present work with slight modifications. The modifications of the original HYSYS case study consist in the total number of stages, position of feed tray, the value of reflux ratio, and constant pressure inside the column. The steady-state simulation of propene-propane separation with HYSYS is generally an easy column to converge. The critical factor in producing good results is not the ease of solution, but the accurate prediction of the relative volatility of the two components. Special * To whom correspondence should be addressed. Tel.: +40-214023902. Fax: +40-21-3185900. E-mail:
[email protected].
Figure 1. HYSYS flowsheet of propene-propane separation.
10.1021/ie801490r CCC: $40.75 2009 American Chemical Society Published on Web 03/20/2009
3874 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 Table 1. General Data data
value
number of plates feed plate position feed flow rate (kmol s-1) feed specification feed composition (mole fractions) feed temperature (K) column pressure (atm) type of condenser reflux ratioa type of reboiler reboiler heat dutya (kW)
150 (89 rectifier + 61 stripper) 99 (10 for stripper) 0.165 saturated vapor propene 0.6, propane 0.4 325.74 20 total 18 total 20153
Additional Data for DDM plate type plate diameter (m) plate area (m2) weir length (m) weir height (m) a
sieve with lateral downcomers 1.5 1.264 1.2 0.05
Prescribed values for the desired stationary regime.
Figure 2. Comparison between liquid concentrations (s DDM; - - HYSYS).
Startup Distillation Operational Procedure Because of the vapor state of the feed, the startup procedure formulated by Ruiz et al.3was slightly modified as follows: Step 1: The vapor feed is introduced on the feed plate and reaches the top of the column. Step 2: The condenser starts to operate as the vapor phase reaches the top plate and the reflux drum starts to fill up. Step 3: The reflux is introduced into the column and operation at total reflux begins. Step 4: The liquid starts to weep to the plate below through the plate holes rather than through the downcomer. The liquid reaches the bottom of the column (or reboiler), thus increasing the liquid level there. Step 5: Heat is introduced in the reboiler and vapors start to go up. Step 6: The vapor flow through the plate holes seals the plates in terms of liquid weeping through the plate holes. This starts to increase the liquid holdup on the plates. Step 7: All plates have enough liquid holdup so that the liquid can start to fall down the downcomers. The downcomers are sealed and vapors cannot pass upward through them. Step 8: Column operation is changed from the total reflux. Distillate, bottom product, and the possible side streams are taken out. In the frame of the present work, at the end of step 8 begins the effective startup transient operating regime procedure, which will be the object of the simulation. For the initial state, the liquid composition is considered the same on all plates, being in thermodynamic equilibrium with the vapor feed composition. Traditionally, the column is operated at constant values of control parameters (feed rate, thermal state of the feed, reflux ratio, reboiler heat duty, bottom pressure, etc.), and they have during the startup transient regime the same values with those corresponding to the final desired stationary regime. All values
of state variables (concentrations, temperatures, pressures on plates and in the output material streams) vary in time. The amplitude of these fluctuations decreases in time and will vanish at steady-state regime. Up to the stationary regime the compositions of the products are, of course, different from the desired values, the products being collected, mixed, and recycled to the feed. DDM Particularization for Propene-Propane Separation The following simplification assumptions are present in this model: (i) The molar vapor holdup is negligible compared to the molar liquid holdup. (ii) Interphase heat transfer is considered much more intense than interphase mass transport; consequently, the liquid and vapor leaving each plate are in thermal equilibrium at the boiling temperature, corresponding to the liquid composition. (iii) On each plate the liquid and vapor compositions correspond to the thermodynamic equilibrium (theoretical plates). (iv) Entrainment and weeping rates, flooding of the plates, downcomer holdup, and delay time between plates are neglected. (v) The variation of total pressure in the column is around 5%, and therefore the pressure is considered constant. For a general tray j, neglecting the molar vapor holdup, massbalance equations are thus defined: Total mass balance around plate j: dNj ) Lj-1 + Vj+1 - Lj - Vj ( FL,j ( FV,j dt
(1)
Mass balance around plate j for component i:
Table 2. Selection of Liquid Propene Concentration Profile at Final Desired State plate j x1,j plate j x1,j plate j x1,j
1 0.9753 55 0.7853 110 0.5373
5 0.9668 60 0.7617 115 0.4861
10 0.9549 65 0.7382 120 0.4185
15 0.9416 70 0.7151 125 0.3366
20 0.9268 75 0.6926 130 0.2479
25 0.9105 80 0.6715 135 0.1635
30 0.8927 85 0.6513 140 0.0930
35 0.8736 90 0.6331 145 0.0408
40 0.8530 95 0.6162 150 0.0054
45 0.8313 100 0.5997
50 0.8086 105 0.5741
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3875
dxi,j 1 ) [Lj-1(xi,j-1 - xi,j) + Vj+1(yi,j+1 - xi,j) dt Nj Vj(yi,j - xi,j) ( FL,j(xF,i,j - xi,j) ( FV,j(yF,i,j - xi,j)] (2) The vapor composition, according to assumption (iii), is yi,j ) κi,jxi,j
(3)
The equilibrium data for the mixture propene-propane at a constant column pressure of 20 atm were obtained by regression of HYSYS application using the Soave, Redlich, and Kwong equations of state: -6
-3
κ1,j ) 1.9636583 × 10 Tj - 2.5974968 × 10 Tj + 1.2884311Tj2 - 284.01531Tj + 23474.885 (4) 4
3
κ2,j ) 4.2997564 × 10-6Tj4 - 5.6805395 × 10-3Tj3 + 2.8144010Tj2 - 619.74136Tj + 51177.216 (5)
application: hj ) (-1.5368181 × 108 + 5.0848 × 105Tj)x1,j + (-1.2321179 × 108 + 4.0872 × 105Tj)x2,j (12) Hj ) (-1.2869932 × 107 + 1.0913 × 105Tj)y1,j + (-1.3230183 × 107 + 1.1177 × 105Tj)y2,j (13) and dhj ) 5.0848 × 105x1,j + 4.0872 × 105x2,j dTj
The original feature of the model1(because iterative solving of nonlinear algebraic equations is avoided) consists in the calculation of the temperature on each plate with the equation m
For eq 4, the average residual is 1.77 × 10-11 and the standard error is 4.49 × 10-5, and for eq 5, the average residual is -1.91 × 10-11 and the standard error is 2.83 × 10-4. The liquid flow rate is obtained on the basis of Francis’ correlation for a plate weir, and the resulting equation1 is
Lj )
1.84lwFL,j m
∑Mx
i i,j
(
m
Nj
∑Mx
i i,j
i)1
εL,jApFL,j
- zw
)
1.5
(6)
i)1
The liquid density in the above equation is computed using the regression of data from HYSYS application: FL,j ) (2.2126 + 0.0315Tj)x1,j + (2.0837 + 0.284Tj)x2,j (7) Vapor flow rate is obtained from total energy balance: d(Njhj) ) Lj-1hj-1 + Vj+1Hj+1 - Ljhj dt VjHj ( FL,jhF,j ( FV,jHF,j - qj (8) where the left side is developed to d(Nj xi,j) dxi,j dNj ) Nj + xi,j dt dt dt
(14)
dTj )dt
∑κ
i,j
i)1 m
dxi,j dt
dκi,j xi,j dTj i)1
∑
(15)
where dxi,j/dt is computed from eq 2, and dκi,j/dTi,j is obtained by derivation of eqs 4 and 5. In this way, knowing the liquid composition and temperature on each tray, we can obtain the vapor composition by straightforward computation, without iterative calculations. Equation 15 corresponds to eq 22 from the mentioned article1 with activity coefficients γi,j ) 1. The core of the model consists of n equations 1, n · (m-1) equations 2, and n equations 15, and therefore n · (m+1) differential equations. To avoid iterative calculations, the vapor flow rates and compositions are computed in the order j ) n, n - 1, ... 1, wheras the liquid flow rates are calculated in the order j ) 1, 2, ... n. The DDM was used for a dynamic simulation at traditional startup procedure (constant control parameters were R ) 18 and qB ) 20153 kW) using the same data as for HYSYS steadystate simulation and several additional data needed by DDM (Table 1). A comparison of the final, stationary liquid plate concentrations given by DDM and HYSYS corresponding
(9)
and from eqs 8 and 9 Vj )
(
1 L h + Vj+1Hj+1 - Ljhj ( FL,jhF,j ( Hj j-1 j-1 dhj dTj dNj - hj (10) FV,jHF,j - qj - Nj dTj dt dt
)
After replacing dNj/dt from eq 1 into eq 10, the vapor flow rate is Vj )
[
1 L (h - hj) + Vj+1(Hj+1 - hj) ( Hj j-1 j-1 FL,j(hF,j - hj) ( FV,j(HF,j - hj) - qj - Nj
]
dhj dTj (11) dTj dt
The liquid and vapor enthalpies in the above equation are computed using the regression of data from the HYSYS
Figure 3. Space-time propene liquid mole fraction evolution during the startup period.
3876 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009
Figure 6. Optimal reboiler heat duty control.
Figure 4. Optimal reflux control (S ) 5).
Figure 7. Optimal reflux ( · · · ) and reboiler heat duty (s) control.
To minimize the final time tf subject to the constraints in eq 16, the penalty performance index was formulated as
Figure 5. Optimal reflux control (S ) 10).
results is presented in Figure 2. The good agreement between these results can be noticed.
∑ I ) tf + ω
Numerical Procedure The TOC problem is to minimize the final time tf and to determine the corresponding control variables ui (τk), i ) 1, 2, ..., nu, k ) 1, 2, ..., S, and length of time stages V(k), k ) 1, 2, ..., S restricted to si(tf) ) ssi
ns
i ) 1, 2, ..., ns
(16)
where the desired final stationary state values ssiare computed from the differential equations of the continuous dynamic system dsi ) V(k)Sfi(s1, s2, ..., sns, u1, u2, ..., unu) dτ
i ) 1, 2, ..., ns; k ) 1, 2, ..., S (17)
Practically, eq 17 corresponds with eqs 1, 2, and 15 after introducing normalized time variable τ, such that τ ) t/tf, with discrete values τk ) k/S, k ) 0, 1, ..., S. In the transformed time domain, each stage is of equal length 1/S.
i)1
|
1ns
si(tf) ssi
|
(18)
where ω > 0 is a penalty coefficient. The application of TOC is the same as in the first article,1 using the algorithm proposed by Bojkov and Luus4 based on the IDP procedure given by Bojkov and Luus.5,6 The system of differential equations was numerically integrated by the fourth-order Runge-Kutta-Gill method. Because of the strong nonlinearity of eq 15, a small initial integration step must be used, 0.5 s for the first 10 min of time integration, and after that the integration step was increased to 3 s. In fact, the reliability of the entire numerical procedure depends only on the value of this integration step. These integration steps were established empirically, after several attempts. For better reliability and perhaps computational efficiency also, an integration procedure that automatically chooses the time step to satisfy local error tolerance can be used (e.g., the IMSL subroutine DVERK7). As in the first article, it can be appreciated that the corresponding increase of the computer time due to the small values of the integration step is justified by avoiding the iterative
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3877
Figure 10. Evolution of propene mole fraction in distillate. Figure 8. Suboptimal reflux control.
Figure 9. Evolution of propene mole fraction in distillate.
solution of algebraic equations, which would require more computational effort. The computer programs were coded in FORTRAN. Results and Discussion The column design and operating parameters are indicated in Table 1. The initial propene liquid mole fraction is 0.5739 on all plates. A selection of liquid propene concentration profiles at final desired state is indicated in Table 2. For traditional procedure at constant control parameters (R ) 18 and qB ) 20153 kW), the startup time is 560 min and energy consumption is 188 094 kWh. Space-time propene liquid mole fraction
evolution during the startup period for the traditional procedure is presented in Figure 3. Figures 4 and 5 present the optimal reflux control for 5 and 10 time stages. The solution for 10stage optimal reflux control is not a refined result of 5-stage optimal reflux control, which is quite different. The explanation of this fact consists in the multiplicity of optimal control solutions and in the use of the IDP procedure which employs randomly chosen candidates for the admissible control. Figure 6 presents the optimal reboiler heat duty control, and Figure 7 presents the combination of optimal reflux and reboiler heat duty control (all for 5 time stages). Figures 5 and 6 present, at the end of the respective optimal startup controls, a very small part of the stationary regime. The attempts to establish piecewise linear control policies failed, since the corresponding optimal startup periods are significantly greater than the above piecewise constant control policies. A possible explanation can be given by the use (according to the procedure proposed by Luus8) of a single point for the s-grid, in the frame of this very high dimensional application. Because, sometimes, a continuous control policy is preferred rather than a policy that requires sudden switching from one level to another, a suboptimal reflux control was proposed: for the 5-stage optimal reflux control a linear transition from one level to the next was imposed. The durations of these transitions were arbitrarily imposed between 5 and 20 min, proportionally with the magnitude of transition. This suboptimal policy is presented in Figure 8 and gives good results as indicated in Table 3. Figures 9 and 10 present the evolutions of propene mole fraction in distillate. Table 3 gives the performances of the optimal and suboptimal policies. The greatest reduction of the startup time is given by 10-stage optimal reflux control (31.15%), and the most important energy saving corresponds to the combination of optimal reflux and reboiler heat duty
Table 3. Performances of the Optimal and Suboptimal Control Policies control policy
startup time (min)
decrease of startup time (%)
reboiler energy saving (kW)
reboiler energy saving (%)
propene supplementary production (t)
propane supplementary production (t)
5-Ra 10-Ra qBa R&qBa s-Ra
388 384 406 418 427
30.77 31.15 27.45 25.28 23.75
57772 59115 63532 64190 52495
30.71 31.42 33.72 34.12 27.91
36.54 37.39 32.71 30.17 29.53
36.64 37.50 32.81 30.25 29.61
a
5-R: 5-stage optimal reflux; 10-R: 10-stage optimal reflux; qB: optimal reboiler heat duty; R&qB: combination of optimal reflux and reboiler heat duty; s-R: suboptimal reflux.
3878 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 Table 4. Computing Data control policy
IDP iterations number to obtain optimal solution
total number of integrations of differential equation system for 20 IDP iterations
computer execution time (h)
5-Ra 10-Ra qBa R&qBa
17 13 17 7
89605 381710 89605 138005
11.68 21.91 13.44 15.03
a 5-R: 5-stage optimal reflux; 10-R: 10-stage optimal reflux; qB: optimal reboiler heat duty; R&qB: combination of optimal reflux and reboiler heat duty.
Figure 12. Convergence profiles of performance index for optimal reflux and reboiler heat duty control.
Figure 11. Convergence profiles of performance index for optimal reflux control and optimal reboiler heat duty control.
control (64 190 kWh, respectively 34.12%). Because of the small relative volatility of this mixture, respective high reflux ratio, the decrease of the reboiler heat consumption corresponding to the minimization of the startup time is significant. The solution for 10-stage optimal reflux control is a little better than that corresponding to 5-stage optimal reflux control, but the computer effort is almost twice as big (Table 4). Table 3 also indicates the supplementary productions of propene and propane obtained because of the decrease in startup time. The economic significance of these results can be considered appreciable. The number of grid points for optimal reflux control and for optimal reboiler heat duty control was 5 for s-grid, 9 for u-grid, and 9 for V-grid. For combined optimal reflux and reboiler heat duty control, because of the increased dimensionality, the number of grid points was 5 for s-grid, 5 for u-grid, and 5 for V-grid. These decreasing grid points give a higher startup time (Table 3) than the use of each of the controls individually. In all applications, according to Bojkov and Luus4 recommendations, the region contraction factor was set at 0.8. The total number of IDP iterations was 20. The convergence profiles of performance index are presented in Figures 11 and 12, and the convergence profiles of final time are given in Figure 13. The applications were executed on an IBM Intel Pentium IV computer (2.66 GHz processor, 1 GB memory, and Microsoft Windows XP Professional SP2 operating system). The computing data are given in Table 4. The computer effort and, consequently, computer execution time are appreciable: for example, in the case of the 10-stage optimal reflux control
Figure 13. Convergence profiles of final time.
171 796 500 differential equations integrations are involved (381 710 integrations multiplied with 450 differential equations). Conclusions In this work, the ability of IDP to solve high dimensional TOC problems, involving complex models, was demonstrated once more. The dynamic distillation model proposed and applied to ideal and azeotropic mixtures by the author in the previous work1 was extended here to nonideal mixtures. This extension is based on a preliminary computation of suitable correlation relations for equilibrium constants, liquid and vapor enthalpies, and liquid density using a HYSYS simulator. A proposed suboptimal policy, which avoids sudden switching from one level to another, gives good results. For the separation of the propene-propane mixture at plant scale, an important industrial application, the established optimal control policies have an appreciable economical significance. Acknowledgment The research presented here had no kind of financial support. This article is dedicated to my professor, collaborator, and friend, Prof. Dr. Octavian Smigelschi.
Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3879 ω ) penalty coefficient
Nomenclature -1
F ) feed or output side streamflow rate (kmol s ) H ) vapor enthalpy (J kmol-1) h ) liquid enthalpy (J kmol-1) I ) performance index L ) liquid flow rate (kmol s-1) lw ) weir length (m) m ) number of components N ) liquid holdup per plate (kmol) n ) number of plates ns ) dimension of state vector ν ) dimension of control vector q ) heat rate losses per plate (W) qB ) reboiler heat duty (W) R ) reflux ratio S ) number of time stages s ) state vector T ) temperature (K) t ) time (s) u ) control vector V ) vapor flow rate (kmol s-1) V ) time stage (s) x ) liquid mole fraction on plate y ) vapor mole fraction on plate zw ) weir height (m) Greek Letters γ ) activity coefficient ε ) volumetric fraction κ ) equilibrium constant F ) density (kg m-3) τ ) normalized time variable
Subscripts 1 ) propene 2 ) propane F ) feed f ) final i ) component j ) plate k ) time stage L ) liquid s ) stationary V ) vapor
Literature Cited (1) Woinaroschy, A. Time-optimal control of startup distillation columns by iterative dynamic programming. Ind. Eng. Chem. Res. 2008, 47, 4158. (2) HYSYS 3.1 Tutorials & Applications; Hyprotech: Calgary, Alberta, Canada, 2002. (3) Ruiz, C. A.; Cameron, I. T.; Gani, R. A. Generalized model for distillation columns III. Comput. Chem. Eng. 1988, 12, 1. (4) Bojkov, B.; Luus, R. Time-optimal control by iterative dynamic programming. Ind. Eng. Chem. Res. 1994, 33, 1486. (5) Bojkov, B.; Luus, R. Use of random admissible values for control in iterative dynamic programming. Ind. Eng. Chem. Res. 1992, 31, 1308. (6) Bojkov, B.; Luus, R. Evaluation of the parameters used in iterative dynamic programming. Can. J. Chem. Eng. 1993, 71, 451. (7) Luus, R. IteratiVe Dynamic Programming; Chapman & Hall/CRC: London, 2000. (8) Luus, R. Piecewise linear continuous control by iterative dynamic programming. Ind. Eng. Chem. Res. 1993, 32, 859.
ReceiVed for reView October 2, 2008 ReVised manuscript receiVed February 18, 2009 Accepted February 23, 2009 IE801490R
