I n d . Eng. Chem. Res. 1987, 26, 1092-1099
1092
ating mode which is characterized by the virtual absence of foam (lower holdups) is called the turbulent bubbling regime. This type of behavior has been observed for the first time in a system with molten paraffin wax as the liquid medium. The start-up procedure determines which one of the two regimes will be attained. A transition from the foamy to the turbulent bubbling regime occurs when ugexceeds a certain critical value, and transition from the turbulent bubbling to the foamy regime occurs as the superficial gas velocity drops below a certain critical value. These critical velocities are influenced by the temperature (i.e., the liquid viscosity) and the gas distributor type. Other factors which might influence these transitions, are column diameter, liquid static height, and impurities, but their effect has not been determined. Some of previous results for gas holdups that either were not well understood (e.g., the foam breakup reported by Farley and Ray (1964)) or appeared to be in a poor agreement (i.e., small vs. high gas holdup values obtained in different studies with different types of gas distributors) can be explained in terms of the existence of these two flow regimes. Bach and Pilhofer’s (1978) empirical correlation provides a very good fit for the gas holdups obtained in studies with orifice plate distributors in the absence of foam.
Acknowledgment This work was supported by DOE under Contract DEAC22-84PC 70027. We are grateful to William Deutchlander, Robert Drummond, Mike Noak, and Dr. Khan Nguyen-tien for their help with design and construction of the experimental apparatus.
Literature Cited Akita, K.; Yoshida, F. Ind. Eng. Chem. Process Des. Dev. 1973, 12, 76. Anderson, J. L.; Quinn, J. A. Chem. Eng. Sci. 1970, 25, 373.
Bach, H. F.; Pilhofer, T. Ger. Chem. Eng. 1978, I , 270. Calderbank, P. H.; Evans, F.; Farley, R.; Jepson, G.; Poll, A. Catalysis in Practice; Insitute of Chemical Engineers: London, 1963; p 66. Deckwer, W. D.; Louisi, Y.; Zaidi, A,; Ralek, M. I n d . Eng. Chem. Process Des. Dev. 1980, 19, 699. Deckwer, W. D.; Serpemen, Y.; Ralek, M.; Schmidt, B. Ind. Eng. Chem. Process Des. Deu. 1982, 21, 231. Farley, R.; Ray, D. J. Br. Chem. Eng. 1964, 9, 830. Gray, D.; Lytton, M.; Neworth, M.; Tomlison, G. Final report to the Department of Energy on Contract EE-77-‘2-01-2738, 1980. Hikita, H.; Asai, S.; Tanigawa, K.; Segawa, K.; Kitao, M. Chem. Eng. J . 1980, 8 , 191. Kolbel, H.; Ackerman, P.; Engelhardt, F. Proceedings of the Fourth World Petroleum Congress; Carlo Columbo Publishers: Rome, 1955; Section IV/C, p 227. Kuo, J. C. W. Final Report to the Department of Energy on Contract DE-AC-22-80 PC 30032, 1983. Kuo, J. C. W.; Di Sanzo, F. P.; Garwood, W. E.; Gupte, K. M.; Leib, T. M.; Smith, J. Quarterly report (1 Jan-31 March) to the Department of Energy on Contract DE-AC-22-83 PC 60019, 1984a. Kuo, J. C. W.; Di Sanzo, F. P.; Garwood, W. E.; Gupte, K. M.; Leib, T. M.; Smith, J. Quarterly report (1 Oct-31 Dec) to the Department of Energy on Contract DE-AC-22-83 PC 60019, 198413. Maruyama, T.; Yoshida, S.; Mizushina, T. Chem. Eng. Jpn. 1981,14, 352. Mersmann, A. Ger. Chem. Eng. 1978, 1 , 1. Pilhofer, T. Chem. Eng. Commun. 1980, 5, 69. Quicker, G.; Deckwer, W. D. Chem. Eng. Sci. 1981, 36, 1577. Riquarts, H. P. Ger. Chem. Eng. 1979, 2, 268. Schlesinger, M. S.; Crowell, J. H.; Leva, M.; Storch, H. H. Ind. Eng. Chem. 1951,43, 1474. Shah, Y. T.; Kelkar, B. G.; Godbole, S. P.; Deckwer, W. D. AIChE J . 1982, 28, 353. Smith, J.; Gupte, K. M.; Leib, T. M.; Kuo, J. C. W. Presented at the AIChE Summer National Meeting, Philadelphia, 1984. Thompson, G. J.; Riekena, M. L.; Vickers, A. G. Final report to the Department of Energy on Contract DE-AC-01-78ET 10159, 1981. Wallis, G. B. One Dimensional Tuo-Phase Flow; McGraw-Hill: New York, 1969. Zahradnik, J.; Kastenek, F. Chem. Eng. Commun. 1979. 3, 413.
Received for review February 3, 1986 Revised manuscript received February 6, 1987 Accepted February 28, 1987
Simulation of Continuous-Contact Separation Processes: Unsteady-State, Multicomponent, Adiabatic Absorption David M. Hitch,’ Ronald W. Rousseau,* and James K. Ferrell Department of Chemical Engineering, North Carolina State University, Raleigh, North Carolina 27695- 7905
A rigorous algorithm has been developed for the unsteady-state simulation of multicomponent, adiabatic absorption in packed columns. The simulation uses an implicit integration scheme to solve the partial differential model equations which may include both nonideal vapor-liquid equilibrium relationships and nonlinear mass-transfer expressions. All physical parameters used in the model are obtained from empirical correlations available in the literature. Simulation predictions, both dynamic and final steady-state, compared favorably with experimental results. The numerical simulation of continuous-contact separation devices is a topic that has long been neglected. In a previous article (Hitch et al., 1986), we pointed out that until recently most research efforts in the area of the design and analysis of separation processes have concentrated on *Author t o whom correspondence should be addressed a t School of Chemical Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0100. Current address: Eastman Kodak Company, Eastman Chemicals Division, Bldg. 152A, Kingsport, T N 37662.
staged operations. Coincident with the development of our SIMCAL algorithm, which is described by Hitch et al. (1986),other workers also produced numerical simulations that can handle rigorously the process of steady-state, adiabatic absorption in a packed column. Chief among these are the “nonequilibrium stage model” algorithms of Krishnamurthy and Taylor (1985). In fact, the nonequilibrium stage approach has been shown to be highly effective for steady-state simulation of a variety of both staged and continuous-contact separation processes. In contrast to the above-referenced work, the present article
0888-5885/87/2626-1092$01.50/0 0 1987 American Chemical Society
Ind. Eng. Chem. Res., Vol. 26, No. 6, 1987 1093 Table I. Model Variables 1. nc liquid component flow rates 2. nc gas component flow rate 3. nc liquid interfacial mole fractions 4.nc gas interfacial mole fractions 5. 1 liquid-phase temperature 6. 1 gas-phase temperature 7. 1 interfacial temperature total number of variables = 4nc + 3
k'yi
= kyi/(yf)i
1, g, x,* Y?
TL TG
t, =
discusses the development of a dynamic, packed-column simulator known as DYNCAL; this simulator should be of use in describing operating characteristics of packed columns of all types. It does not rely on any adjustable parameters, using instead available correlations for physical properties and packing characteristics.
Development of Descriptive Equations The mathematical model chosen to describe unsteadystate heat and mass transfer in an adiabatic, packed absorption column is substantively different from the model previously used to depict steady-state absorption. The main difference is in the incorporation of transient terms that describe the rates of accumulation and/or depletion of a component within gas and liquid phases. An analogous term was added to the overall energy balance to account €or the accumulation of energy within the liquid phase but was omitted for the vapor phase because vapors and gases have an exceedingly low heat capacity relative to the liquid phase. It is important to stress, however, that the commonly adopted assumption of negligible gas-phase holdup was not made when deriving the unsteady-state material balances; it will be shown that this term plays a significant role in describing the time-dependent behavior of the experimental system chosen for simulation. In addition, earlier model assumptions of Hitch et al. (1986) were relaxed so that the effect of nonideal vaporliquid equilibrium relationships and nonlinear masstransfer expressions on simulation results could be observed. The only remaining major restrictions imposed on the system equations are as follows: (1)The absorption column is adiabatic. (2) Pressure drop through the column is negligible compared to the overall column pressure. (3) Gas and liquid streams are in plug flow so that there are no radial temperature or concentration gradients. Given these restrictions, the mathematical formulation of the physical model may be stated in the following manner: unsteady-state individual comDonent material balances
ai, ag, (nc equations) at az uz at mass-transfer relationships liquid side -abLi + - = -abci --
(nc equations) gas side (nc equations) where k iiand k kiare mass-transfer coefficients "corrected through the use of the Wilke film factor (Wilke, 1950), which is intended t o account for the effect of net bulk movement of all transferring species on the mass-transfer rate of each individual component. For the gas phase, the Wilke film factor, yf, is computed for each component i as follows:
CN,/N,
,=I
Tc
3. equilibrium relationships yi* = &xi*
(nc equations)
where
Ki = f(T*,PJi*,Yi*) 4. unsteady-state, adiabatic energy balance
-~- & _ _H_ L-- ~ ( L H L ) ~ ( G H G )(1equation) at
az
az
where nc
HG
= EY~[CG~(TG - To) + Xi] i=l
nr
nc
G = Cgi i=1 nc
L = Eli i=l
i=l
5. heat-transfer relationships a. liquid side nc d(LHL) -- ~ L ~ [- T T*] L- iCN,ahi* =l
az
(I equation)
b. gas side
a(GHd --
az
nc
- h&[T* - TG]- CNiaHi* (1 equation) i=l
where N , is the mass flux of component i across the interface and is taken to be positive for absorption, i.e., transfer from the gas to the liquid phase. total number of equations = 4nc + 3 The above equations must be solved simultaneously to determine values of the dependent variables, listed in Table I, as functions of time and position uithin the column. In this study, equilibrium distribution coefficients, K,, were calculated either from solubility data presented by Lazalde-Crabtree et al. (1980) and Landolt-Bornstein (1976) or by using equations of state for both the liquid and vapor phase. The equations of state used in this work were the Soave-Redlich-Kwong (SRK) equation and the Peng-Robinson (PR) equation: binary interaction parameters were obtained from Chang (1984). Mass-transfer coefficients, k,, and kyl, and the specific interfacial area, a, were estimated by using co lations developed by Onda et al. (1968) or by Shuln 1 et al. (1955). Finally, heat-transfer coefficients, hL and ~ Z Gwere , determined by using the j-factor analogy (i.e., jH = j , (Treybal, 1980)), with the assumption that the same interfacial area applies for both heat and mass transfer. Heat-transfer coefficients estimated in this way are strictly valid only in the absence of mass transfer. For the gas
-
1094 Ind. Eng. Chem. Res., Vol. 26, No. 6, 1987
phase, a correction factor has been used to account for the effects of mass transfer (Ackerman, 1937). A similar factor does not exist for the liquid phase, but the magnitude of hLa is generally such that the effect of a correction for mass transfer would be negligible. Correlations giving both static liquid holdup (aLs) and operating liquid holdup (ah)as a function of packing size and type and the flow rate and physical properties of the solvent are given by Treybal (1980). The validity of these correlations for the particular type of packing material used in this study was verified through the use of data supplied by Kelly (1985). The total liquid holdup (@Lt) is defined as @Lt E @Lo
+
@La
(1)
The total volumetric vapor holdup may then be assumed to be the difference between the dry packing void fraction (c) and the total liquid holdup; i.e., @Gt
=
6
- @Lt
(2)
The value of E for a number of packing types and sizes may also be found in Treybal (1980). Once the total volumetric liquid and vapor holdups are known, the corresponding total molar holdups (BL and BG) are found through the relationships
(3)
time-dependent derivative terms
where C k again represents some dependent variable and the designation (N) represents an arbitrary "position" in the finite difference grid associated with the time domain. Therefore, the rate of change of variable Ck with time is approximated by the difference between the value of Ck at some time N + 1 and its value at time N divided by the time step At. This expression is analogous to the forward difference approximation for spatial derivatives. For example, by use of finite difference formulae to approximate both time and spatial derivatives in the unsteady-state material balance for component 1, bLi(J,N+1) - bLi(JtN) bGi(J,N+1) - b ~ i ( J f l ) At At L,(J+l,N+l)- l,(J,N+l) g,(J,N+l) - g,(J-l,N+l) +
-
AZ
(8) Approximating the solution to a system of partial differential equations in this fashion is generally referred to as an implicit integration method. This implicit scheme, also known as the backward Euler method, generates a system of algebraic equations that form a block-tridiagonal matrix:
(4)
. . .
where pL is the density and MLis the average molecular weight of the liquid solvent including any absorbed species. Equation 4 is based on the assumption of ideal gas behavior. (Although it was not done in the present paper, a compressibility factor could be estimated from an equation of state, SRK for example, and included in eq 4.) The individual component liquid and vapor molar holdups (bLi and bGJ may now be computed as bL,
1, = x~BL= ZBL
(5)
AZ
A(J,N)B(J,N)D(J,N 1
or
MxT=Z
Unsteady-State Simulation Algorithm Relaxation methods have been used extepsively in the analysis of unsteady-state staged columns. Such methods are highly stable because, in principle, they follow the actual transient behavior of a physically realizable start-up process. In addition to providing information about the transient outlet response to step changes in column feed composition, relaxation is a reliable way to converge to the final steady state. An important added benefit derived from the increased stability of the relaxation procedure is the ability of the algorithm to handle more complex mathematical models, including nonideal vapor-liquid equilibrium relationships and nonlinear mass-transfer expressions, without the convergence problems often encountered with other methods (Orbon, 1984). The numerical method chosen for the unsteady-state simulation uses finite difference expressions to approximate the derivatives in the differential model equations, thereby reducing the system of partial differential equations to a coupled system of algebraic equations. Forward and backward difference approximations to the spatial derivative terms were described earlier (Hitch et al., 1986); a similar expression may be used to approximate the
(10) where each of the terms in the coefficient matrix, M, is itself an NV X NV submatrix, NV being the number of dependent variables in the system of model equations. Similarly, each term in the solution vector, T , is itself a vector of dimension NV as is each term in the right-hand side vector, 2. The numerical procedure used to solve this block-tridiagonal matrix is analogous to that given previously for the steady-state simulation (Hitch et al., 1986) except that now the dependent variables are functions of both time and space. Therefore, one must decompose this matrix once for each time step to achieve an approximation to the actual individual component concentration and bulk temperature profiles over the entire column at that point in time. As the number of time steps taken increases, the solution relaxes toward the steady-state solution. This solution scheme is summarized in Figure 1. The accuracy of the solution, at any given time, depends on the magnitudes of both At and Az; smaller step sizes, of course, give the more accurate solution. The final steady-state solution for a given value of Az, however, is independent of the time step used. Therefore, if one is not interested in the time-dependent behavior of a particular problem, an arbitrarily large value of At may be used to converge quickly to the final solution, provided the value of At
Ind. Eng. Chem. Res., Vol. 26, No. 6, 1987 1095
0
Table 11. DAVO Absorber Operating Conditions
START
run
solvent temp, O F
DAVO-2”
11.4
51.5
76.0
NP,63
DAVO-3
16.3
58.3
75.0
CO,, 26 N2,66
DAVO-4n
13.9
64.9
77.5
Nz, 65
I I1 INPUT INITIAL CONDITIONS AT EACH END OF COLUMM AND SPECIFY COLUMN HEIQHT
CALL SUBROUTINE TO COMPUTE PHYIlCAL AND TRANSPORTPROPERTY PARAMETER VALUES AT TIME N
t
1
COz, 30
l.OO A
I
f
$
I
0.75
CO, 8
co, 10
~
-
-
e
Y
5
0.50
0.25
0”
-
-
-m
G
1 CONVERGENCE ACHIEVED?
-
logic flow diagram.
chosen is not so large as to introduce numerical instability into the calculation.
Experimental System The packed-absorber simulation program, DYNCAL, has been tested by comparison with data from the NCSU Coal Gasification Pilot Plant. The operation of the pilot plant and, in particular the Acid Gas Removal System (AGRS), has been described in detail elsewhere (Kelly, 1981; Staton, 1983; Hitch e t al., 1986), and a schematic diagram of the AGRS was given by Hitch et al. (1986). There are two possible modes of operation for the AGRS: (1)receiving the product gas from the fluidizedbed coal gasifier as a feed gas (referred to as integrated runs), or (2) using a specified syngas mixture as a feed (syngas runs). Both of these modes were used to obtain experimental data for the purpose of verifying the accuracy of DYNCAL predictions. Steady-state data were obtained from the seven integrated runs that had the best mass balance closures around the AGRS, Le., they were all within f2% of 100%. In each of these runs, refrigerated methanol was used as the solvent for absorption. Several syngas runs were made to provide transient data for use in model comparisons. These syngas runs are referred to as the Dynamic Algorithm verification Operations (DAVO runs). In these runs steady-state conditions were achieved within the AGRS using a feed gas of pure nitrogen, and then a step change in feed composition was effected by replacing the pure nitrogen feed with a syngas mixture. The absorber operating conditions and final feed compositions for the successful DAVO runs are given in Table 11. Both the total gas feed rate of -9 scfm and the liquid feed rate of -0.9 gpm were kept constant during these runs. The DAVO runs were purposely made at conditions designed to give significant amounts of C02 in the absorber product gas (sweet gas) so that experimental results could be compared more easily with DYNCAL predictions. It should be noted that column pressures were much lower and solvent temperatures considerably higher
+
t 0.00
DYNCAL
co. 7
-
E
0
Figure 1.
O F
‘Data Acquisition System Down.
INVERT COEFFICIENT MATRIX V I A NEWMAN’S METHOD
C
gas temp,
COZ, 25
CALCULATE ELEMENTS OF BLOCK TRIDIAQONAL COEFFICIENT MATRIX
COMPOSITION PROFILES AT TIME N + 1
gas feed composition, %
pressure, atm
I:
0.25
0.50
N2 flow rete CH4 flow rete 0.75
1.00
Measured Flow (Ib mole/h) Figure 2. DYNCAL vs. measured sweet gas flow for VLE behavior described by Henry’s law.
than those used in the integrated runs. The measured response in these runs was the composition of the sweet gas leaving the absorber. The increase in CO and C 0 2 in this stream from zero to steady-state values followed what are referred to here as breakthrough curves.
Simulation Predictions vs. Experimental Results Steady-State Results. The first comparisons between experimental results and DYNCAL predictions were of the amounts of certain components appearing in the absorber product gas at steady state in each of the seven integrated runs. Of the nine components typically present in the feed gas, only four were found in relatively large concentrations in the sweet gas. These four compounds (H2,CO, N2, and CH$ are only slightly soluble in refrigerated methanol and they typically comprise over 99% of the total sweet gas volume. Because of their very low solubility, a “pinched” condition exists a t the bottom of the column for these gases. Therefore, only equilibrium relationships have any impact on the predicted steady-state outlet concentrations of these components. Figure 2 shows the relationship between DYNCAL predictions and measured sweet gas flow rates of the four slightly soluble components for the seven integrated runs. Henry’s law was used to compute the equilibrium distribution coefficients (KJin these computations; agreements similar to those shown in Figure 2 were obtained when the SRK equation of state was used to estimate distribution coefficients. Because the integrated runs were made at favorable absorption conditions (high pressure and low temperature), only minute amounts of C02, H2S, and COS were detected in the sweet gas, and the measurements agreed with DYNCAL predictions of almost zero sweet gas concentrations for these three components. However, even though DYNCAL correctly simulates a situation where essentially none of the highly soluble acid gases appear in the sweet gas, this
1096 Ind. Eng. Chem. Res., Vol. 26, No. 6, 1987 Table 111. Measured vs. DYNCAL-Predicted Steady-State CO and N, Sweet Gas Flow Rates DAVO-3 DAVO-4 component VLE" type sweet gas flow rate, lb.mol/h % dev. sweet gas flow rate, lb.mol/h 0.142 f -3% co measd 0.112 f -3% -3.6 0.136 co Henry 0.108 0.0 0.141 co SRK 0.112 0.112 0.0 0.140 co PR N2 measd 0.927 f -3% 0.933 f -3% -2.0 0.899 N2 Henry 0.908 N2 SRK 0.891 -3.9 0.885 N:! PR 0.878 -5.3 -0.874
% dev.
-4.2 -0.7 -1.4 -3.6 -5.1 -6.3
Vapor-liquid equilibrium relationship used (measd = pilot plant data, Henry = Henry's law, SRK = Soave-Redlich-Kwong equation of state, PR = Peng-Robinson equation of state).
Table IV. Measured vs. DYNCAL-Predicted Steady-State Rate of C O , Absorption DAVO-3 Wilke film VLE" type MTCb type factor used? COP absorbed, lb.mol/h % dev. measd measd measd 0.326 & -3% +4.2 Onda no 0.340 Henry +4.5 Onda yes 0.341 Henry +0.8 Onda no 0.329 SRK +1.1 Onda yes 0.330 SRK +4.8 Onda no 0.342 PR +5.4 Onda yes 0.344 PR -1.1 Shulman no 0.323 Henry +0.2 Shulman yes 0.327 Henry -4.1 Shulamn no 0.313 SRK -2.9 Shulman yes 0.317 SRK -0.7 Shulman no 0.324 PR +0.5 Shulman Yes 0.328 PR
DAVO-4 CO, absorbed, lb.mol/h 0.258 f -3% 0.277 0.277 0.270 0.270 0.286 0.286 0.269 0.271 0.262 0.264 0.276 0.278
70dev. f7.4 +7.4 +4.7 f4.7
+10.9 +10.9 +4.3 +5.1
+1.6 12.4 +7.0
+7.8
" Vapor-liquid
equilibrium relationship used (Henry = Henry's law, SRK = Soave-Redlich-Kwong equation of state, PR = Peng-Robinson equation of state). Mass-transfer coefficient correlation used.
is no guarantee that the algorithm would predict accurately the sweet gas flows of these constituents under less favorable absorber conditions. Therefore, in order to establish the validity of the simulation under conditions where there is incomplete acid gas absorption and, also, to judge the accuracy of DYNCAL predictions of the transient response to a step change in feed composition, the DAVO syngas runs were made. As with the integrated pilot plant runs, operating conditions for the DAVO syngas runs were such that the absorption of the two slightly soluble gases, Nzand CO, was equilibrium-controlled. Therefore, the DYNCAL predictions of the rate of absorption of these two components did not depend on the type of mass-transfer coefficient correlation used. Table I11 gives a comparison of the measured steady-state sweet gas flow rates for Nzand CO with the flow rates predicted by DYNCAL, for runs DAVO-3 and DAVO-4 using three different types of vapor-liquid equilibrium relationships. The predicted flow rates differ from the measured results by an amount that is generally less than the experimental error of about 3% associated with pilot plant flow measurements. In contrast to the slightly soluble gases, C 0 2 flow rates predicted by DYNCAL were highly dependent on the mass-transfer expression. In fact, the sensitivity of the predicted COz flow to changes in the mass-transfer coefficient is as great or greater than it is for changes in the type of VLE relationship used, thus confirming the absence of a pinch region in the column for this highly soluble component. Table IV gives a comparison of measured C 0 2 absorption rates with those predicted by DYNCAL using a variety of vapor-liquid equilibrium relationships and mass-transfer coefficient estimation procedures. Several features are evident from these results, with the most obvious being the consistent slight overprediction of C 0 2 absorption rate. This might be due to small deviations from adiabatic conditions in the absorber (whose tem-
perature is below ambient) which would result in slightly higher column temperatures than are predicted and, concomitantly, lower solubilities and absorption rates in the experimental system. Aside from the tendency to overpredict the rate of C 0 2 absorption, there are a number of trends in the DYNCAL predictions that should be mentioned. First, the COz absorption rates predicted by using the Shulman correlation to estimate mass-transfer coefficients were always considerably lower than the rates calculated using masstransfer coefficients estimated via the Onda correlation. Gas-phase mass-transfer coefficients estimated from the Onda correlation were a factor of 3 larger than those obtained from the corresponding Shulman correlation; on the other hand, liquid-phase mass-transfer coefficients obtained from the two correlations were about the same. Second, differences in the predicted steady-state C 0 2 absorption rate resulting from a change in VLE relationship were also quite significant; the SRK equation of state always gave the lowest rates of COz absorption, for a given type of mass-transfer coefficient, followed by Henry's law and then the PR equation of state. Finally, incorporation of the Wilke film factor into the mass-transfer rate expressions increased the estimated rate of COz absorption. This behavior is expected because the Wilke factor accounts for the effect of the net molar transport of materials between phases on the concentration dependence term of the mass-transfer rate equations; since the net bulk movement of components for absorption is from the gas phase to the liquid phase, the predicted absorption rate of C 0 2 should always be enhanced by the inclusion of the film factor term. It is impossible with such limited data to determine which is the most general and appropriate approach for estimating vapor-liquid equilibrium relationship and mass-transfer coefficient correlations. In any case, as shown in Table IV, deviations between predicted and
Ind. Eng. Chem. Res., Vol. 26, No. 6, 1987 1097 1.oo
s
0
5
0.75
$ v)
I
t
C
m
0
0.0
10.0
20.0
Experimental data
30.0
j
. -
.-!A
Ti/
I
0 0.25
8
40.0
Time ( m i d Figure 3. Input CO, concentration fit by COP = Bo exp(-E,/t2).
0.50
0.0
10.0
0
DAVO-2,
P = 11.14 atm
A
DAVO-3.
P
DAVO-4,
P = 13.93 atm
20.0
= 16.33 atm
30.0
40.0
Time ( m i d
Figure 5. Normalized experimental breakthrough curves for CO in sweet gas.
6
-c
" t 0.15
P
,No
vapor holdup
i
Y
/ /
.-
=
a
.-
8
00 .
10.0
20.0
30.0
t
050
IT 1
1
ii
Measured Data
tI /
i
/
ON
-1
P
40.0
0
Time ( m i d
F i g u r e 4. Effect of neglecting vapor-phase holdup on DYNCAL breakthrough predictions.
measured rates of C 0 2 absorption were small (less than 11%in all cases, and less than 5% for most cases). The combination that provided the best agreements between measured and predicted values for these two runs was the SRK equation of state matched with Shulman masstransfer coefficients which, when coupled with the Wilke film factor, resulted in less than 3% deviation for both runs. Unsteady-State Results. Although the experimental plan called for introduction of a step change in the concentration of COz in the feed gas, the actual change, shown in Figure 3 for DAVO-4, was somewhat different. The data shown were fit by a nonlinear two-parameter equation that was used as a forcing function in the DYNCAL algorithm to simulate the time-dependent changes in feed gas composition. The transient response of the absorber to a change in feed composition can be detected readily by observing the composition of sweet gas as a function of time. Among the factors affecting the transient sweet gas composition are solute solubilities, rates of mass transfer, magnitudes of the liquid and gas flow rates and holdups within the column, and axial dispersion. As liquid and gas feed rates were the same for all DAVO runs, these factors are to be excluded from the following analysis. A key factor that has been neglected generally in modeling gas-liquid separation processes is the effect of gasphase holdup. The importance of including a gas-phase holdup term in the system of model equations that describe our pilot-scale apparatus is demonstrated in Figure 4, which shows the DYNCAL predictions of the breakthrough curve for CO in run DAVO-3 compared to the experimental breakthrough curve and a curve produced by setting the gas-phase holdup equal to zero. Clearly, gas holdup plays a significant role in the dynamics of this
0.25
0"
0
DAVO-2,
A
DAVO-3.
+ DAVO-4,
0.0
10.0
P = 11.14 atm
P = 16.33 a t m
30.0 P = 13.93 a i m 40.0
20.0
Time ( m i d Figure 6. Normalized experimental breakthrough curves for CO, in sweet gas.
relatively high-pressure system. Since the total molar gas-phase holdup is directly proportional to the column operating pressure, a much lower pressure would reduce the effect of gas-phase holdup, and this term might be omitted from the descriptive equations without loss of accuracy. The normalized CO and COz breakthrough curves in Figures 5 and 6 were generated to assist in evaluating the influence of several of the other above-mentioned factors on the transient response of packed absorption columns. As shown in Figure 5 , there was little difference in the CO breakthrough curves for the three DAVO runs, which is expected for a slightly soluble component. Additionally, the sharpness of the breakthrough indicates that axial dispersion is not a significant role in the dynamics of the process. There is, however, a slight but noticeable lag between the step change ( t = 0.0) and the initial breakthrough of CO into the sweet gas. Because of the very slight solubility of CO in methanol, this lag time must be associated with the gas-phase holdup alone. In contrast to the uniformity of the CO breakthrough curves, Figure 6 clearly shows that as the pressure in the absorber increases, the lag to initial C02 breakthough increases and the slope of the breakthrough curve decreases. Such behavior means that as the absorber pressure increases, return of the column to steady state after an upset requires a longer period of time. This behavior occurs because the solubility of C 0 2 in methanol increases dramatically as the operating pressure of the column is raised. Since the liquid- and gas-phase driving forces are greatest immediately after a step increase in COz feed rate, masstransfer rates from the gas to the liquid are very high. This, coupled with the enhanced capacity of the solvent to absorb the COi a t higher pressures, results in a longer
1098 Ind. Eng. Chem. Res., Vol. 26, No. 6, 1987 0.12
.c
1
Henry. Shulman. Wilke
5 -
E
0.09
I
I
I
I
I
DAVO-4 Step Change
c
' 7 1
-I
, i I
0.0
Time ( m i d
g
-
r
0.09 Jb.K,
I
I
S
h
u
I
I
l
m
:
I
I
I
,
1
0
Y
/
(II
m r
Henry, Shulman
I
I
SRK. Onda
j
0.0
10.0
20.0
I
1
I
40.0
I
60.0
d 80.0
Time ( m i d
Figure 7. DYNCAL-predicted vs. experimental CO, breakthrough curves for run DAVO-3. 0.12
20.0
Measured Data
30.0
1 40.0
Time ( m i d
Figure 8. DYNCAL-rredicted vs. experimental CO, breakthrough curves for run DAVO-4.
lag time required before COz begins to appear in the sweet gas. As time passes, and more C 0 2 is absorbed, the liquidand gas-phase driving forces decrease, as does the rate of mass transfer, thereby slowing the rate of saturation of the liquid holdup within the column. Consequently, the time required to reach a steady-state sweet gas composition is lengthened. At lower pressures, the liquid holdup becomes saturated more quickly and steady state is achieved more rapidly. The parameters used in DYNCAL to compare the predicted transient with experimental data were those that gave the best steady-state COz absorption rates. For DAVO-3 these were the Henry's law VLE relationship and Shulman mass-transfer coefficients with the Wilke film factor. Another combination that gave good results for this run was the SRK equation of state and Onda mass-transfer coefficients without the Wilke film factor. Figure 7 compares measured COz sweet gas flow rate for DAVO-3 with the DYNCAL predictions. Although both sets of parameters gave a satisfactory comparison for the steady-state sweet gas flow of COz, the Onda mass-transfer coefficients gave a better match to the experimental transient of the COz breakthrough curve than did coefficients from the Shulman correlation. For run DAVO-4, the best combination of parameters for fitting steady-state behavior was the SRK equation of state and the Shulman correlation without the Wilke film factor. The DYNCAL prediction of the C 0 2 breakthrough curve for DAVO-4 using this set of parameters is shown in Figure 8 along with the measured response. Also presented in this figure are plots of the DYNCAL-predicted breakthrough curves obtained by using the SRK equation of state combined with Onda coefficients and Henry's law combined with Shulman coefficients. Note that even
Figure 9. DYNCAL simulation of C 0 2 breakthrough curve for hypothetical process upset.
though two different types of mass-transfer coefficients were used, both the transient and steady-state portions of these two curves have approximately the same appearance, in contrast to the two breakthrough curves plotted in Figure 7. As a possible explanation of this behavior, recall that the gas-phase mass-transfer coefficients calculated by using the two empirical relationships differ by a factor of roughly 3 to 1,while the liquid-phase coefficients are approximately the same. In run DAVO-3, the pressure and, consequently, the solubility of COSin methanol were higher than in run DAVO-4. When the solubility of a solute in a solvent increases, the mass transfer between the gas and the liquid becomes more gas-phase-controlled; i.e., the rate of mass transfer is affected more by the magnitude of the gas-phase mass-transfer coefficient. Since the gas-phase coefficients computed by the two correlations are significantly different, we would expect different predictions of the breakthough curve for high pressure runs. In contrast, at lower pressure the mass-transfer rates become more liquid-phase-controlled and the magnitude of the liquidphase mass-transfer coefficient becomes more significant. Because there is good agreement between the liquid-phase coefficients for the two correlations, the predicted breakthrough curves should be more similar in appearance for lower pressure runs such as DAVO-4. The point should be made again that no adjustable parameters have been used in the model equations to force the DYNCAL-predicted transient response to fit the measured data points. Instead, all model parameters, including vapor and liquid holdups, have been estimated through the use of correlations available in the literature. Taking this into consideration, the DYNCAL predictions of the dynamic response of the sweet gas flow, depicted in Figures 7 and 8, to the step change in feed composition are very good approximations to the measured response. Finally, to illustrate the ability of the algorithm to react to rapid changes in absorber conditions in a stable fashion, DYNC& was used to simulate a hypothetical process upset. Using the SRK equation of state and Onda mass-transfer coefficients, DYNCAL first simulated the step change in feed conditions experienced in run DAVO-3. After the algorithm had converged to the steady-state for these conditions, the feed composition and temperature as well as the column pressure were abruptly altered to match those of run DAVO-4, and the simulation was continued until a new steady-state was achieved. The DYNCAL-generated C 0 2 breakthrough curve for this hypothetical situation is shown in Figure 9. Note that the final steady-state COz sweet gas flow rate predicted by DYNCAL for the upset from DAVO-3 to DAVO-4 conditions is the same as the steady-state flow predicted for a step change to DAVO-4
Ind. Eng. Chem. Res., Vol. 26, No. 6, 1987 1099 conditions from start-up. The rapid jump of the sweet gas flow rate occurring immediately after the upset is consistent with the flashing of COz from the solvent which would be expected in the event of a sudden drop in column pressure. Conclusions The development of DYNCAL represents a considerable improvement in the ability to simulate both steady-state and transient absorption in packed columns. The algorithm accommodates a rigorous set of model equations that may include both nonideal phase equilibrium relationships and nonlinear mass-transfer expressions. In addition, the relaxation procedure used to converge to steady-state is very stable and convergence is assured for any physically realizable situation. A very good approximation to the transient behavior of the system may be obtained by using relatively small time steps, At, in the implicit integration procedure; on the other hand, if knowledge of the transient behavior is not required, rapid convergence to steady-state conditions may be achieved by using relatively large values of At. Accuracy of the simulation depends on the validity of correlations used in phase equilibrium computations and in estimating transport properties. No system-specific adjustable parameters were used in the present work; however, adjustable parameters may be readily included should a higher level of accuracy be required than can be realized through the use of general empirical correlations. Although development and testing of the algorithm focused on absorption, the approach should be valid in simulating other continuous-contact operations such as distillation, stripping, and extraction. Several aspects of modeling absorbers of the type used in the experimental portion of this study also should be mentioned here: (1) The effect of gas holdup on system dynamics cannot be ignored at pressures in the range of the present study. (2) Differences between the empirical correlations of Shulman et al. (1955) and Onda et al. (1968) become more significant at conditions under which the gas phase is controlling the rate of mass transfer. (3) There was generally good agreement between steady-state experimental data and model predictions regardless of the correlations used in estimating transport rates; on the other hand, the accuracy of the predictions of transient behavior depended strongly on which correlations were used. Acknowledgment We gratefully acknowledge support of this research effort by the Environmental Protection Agency under Cooperative Agreement CR-809317. Nomenclature a = interfacial area per unit volume of packing, area/volume bG, = molar gas-phase holdup of component i, moles/volume bL = molar liquid-phase holdup of component i, moles/volume BG = total molar holdup of gas phase, moles/volume B t = total molar holdup of liquid phase, moles/volume CG1= heat capacity of i in gas phase, energy/(mole-degree) CL, = heat capacity of i in liquid phase, energy/(mole-degree) g, = molar flow rate of i in gas phase, moles/(time-area) G = total molar flow rate of gas phase, moles/(time-area) h,* = partial molar enthalpy of component i in the liquid at concentration x,* and temperature T* h ~= ’ gas-phase heat-transfer coefficient corrected for mass transfer, energy/(time-area-degree) hL = liquid-phase heat-transfer coefficient, energy/ (timearea-degree)
Hi* = partial molar enthalpy of component i in the gas at concentration yi* and temperature T* HG = enthalpy of the gas phase, energy/mole HL = enthalpy of the liquid phase, energy/mole AHsi = heat of solution of i, energy/mole in solution kZi= mass-transfer coefficient for component i in liquid phase, moles/ (time-interfacial area) kYi= mass-transfer coefficient for component i in gas phase, moles/(time-interfacial area) k i i = Wilke film factor corrected mass-transfer coefficient for component i in liquid phase, moles/(time-interfacial area) k = Wilke film factor corrected mass-transfer coefficient for component i in gas phase, moles/(time-interfacial area) Ki = Henry’s law constant li = molar flow rate of i in liquid phase, moles/(time-area) L_ = total molar flow rate of liquid phase, moles/(time-area) ML = average molecular weight of liquid phase, mass/mole Ni= molar flux of i across the gas-liquid interface, moles/ (time-interfacial area) nc = number of components in system P = total pressure, atm R = ideal gas constant, volume-atm/ (mole-degree) TG = bulk temperature of gas phase, degrees TL = bulk temperature of liquid phase, degrees To = reference temperature, degrees T* = interfacial temperature, degrees x i = mole fraction of i in liquid phase xi* = mole fraction of i in liquid-phase interface yi = mole fraction of i in gas phase yi* = mole fraction of i in gas-phase interface z = position in column, length
;i
Greek Symbols
= dry packing void fraction, volume/volume Xi = latent heat of vaporization of i, energy/mole @‘ct = total volumetric holdup of gas phase, volume/volume aL0= operating liquid-phase holdup, volume/volume @La = static liquid-phase holdup, volume/volume QLt = total volumetric holdup of liquid phase, volume/volume pL = average density of liquid phase, mass/volume t
Literature Cited Ackerman, G. Forschungsheft 1937,382,l. Chang, T. Ph.D. Thesis, North Carolina State University, Raleigh, 1984. Hitch, D. M.; Rousseau, R. W.; Ferrell, J. K. Znd. Eng. Chem. Process Des. Dev. 1986,25,699. Kelly, R. M. Ph.D. Thesis, North Carolina State University, Raleigh, 1981. Kelly, R. M., The Johns Hopkins University, Baltimore, personal communication, 1985. Kelly, R. M.; Rousseau, R. W.; Ferrell, J. K. Znd. Eng. Chem. Process Des. Dev. 1984,23,102. Krishnamurthy, R.; Taylor, R. AZChE J. 1985,31, 449. Landolt-Bornstein, Zahlenwarte und Funktionen 6. Auflage ZV. Band Technik, 4. Teil Warmetechnik Bestandteil C I , Absorption in Flussigkeiten mit niedrigem Dampfruck.; Springer-Verlag: Berlin-Heidelberg-New York, 1976. Lazalde-Crabtree, H.; Breedveld, G. J. F.; Prausnitz, J. M. AZChE J . 1980,26,462. Onda, K.; Takeuchi, H.; Okumoto, Y. J. Chem. Eng. J p n . 1968,1, 56. Orbon, J. J . M.S. Thesis, North Carolina State University, Raleigh, 1984. Shulman, H. L.; Ullrich, C. F.; Proulx, A. A,; Zimmerman, J. 0. AZChE J. 1955,1,253. Staton, J. S. M.S. Thesis, North Carolina State University, Raleigh, 1983. Treybal, R. E. Mass-Transfer Operations, 3rd ed.; McGraw-Hill: New York, 1980. Wilke, C. R. Chem. Eng. Prog. 1950,46, 95. Received for review February 24, 1986 Revised manuscript received January 27, 1987 Accepted February 28, 1987
