Ind. Eng. Chem. Res. 1998, 37, 1063-1070
1063
SEPARATIONS Rigorous Simulation of Gas Absorption into Aqueous Solutions R. Cadours† and C. Bouallou* Ecole Nationale Supe´ rieure des Mines de Paris, Centre Re´ acteurs et Processus, 32 Bd Victor, 75015 Paris, France
A general model based on the film theory is used to set up the diffusion/reaction partial differential equations which describe the absorption of gases with multiple reversible reactions. This model includes the electrostatic terms that arise from differences in diffusivities. The system of equations made up of a set of partial differential equations and an algebraic equation is solved by finite-difference iterative method. The model developed in this paper, combined with appropriate literature parameters, will be used to model acid gas treating systems in the general case. Optimization using a Levenberg-Marquardt algorithm provides the rate coefficients for different kinetic mechanisms, from experimental data. The utility of this model is illustrated by applying it to the absorption of carbon dioxide into MDEA aqueous solutions. Comparison of different kinetic mechanisms shows a significant influence of the reaction between CO2 and OH-. The contribution of this reaction seems to be generally underestimated in the previous works using analytical or numerical methods. Introduction Processes in which mass transfer is enhanced by chemical reaction are frequently encountered in the process industry. For both energy and environmental purposes, acid gases must be removed from natural gas, as well as in petroleum refining, coal gasification, and hydrogen production. For such a purification, absorption with aqueous solutions of alkanolamines remains the dominant industrial technology. The design of absorbers and strippers requires information about transfer and kinetic mechanisms, and the related rate expressions. Much work has been done with respect to modeling and calculating the rate of mass transfer in absorption of gases followed by complex chemical reactions. In the cases of single reversible or irreversible firstorder reaction, instantaneous or not, analytical expressions based on film or penetration theory have been presented for the mass-transfer rate (Secor and Beutler, 1967). The major drawback of these models is that their validities are limited due to the assumption involved. When no analytical or approximate solution is available, it is necessary to use numerical methods to estimate the effect of reaction rate on absorption or to estimate a rate constant from absorption data. Cornelisse et al. (1980) have studied the simultaneous absorption of CO2 and H2S in a secondary amine. The CO2 absorption rate was enhanced by reaction of amine to form carbamate, and the H2S reaction with the amine was assumed to be instantaneous. They used a Newton-Raphson method to linearize the nonlinear reaction * Author to whom correspondence is addressed. Phone: 33 (1) 45 52 54 70. Fax: 33 (1) 45 52 55 87. E-mail: bouallou@ paris.ensmp.fr. † E-mail:
[email protected].
terms. This technique decoupled the simultaneous differential equations and allowed them to be solved as a set of linear equations. Littel et al. (1991) have used fundamentally the same method to describe the simultaneous absorption of H2S and CO2 in a solution of a primary or secondary amine and a tertiary amine. In the model presented by Glasscock and Rochelle (1989) for the absorption of CO2 into MDEA aqueous solutions, multiple reactions and the diffusion of ionic species were taken into account. The coupling between positive and negative ions was described by the NernstPlanck equation. For gas absorption accompanied by second-order reversible reaction and the CO2 absorption in MDEA, comparisons were made between the enhancement factors obtained for film theory, Higbie’s penetration theory, Danckwert’s surface renewal theory, and simplified eddy diffusivity theory. For the steadystate theories, the equations are transformed into a large set of coupled, nonlinear algebraic equations. For the unsteady-state theories, a spatial discretization results in a set of coupled ordinary differential/algebraic equations which are integrated through time by DDASSL (Petzold, 1983). Rinker et al. (1995) presented a study in order to determine the significant reactions involved in the absorption of CO2 into MDEAaqueous solutions, in particular the effect of including or neglecting the CO2/ OH- reaction on the estimation of the forward rate coefficient of the MDEA-catalyzed hydrolysis reaction. The method of lines was used to transform each partial differential equation into a system of ordinary differential equations which were then solved by DDASSL (Petzold, 1983). The same numerical method was used by Hagewiesche et al. (1995) to estimate the rate coefficient of the reaction between CO2 and MEA at 313 K from experimentally measured absorption rates of CO2 into blends of MDEA and MEA.
S0888-5885(97)00546-0 CCC: $15.00 © 1998 American Chemical Society Published on Web 01/24/1998
1064 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998
We developed a general model for the diffusion/ reaction processes. We have considered the most rigorous approach of taking into account the electrostatic potential gradient terms with unequal diffusion coefficients for the ionic species even if much more computational effort is required. This model allows us to consider reactions with finite rate and instantaneous reactions with respect to mass transfer. We assume reversibility, generalized reaction rate expressions, and generalized stoichiometry. The applicability of our model is not limited by the type of reactions modeled. The main advantage lives in the simpleness of the numerical treatment because no numerical method is needed to transform each partial differential equation into a system of ordinary differential equations. Optimization using a Levenberg-Marquardt algorithm provides the rate coefficients for different kinetic mechanisms, from experimental data. In this work, the model has been applied to investigate the discrepancies between different works reported in the literature for the CO2 absorption into MDEA aqueous solutions, especially at high temperatures.
Together with appropriate initial and boundary conditions, the set of equations gives a unique solution. Initial Condition. For the volatile species, 1 e i e NG, we consider linear profiles, taking into account the boundary conditions at the gas-liquid interface and their bulk concentrations. For the nonvolatile species, the concentrations are equal to their bulk concentrations.
(
Ci(x,0) ) 1 -
For volatile species
In eq 1, φ represents the electrostatic potential gradient, which couples the diffusion of ionic species. Under the assumption of dynamic electroneutrality and the Nernst-Einstein equation, the electrostatic potential can be expressed as a function of ion concentrations and ion diffusivities: NC
φ(x,t) )
∑ zqDq RT q)1 F
∂Cq(x,t) ∂x
NC
(2)
∑ zq2DqCq(x,t)
Boundary Conditions. (a) x ) 0. (i) For the volatile species, the boundary conditions at the gasliquid interface in the isothermal case were obtained by considering the mass-transfer rate in the gas near the interface equal to the mass-transfer rate in the liquid near the interface. This leads to the following boundary condition:
Ri is the production term: NR
Ri(x,t) )
|
Di ∂Ci kGiHi ∂x
q,j
x)0
Pi Hi
(6)
x)0
∀t > 0
(7)
(ii) For the nonvolatile species, the fluxes are equal to zero:
|
∂Ci ∂x
)0
(8)
x)0
(b) x ) δ. (i) For all chemical species i, the concentrations are equal to their liquid bulk concentrations:
Ci(δ,t) ) Ci,bulk
(9)
There are a great number of advantages in using nondimensional variables in numerical work, though it is not essential to do so. A good deal of arithmetic is involved in our problem. A whole set of solutions with different physical parameters can be obtained from one basic solution in nondimensional variables with considerable economy of computer time. The fundamental parameters are highlighted, and analogies with physically different systems become clearer. To obtain a dimensionless system, we consider the following reduced variables:
x δ
˜t )
tD1 2
δ
D ˜i )
Di D1
φ˜ )
φFδ RT
(10)
(3) For volatile species
For one of the ionic components, the partial differential equation is replaced by the static electroneutrality condition in order to maintain electroneutrality throughout the mass-transfer zone:
C ˜i )
For nonvolatile species C ˜i )
NC
ziCi(x,t) ) 0 ∑ i)1
+
|
∂Ci ∂x
where D1 is the diffusion coefficient of the first volatile component.
NC
λi,jkj∏CRq ∑ j)1 q)1
Ci(0,t) )
x˜ )
q)1
}
0 e x e δ (5)
Numerical Model
∂Ci(x,t) ) ∂t ∂Ci2(x,t) F ∂(φ(x,t) Ci(x,t)) + Ri(x,t) (1) - ziDi Di 2 RT ∂x ∂x
)
For nonvolatile species Ci(x,0) ) Ci,bulk
kGi(Pi - HiCi) ) -Di
In this work, the phenomenon of mass transfer accompanied by chemical reaction was modeled according to the film theory. We use a fictitious time to find a steady-state solution corresponding to the film theory. The effect of the electrostatic potential gradient on the diffusion of ionic species was also taken into account. This yields for each component i a partial differential equation:
x Pi + Ci,bulkx δ Hi δ
Hi Ci Pi Ci Ci,bulk
R ˜ i ) Ri
δ2Hi D iP i
δ2 R ˜ i ) RiD C i
i,bulk
}
(11)
(4) Equations 1 and 4 become
Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 1065
∂C ˜i ∂C ˜i ∂2C ˜i ∂φ˜ ˜ iφ˜ ˜ iC ˜i )D ˜ i 2 - ziD - ziD +R ˜ i (12) ∂t˜ ∂x˜ ∂x˜ ∂x˜ NC
ziC ˜i ) 0 ∑ i)1
(13)
Initial Condition. C ˜ i(x˜ ,0) ) 1 +
For volatile species
˜ i(x˜ ,0) ) 1 For nonvolatile species C
(C P H - 1)x˜ i,bulk
C ˜ i(0,t˜) )
|
∂C ˜i For nonvolatile species ∂x˜
x˜ )0
|
˜i Di ∂C kGiHiδ ∂x˜
)0
x˜ )0
+1
}
∀t˜ > 0 (15)
(b) x˜ ) 1.
Ci,bulkHi Pi For nonvolatile species C ˜ i(1,t˜) ) 1
For volatile species
}
∀x˜ g 0 (14)
Boundary Conditions. (a) x˜ ) 0. For volatile species
i
i
C ˜ i(1,t˜) )
}
max
1eieNC 1ejejmax
|
|
n-1 -C ˜ i,j n C ˜ i,j
e 5 × 10-6
solution of MEA at 298.15 K. The value of the equilibrium constant K is 275. The solubility of H2S is taken as 0.10 gmol‚L-1‚atm-1. The diffusivity of monoethanolamine is 0.64 times that of H2S. The reaction between H2S and MEA is instantaneous.
H2S + MEA S MEAH+ + HS-
∀t˜ > 0 (16)
Numerical Treatment. The production terms are nonlinear. Consequently, equations render a set of nonlinear coupled partial differential equations. An implicit centered finite-difference iterative method was employed to solve the system of coupled differential equations subject to the boundary and initial conditions given above. The implicit method ensures stability for the system of equations. From the initial concentration field specified above, the concentration gradients are obtained by numerical derivation. Then the electrostatic potential and its gradient are deduced. They enter as source terms in the transfer equations which are solved using the finite-difference iterative method, giving the solution for the next iteration. The calculations are repeated until the convergence criterion n C ˜ i,j
Figure 1. Concentration profiles for H2S absorption into MEA aqueous solution.
(17)
n is the concentration of species i at is satisfied. C ˜ i,j point j and at iteration n. The fictitious adimensional time was chosen in order to realize a compromise: it must be large enough to limit roundoff error and also sufficiently small for convergence of the small grid size. The steady-state solution is found by iterating upon an index (n). In this way, an efficient and consistent algorithm is obtained for the set of equations. The model developed in this work, combined with appropriate literature parameters, will be used to model acid gas treating systems in the general case. A Levenberg-Marquardt fitting procedure was used to infer kinetic rate constants from the experimental data.
Validation A first validation was established with the case of H2S absorption by monoethanolamine presented by Danckwerts (1970). H2S at 0.1 atm is absorbed in a 1 M
The equilibrium is available everywhere in the film, and the concentrations of the products HS- and MEAH+ are taken to be equal to one another everywhere. In this case, it is possible to get an analytical solution for the absorption rate:
φ ) kL(CH2S,int - CH2S,bulk) ×
[
1+
]
DMEAH+ β (18) DH2S CH2S,int - CH2S, bulk
with
β)
{[(
1 2
DHS-1 DMEAH+
]
)x
KCMEA,bulkCH2S,bulk +
[ ]} {(
DHS- DMEA + 4CH2S,intK C + DMEA DMEAH+ MEA,bulk 0.5 DHS1 KC C +1 × MEA,bulk H S,bulk x 2 2 DMEAH+ DHSKC C + C K (19) MEA,bulk H S,bulk x 2 DMEA H2S,int
DHSC K DMEA H2S,int
2
) }
We assume that the bulk liquid is free from dissolved H2S. We also assume that the diffusivities of the products are equal to that of the amine. Using the same parameters for numerical and analytical computations, we obtained the same result within a deviation of 0.6%: we obtain 4.97 × 10-2 mol‚m-2‚s-1 with the numerical model even though the analytical solution proposed by Danckwerts (1970) gives 5 × 10-2 mol‚ m-2‚s-1. As suggested by Glasscock and Rochelle (1989), we treated the instantaneous reaction between H2S and MEA by making rate constants large enough so that the reaction is in equilibrium everywhere in the boundary layer. The concentration profiles for each species are represented in Figure 1. These profiles clearly show the depletion of the rapid-reacting MEA near the gas-liquid interface.
1066 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 Table 1. Parameters for the MDEA System at 318 K Used by Glasscock and Rochelle (1989) parameter
value
reference
KI KII KIII KIV kI,f kII,b HCO2 DCO2 DMDEA DHCO3DCO3DOHloading [MDEA]T
132 4.18 × 10-8 kmol‚m-3 5.96 × 10-4 kmol‚m-3 nonindependent 10 m3‚kmol-1‚s-1 3.47 × 104 m3‚kmol-1‚s-1 50 atm‚m3‚kmol-1 1.62 × 10-9 m2‚s-1 0.75 × 10-9 m2‚s-1 0.94 × 10-9 m2‚s-1 0.7 × 10-9 m2‚s-1 4.5 × 10-9 m2‚s-1 0.005 molCO2/molMDEA 2.0 kmol‚m-3
Critchfield and Rochelle (1987) Critchfield and Rochelle (1987) Critchfield and Rochelle (1987) Critchfield and Rochelle (1987) Critchfield and Rochelle (1987) Astarita et al. (1983) Critchfield and Rochelle (1987) Versteeg (1986) Kigoshi and Hashitani (1963) Kigoshi and Hashitani (1963) Newman (1973)
Figure 3. Comparison between the concentration profiles from this work ([-]) and those from Glasscock and Rochelle (1989) ([-]*) for the MDEA system at high CO2 partial pressure (CCO2,int ) 10-2 kmol‚m-3).
model would give an underestimation of the influence of reaction II. With this mechanism, we checked that we did not introduce any perturbation by replacing for one of the ionic species the differential equation of diffusion with the algebraic equation of static electroneutrality. Results and Discussion
Figure 2. Comparison between the concentration profiles from this work ([-]) and those from Glasscock and Rochelle (1989) ([-]*) for the MDEA system at low CO2 partial pressure (CCO2,int ) 10-4 kmol‚m-3).
A second verification of the model was carried out by comparison with the numerical solution given by Glasscock and Rochelle (1989) for a four-reaction mechanism involving two finite rate reactions, reactions I and II, and two instantaneous equilibria involving a proton transfer, reactions III and IV.
CO2 + MDEA + H2O S MDEAH+ + HCO3- (I) HCO3- S CO2 + OH-
(II)
CO32- + H2O ) HCO3- + OH-
(III)
MDEA + H2O ) MDEAH+ + OH-
(IV)
At the conditions of the simulation, where the concentration of hydronium ion is very low, we have neglected the water dissociation reaction:
2H2O ) OH- + H3O+
(V)
We use the same parameters as Glasscock and Rochelle (1989) (Table 1). The relationship between the physical mass-transfer coefficient and the film thickness is
δ ) D1/kL
(20)
The results in terms of concentration profiles are in good agreement in the case of a low CO2 partial pressure (Figure 2). Some differences are observed between the two models at high CO2 partial pressure (Figure 3). They are explained by the choice of the convergence criterion (eq 18), which seems to be harder in our case. The hydroxide depletion is also overestimated by the model of Glasscock and Rochelle (1989), and then their
The emission of carbon dioxide is considered one of the biggest causes of global warming. Reduction of CO2 emitted as flux gases from a large-scale stationary source has become a worldwide problem. Removal of CO2 from process gas has been achieved in industry by gas absorption processes employing alkanolamines. From the viewpoint of economics, high-level regeneration of the absorbent is indispensable. Since the regeneration step includes heating of the absorbent, the energy requirement will significantly increase when treating a large volume of absorbent. Therefore, the reduction of heating energy becomes a determining factor for realizing gas absorption processes. CO2 absorption into alkanolamine aqueous solutions (MEA, DEA, MDEA, and their blends) has been extensively studied. However, there has been great confusion and disagreement in the literature concerning the magnitude of the rate coefficient of each reaction involved in the considered mechanisms for absorption in aqueous solutions of alkanolamines. The reaction between CO2 and aqueous MDEA has been the subject of many experimental studies. However, there are widely varying opinions in the literature regarding the mechanism of this reaction. In general, the reaction mechanism between CO2 and MDEA proposed by Donaldson and Nguyen (1980) is considered. This mechanism is a base-catalyzed hydration reaction which implies that MDEA do not react directly with CO2 (reaction I). The first investigations of the kinetics of this reaction were reported in the literature by Blauwhoff et al. (1984) and Barth et al. (1984). Their results were limited at 298 K. The determination of the kinetic parameters at various temperatures have been reported by Yu et al. (1985), Haimour et al. (1987), Critchfield (1988), Versteeg and van Swaaij (1988a), Toman and Rochelle (1989), Tomcej and Otto (1989), Littel et al. (1990), and Rinker et al. (1995), and recently Pani et al. (1997) reported a complete review of the data on the reaction between CO2 and aqueous MDEA. The significant discrepancies in the results presented by these authors can be explained by the significant effect of the reaction between CO2 and OH-, especially when
Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 1067
estimating the rate coefficient of reaction I for unloaded MDEA aqueous solutions. Rinker et al. (1995) have considered over and above reaction I reactions II-V. They have shown clearly that neglecting the CO2/OHreaction can result in large errors in the rate coefficients for the MDEA-catalyzed hydrolysis reaction, especially at high temperatures. Pani et al. (1997) have developed an efficient apparatus to measure absorption kinetics of acid gases into amine solutions. A thermoregular constant interfacial area Lewis-type cell was operated by recording the pressure drop during batch absorption. The experimental data obtained by these authors are well adapted to be treated with the film theory. For the calculation of absorption rate and also of kinetic parameters, the values of the diffusivities of all species and the Henry’s law constant for CO2 are required. The equilibrium constants of all chemical reactions are also needed. The diffusivity of CO2 in amine solutions is estimated with the N2O analogy presented by Versteeg and van Swaaij (1988b):
(DN2Oµ0.8)aq sol ) (DN2Oµ0.8)water
(21)
The diffusion coefficient of CO2 in the amine solution is then calculated using
( ) ( ) DN2O DCO2
)
aq sol
DN2O DCO2
(22)
water
The diffusion coefficients for N2O and CO2 in water are calculated according to the modified Stokes-Einstein relation (Versteeg and van Swaaij, 1988b):
(DCO2)water ) 2.35 × 10-6 exp(-2119/T) (DN2O)water ) 5.07 × 10-6 exp(-2371/T)
( )
DCO2 µaq sol 2.43 µwater
(24)
0.2
(25)
The diffusion coefficients of the ionic species are assumed to be equal to those of MDEA in order to compare our results with the literature results. The viscosities of water and MDEA aqueous solutions and the Henry’s law constant for CO2 are calculated from the correlations of Al-Ghawas et al. (1989). When reaction II is considered, the values of the backward rate were calculated from the correlation proposed by Pinsent et al. (1956) for the temperature range 273-313 K:
log(kII,b) ) 10.635 - 2895/T
4229.195 log(T) + 9.7384T - 0.0129638T2 + 1.15068 × 10-5T3 - 4.602 × 10-9T4 (27) A correlation for the determination of the values of (KIIKV) is obtained from the data reported by Read (1975) for the temperature range of 273-523 K:
log(KV/KII) ) 179.648 + 0.019244T 67.341 log(T) - 7495.441/T (28) The equilibrium constant of reaction III is obtained from the correlation of Danckwerts and Sharma (1966) over the temperature range of 273-323 K:
log(KV/KIII) ) 6.498 - 0.0238T - 2902.4/T
(26)
For all the equilibrium constants, we use the same correlations as previously used by Rinker et al. (1995). Olofsson and Hepler (1975) reported a correlation for the water dissociation constant KV for the temperature range of 293-573 K:
(29)
The correlation of Barth et al. (1984) established over the temperature range of 298-333 K is used for the calculation of the equilibrium constant of reaction IV:
log(KV/KIV) ) -14.01 + 0.018T
(30)
All these correlations are extrapolated up to 343 K, the greatest experimental temperature reached by Pani et al. (1997). The first time, reaction I is treated as a reversible second-order reaction and the effect of reactions II-IV is neglected (mechanism 1). This mechanism, proposed by Donaldson and Nguyen (1980), is generally presented in the literature as reasonable. It is a base catalysis of the CO2 hydration reaction. We have determined with our model the kinetic rate of the forward reaction I for Table 2. Identification of the Forward Rate Coefficients of Reaction I for Mechanisms 1 and 2 from the Data of Pani et al. (1997)
(23)
The diffusion coefficient of MDEA is given by the correlation proposed by Pani et al. (1997):
DMDEA )
log(KV) ) 8909.483 - 142613.6/T -
kI,f, m3‚mol-1‚s-1 T, K 295.95 295.85 296.15 295.95 295.75 295.25 296.45 296.45 296.45 296.45 296.45 296.45 296.45 296.45 296.45 317.75 317.55 318.35 318.15 317.65 343.45 343.45 343.55 342.55 342.65 342.65 342.55 343.55 343.55 342.75 342.95 342.25
CMDEA,T,
mol‚m-3
844.4 1272.6 1984.1 1984.2 2583.2 3477.4 4379.1 4379.1 4379.1 4379.1 4379.1 4379.1 4379.1 4379.1 4379.1 838.0 1966.6 2556.7 3435.8 4324.1 826.7 826.7 1245.0 1939.2 1939.0 1939.0 1939.2 1937.9 1937.9 2519.8 3381.4 4251.2
mechanism 1
mechanism 2
6.33 5.10 3.82 3.96 3.52 3.30 5.92 5.77 4.37 4.43 3.87 3.13 3.62 4.61 3.11 19.93 13.65 13.17 13.59 16.80 51.08 46.81 50.07 58.65 55.14 52.32 49.43 54.27 54.32 47.87 51.66 63.20
4.70 3.85 3.11 3.77 3.37 3.14 5.76 5.52 4.10 4.14 3.60 2.94 3.42 4.35 2.89 6.71 10.25 9.80 10.30 12.78 13.21 11.49 13.45 22.27 18.99 16.02 12.19 17.96 17.96 12.61 27.92 36.28
1068 Ind. Eng. Chem. Res., Vol. 37, No. 3, 1998 Table 3. Comparison with Energies of Activation from Previous Works ref
T, K
CMDEA, 103 mol‚m-3
PCO2, 105 Pa
E a, kJ‚mol-1
Yu et al. (1985) Haimour et al. (1987) Critchfield (1988) Versteeg and van Swaaij (1988a) Tomcej and Otto (1989) Littel et al. (1990) Rinker et al. (1995) Pani et al. (1997) this work
313-333 288-308 282-350 293-333
0.2-2.5 0.85-1.7 1.7 0.17-2.7
1 1 1
