Ind. Eng. Chem. Fund8m. 1982,27, 323-325
323
COMMUNICATIONS Comparison of Energy Models for Dlstlllation Columns Three different forms of the energy equation for a distillation column tray are quantitatively compared. The most rigorous version is shown to give the most realistic responses, particularly for transients in vapor rates throughout the column. The most rigorous model increases computing time only very slightly.
Introduction Many simplifications are made in deriving dynamic models of distillation columns. One of the most common is to use an algebraic form of the energy equation. This eliminates one differential equation per tray. The justification for this simplifying assumption is that the rate of energy transfer is usually much faster than the rate of change of liquid holdup or composition in a given tray. Several authors (McCune and Gallier, 1973; Doukas and Luyben, 1978) have reported difficulties when this approach was taken. Abnormal dynamic responses in flow rates and compoeitions were observed. The purpose of this paper is to provide a quantitative comparison among three energy equation models. The first model completely neglects the energy derivative. The second model partially considers the energy derivative. The third model uses the complete energy derivative. The binary methanol-water system is considered. The open loop dynamic behavior given by each of these models is shown for step disturbances in reflux flow rate, feed flow rate, and feed composition. Dynamic Model Several assumptions were made in deriving the mathematical model. Pressure was assumed constant. The dynamics of the total condenser and partial reboiler were neglected. Feed and reflux were taken as saturated liquid. The trays were assumed to be 100% efficient. Vapor holdup was assumed to be negligible. The vapor-liquid equilibrium was given by simple polynomial relationships (See Appendix A).
Distillate and bottoms flow rates were manipulated by level controllers. D = D(MD) (10)
B = B(MB) (11) There were also energy and mass equations for the reboiler and condenser. The algorithm used at each point in time to solve the system of equations was as follows. Step 1. Calculate Y N and T N (eq 1 and 2). step 2. Calculate p ~ hN, , and HN (eq 3 and 5). Step 3. Calculate LN (eq 6). Step 4. Calculate VN (eq 9). Step 5. Evaluate derivatives (eq 7 and 8). Step 6. Integrate. Step 7. Go back to step 1. The only difference among the three models is how VN is evaluated. For model 1,the energy derivative is made equal to zero and VNcan be obtained from an algebraic equation (12) VN = [LN+lhN+l+ VN-1HN-1 -LNhNl/HN For model 2, the energy derivative is assumed to be given by WN~N " dt h N d t Using the expression for dMN/dt given by eq 7, VNcan be found from
VN = [LN+i(hN+i- h ~ +) VN-I(HN-I - ~ N ) ] / ( H-Nh ~ ) (13)
For model 3, the complete energy derivative is used. The same type of relationships were used for liquid density and liquid and vapor enthalpies (See Appendix B). P N = P(XN,TN) (3) (4) hN = h(xN,TN) HN = H(YN,TN) (5) The liquid flow from each plate was given by the Francis weir formula (Van Winkle, 1967) or in a simplified form LN = UMN) (6) The following differential equations account for mass and energy transfer in each tray. dMN/dt = LN+1 + vN-1- LN - V N (7)
WNXN =L /~ N +~~ x N++VN-IYN-I ~ - L N ~ -NVNYN (8) m&N/dt = LN+ihN+i+ VN-~HN-I - L N ~ -NVNHN(9) 0198-4313/82/1021-0323$01.25/0
W N ~-- N W --
N + MN-dhN (14) dt h N T dt dMN/dt can be found from eq 7, and dhN/dt, following Doukas' suggestion (Doukas and Luyben, 1978),was given by the expression hN(t f At) - hN(t) -dhN = (15) dt At V N is then given by VN = [ L ~ + i ( h ~ -+ h i ~+)VN-~WN-I - h ~- ) M ~ ( d h n / d t ) ] / [ H-~h ~ (16) ]
Open Loop Responses The system considered in this study was a 32-tray column separating a binary mixture of methanol and water. The system exhibits modest nonequimolal overflow be@ 1982 American Chemical Society
324
Ind.
Eng. Chem.
Fundam., Vol. 21, No. 3, 1982
Table 1. Column Steady-StateDesign Values feed rate, kg-mol/min 27.24 feed composition, 0.50 mole fraction methanol distillate rate, kg-mol/min 13.62 distillate composition, 0.999 mole fraction methanol bottoms rate, kg-molimin 13.62 bottoms composition, 0.001 mole fraction methanol operating pressure, atm 1 relative volatility 7.5-2.4 number of trays 32 feed tray location 6 tray spacing, m 0.6 tray efficiency, % 100 column diameter, m 3 reflux ratio 0.974 outlet weir height, cm 3.81 reboiler heat duty, l o 6 kcal/h 13.5
"6
25
....
MODEL I MODEL 2 MODEL 3
lX.
o . 5 5 k x6
045
v 3I 26.5
0 99651
P76j 1
0
30
60
90
TIME (mlnut*S)
394r-7--
Figure 2. Openloop response for a 10% increase in feed flow rate. :,-.. )Io.
.,.,..
' 0°*[ /------x,,
MODEL I 23
- M OO DD EE LL
L 6 390
L 0 9951
0 9 9c
L-__--.--C
30 TIME
60
90
2 5.4
I
(minutes1
Figure 1. Openloop response for a 10% increase in reflux rate.
havior. Numerical values of design parameters are given in Table I. Figure 1 gives the openloop responses of tray 31 composition and vapor rate for a 10% step increase in reflux flow rate. All the models predict the same behavior for liquid flow rates and liquid Composition changes. However, vapor rates are different. Model 1 predicts an initial increase in vapor rates, a subsequent decrease, and later a final increase in vapor rates. This behavior is a consequence of having neglected the mass accumulation term and the enthalpy accumulation term in the energy differential equation. Model 2 predicts an initial decrease in vapor rate followed by an increase. The correct behavior, as given by model 3, should predict an increase in vapor rates. Tray compositions are increasing for the more volatile component, heats of vaporization decrease, and consequently vapor rates must increase. Figure 2 shows the openloop responses of tray 6 and tray 31 for a 10% step increase in the feed flow rate. Model 1 predicts an oscillatory behavior for vapor rates. The effects of these vapor rate changes are propagated to the upper section of the column, as shown by tray 31 vapor rate. Tray composition, responding to vapor rates change, also shows an oscillatory behavior on tray 31. The same comments used to explain the behavior for an increase in the reflux flow rate apply for an increase in the feed flow rate. Model 2 predicts a wrong initial decrease in vapor flow rates since the enthalpy accumulation term was neglected in the energy derivative. Later on, since feed tray composition is slowly increasing for the more volatile component and molar heats of vaporizationare decreasing, vapor flow rates increase. Model 3 predicts a reasonable increase in vapor flow rates.
0
5
5
t
7
0 9955' 0
60 TIME
120
180
(minutes)
Figure 3. Openloop response for an increase in feed composition from 50 to 60 mol % methanol. Notice that tray 31 composition decreases for this increase in feed rate. For equimolar overflow systems the opposite would be expected. Figure 3 gives the openloop responses for a step increase in feed composition from 50 to 60 mol % methanol. Model 1 predicts an initial decrease in vapor rates. Since feed tray composition is increasing and the heats of vaporization are decreasing, an increase in vapor flow rates must be the correct dynamic behavior, as predicted by model 3. Model 2 is a little closer to reality. Conclusions These openloop results show that the small computer time savings realized by using model 1or model 2 do not compensate for the abnormal responses obtained with these models. Appendix A Equilibrium Functions. Equilibrium data taken from the literature were fitted to polynominal functions to give
325
Ind. Eng. Chem. Fundam. 1982, 21, 325-326
HHzo = 10790 + 6.57T + 0.0036P
composition-dependent values of relative volatility and temperature. The following forms were used. a = a1 ugx a$ + u4x3
+
+ T = tl + t2x + t a x 2 + t4x3 + t5x4 + t6x5 CYX
Y = 1 + (a- 1)x
ai, ~
2 ~, 3 ~4,
= 7.450144, -11.22696, 8.477119, -2.27356
t4, t 5 , t 6 = 2 2.0, -285.649, 898.346, -1618.808, 1425.993, -483.521
tl, t 2 , t 3 ,
Vapor-liquid equilibrium data were taken from Van Zan’dicke and Verhoeye (1974), Uchida et al. (1953), Verhoye and Schepper (19731, Ocon and Rebolleda (1958), and Ocon and Taboada (1959). Appendix B B.l. Density Functions. Reference: Yaws (1977). The mixture molar density was found from X 1-x -1 =-
+-
Pmix pHz0
PMeOH
PHzO
= 0.05667 - 0.0000333T
P M ~ O H=
0.0253 - 0.0000313T
where T is in “C and p is in kg-mol/dm3. B.2 Enthalpy Functions. Reference: Touloukian and Makita (1970). The enthalpy of a mixture is obtained from the pure components enthalpy by neglecting the mixing effects. hmix = XhMeOH + (1- X)hHzO Hrnix
= YHMeOH + (l - Y)HHzO
The pure components enthalpy, in the range of 50 to 100 “C is given by h H z o = 17.46T + 0.0036P hMeOH
= 17.28T
HMeOH
= 9680 - 3.52T + 0.0384P
where Tis in “Cand enthalpies are in kcal/kg-mol. The pure components enthalpy was taken as zero at 0 OC. Nomenclature B = bottoms flow rate, kg-mol/min D = distillate flow rate, kg-mol/min hN = liquid enthalpy, kcal/kg-mol HN = vapor enthalpy, kcal/kg-mol LN = liquid overflow rate, kg-mol/min M B = bottoms accumulator holdup, kg-mol M D = reflux accumulator holdup, kg-mol M N = tray liquid holdup, kg-mol TN =.tray temperature, OC t = time, min VN = vapor rate, kg-mol/min x N = liquid composition, mole fraction yN = vapor composition mole fraction Greek Letters a = relative volatility pN
= liquid density, kg-mol/dm3
Literature Cited Doukas, N.; Luyben, W. L. Instrum. Techno/. 1978. 25, 43. McCune, L. C.; Gallier, P. W. ISA Trans. 1973, 72. Ocon, J.; Rebolleda, F. An. R . Soc. ESP. Fls. oukn.Ser. B 1958, 54, 525. Ocon, J.; Taboada, C. An. R . Soc.ESP. Fis. Quim. Ser. B 1959 55, 255. Toubukian, Y. S.; Makite, T. “Thermophysical Properties of Matter”; IFIIPlenum: New York, 1970; Vol. 6. Uchida, S.; Ogawa, S.; Hirata, M.; Shimada, Y.; Schimokawa, S. Kagaku KlkelCt”. Eng. 1953, 17, 191. Van Winkle “Distillation“; McGraw-Hill: New York, 1967. Van Zandljcke, F.; Vethoeye, L. Appl. Chem. Blotechnol. 1974, 24, 709. Vehoeye, L.; Schepper, H. de Appl. Chem. Blotechnol. 1973, 23, 607. Yaws, C. L. “physical Properties”; Mc(Law-Hili: New York, 1977.
Department of Chemical Engineering Carmelo Fuentes Lehigh University William L. Luyben* Bethlehem, Pennsylvania 18015 Received for review November 14,1980 Revised manuscript received August 14,1981 Accepted March 26, 1982
+ 0.0384P
Acentric Factor and the Critical Volumes for Normal Fluids An empirical equation relating the critical volume, V,, the acentric factor, w, the critical pressure, P,, the critical temperature, T,, and the gas constant, R , has been obtained. This equation accurately predicts the values of V , for numerous normal fluids.
Nath et al. (1976) have obtained an empirical equation relating the reduced vapor pressure, PI,the reduced temperature, TI, and the acentric factor, w (Pitzer et al., 1955), for normal fluids in the region of the coexistence of liquid and vapor phases. They have found that this equation predicts the vapor pressure and other equilibrium properties accurately for unassociated polar and nonpolar organic liquids. Further, Nath (1979) has obtained an empirical equation relating AHv/RT (where AHv is the heat of vaporization per mole at the temperature T (K) and R is the gas constant), T,, and w. This equation predicts the values of AHJRT for numerous unassociated polar and nonpolar organic liquids accurately. Recently Nath (1980) has reported an empirical equation which relates BPJRT, 0196-4313/82/1021-0325$01.25/0
(where B is the second virical coefficent, P, is the critical pressure, and T, is the critical temperature in degrees K) with TI and w. This equation has been found to predict quite accurate values of BPJRT, for numerous unassociated polar and nonpolar gases. In the present program an equation is proposed which can accurately correlate values of the critical volumes V , for normal fluids. Hence, the data on critical constants for a number of hydrocarbons, as reported by Kudchadker et al. (1968), have been examined, and the following empirical equation has been obtained
RTC V, = (0.2908 - 0.0825~)-
PC
0
1982 American Chemical Society
(1)
