Equilibrium Approach for CO2 and H2S Absorption with Aqueous

Jun 7, 2013 - Acid gas removal from gas streams is mainly executed industrially via aqueous alkanolamines absorption, whose modeling often involves io...
3 downloads 0 Views 7MB Size
Article pubs.acs.org/IECR

Equilibrium Approach for CO2 and H2S Absorption with Aqueous Solutions of Alkanolamines: Theory and Parameter Estimation José Luiz de Medeiros,* Leandro Chagas Barbosa, and Ofélia de Queiroz F. Araújo Escola de Química, Federal University of Rio de Janeiro, Av. Horácio Macedo, 2030, Centro de Tecnologia, E - Ilha do Fundão, Rio de Janeiro - RJ, 21941-909 Brazil S Supporting Information *

ABSTRACT: Acid gas removal from gas streams is mainly executed industrially via aqueous alkanolamines absorption, whose modeling often involves ionic nonequilibrium chemical reactions and mass/heat interfacial transfers. This demands high computational efforts and often nonavailable physical parameters. Also, the role of ions, created by weak dissociations, is not well understood, and despite the high pressures, concentrations, and loadings, idealities are often recruited, such as ideal laws for gas and mass action. In this work, these issues are circumvented by prescribing only molecular species incorporated into a chemical theory equilibrium framework using full thermodynamics via cubic equations of state. Nonvolatile molecular Complex species are formed via equilibrium reactions involving acid gas, water, and alkanolamines (AGWA) systems. The approach is calibrated with large sets of AGWA data, partial pressures, and loadings, by searching equilibrium constants via maximum likelihood nonlinear implicit parameter estimation. These constants are used to address ideal gas properties of Complex species, which are necessary to model energy effects in AGWA separations.

1. INTRODUCTION The removal of acid gasesH2S and CO2from gas streams is probably the most required supporting operation in modern petroleum and natural gas (NG) industries. In part, this is due to the intensification of environmental constraints, combined with the tendency to process crudes with increasing C/H, S/H ratios1 and NG with high CO2/CH4 and H2S/CH4 ratios. Additionally, acid gas removal is also a critical step in industries of basic chemicals such as H2, SYNGAS, ethylene oxide, NH3, and urea.2 Several classes of technology are available for these purposes, of which the most relevant are the following: (i) absorption; (ii) adsorption on solid; (iii) cryogenic separation; (iv) membrane permeation; (v) membrane contactors; and (vi) chemical conversion. Among these, regenerative absorption is widely applied in high capacity gas treatments. A subclass in this category is absorption with aqueous alkanolamines, which has been used for more than 50 years for H2S/CO2 removal in oil and NG industries.3 Alkanolamines are low cost liquids easily produced in the ethylene oxide industry. This qualifies alkanolamine absorption as ready for large scale application on conventional coal power plants4 in postcombustion CO2 capture. Nevertheless, the escalation of this technology relies on continuous data gathering for systems composed by acid gas + water + alkanolamine, or AGWA, systems. An AGWA system may include one or more acid gases and one or more alkanolamines, but the presence of water, in high molar proportion to the amines, is essential for alkalinity. An AGWA system is normally in a reactive vapor−liquid equilibrium (RVLE) condition, under pressures from a few millimeters of mercury up to 70 bar and temperatures up to 180 °C. Primary [monoethanolamine (MEA)], secondary [diethanolamine (DEA)], tertiary [methyldiethanolamine (MDEA)], and hindered [2-amino-2-methylpropanol (AMP)] are common alkanolamines in AGWA systems. MEA is a benchmark © 2013 American Chemical Society

cosolvent with good characteristics for CO2 capture, such as absorption capacity, fast kinetics, water solubility, and low cost.5 In spite of this, MEA aggregates stripping concerns such as the following: high heat load per unit of gas, corrosiveness, degradation, and vapor losses.1 DEA is less reactive than MEA but demands less heat per unit of gas with less degradation, corrosiveness, and volatility. MDEA has the lowest reactivity with CO2, but it has some comparative advantages: highest equilibrium absorption capacity, lowest degradation, and heat load per unit of gas, H2S/CO2 selectivity, and negligible vapor losses. Alkanolamines with steric hindrance can improve absorption capacity or rate.6,7 This is the case of AMP, with absorption capacity and stripping heat similar to those of MDEA, but faster CO2 reaction.8 Blends of tertiary + secondary/primary alkanolamines for conjugating desirable items (and mitigating undesirable ones) have been practiced.9 Compared to single amines, the blends MDEA + MEA/DEA are better in terms of absorption, corrosiveness, and heat load,10 whereas MDEA + DEA + AMP can improve loadings.11,12 AGWA models in the literature have a tendency to approach absorption-stripping within rate-based formalisms13−31 where the reactive two-phase AGWA system is described via nonequilibrium interfacial mass/heat transfers, coupled to nonequilibrium and chemical equilibrium (ChE) ionic reactions in the aqueous phase, invariably with weak ion speciese.g. MEA·H+, DEA·H+, MDEA·H+, HOC2H4NHCOO−, HCO3−, HS−, H3O+, OH−, and (C2H4OH)2NCOO−. Even the atomic structure of some of these ions and the chemical bonds Received: Revised: Accepted: Published: 9203

September 20, 2012 May 25, 2013 June 7, 2013 June 7, 2013 dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

The point is that engineering and optimization of AGWA operations can be addressed if equilibrium limiting cases are correctly modeled under mass/energy/momentum conservation and thermodynamic consistency. In this context, another lineage of papers addressed RVLE models in AGWA systems32−37 at low pressure, where electrolytes play central roles and ionic excess models are used for activity coefficients.29−31 However, the inexorable deficiency of all RVLE ionic approaches is that ion excess models are not well suited to high concentrations, a basic truth in industrial AGWA separations. Another drawback is the unreliable prediction of thermal properties for the ionic liquid phase, such as enthalpy and reactive thermal effects, which are strong in industrial AGWA operations and totally related to the liquid transformations. Thermodynamic consistency means that the free course of equilibrium shifts is observed in response to changes in variables such as (T, P). As in distillation, this should be modeled within process simulators, which are suitable for solving separation trains. This demands full access to the thermodynamic properties of both phases in wide ranges of (T, P) and compositions. Currently, this has been achieved only for oil, NG, and petrochemical processes, all comprising molecular scenarios where the SRK/PR EOS38 can be used safely. This work intends to fill exactly this lacuna that was left by both castes of models: rate-based and RVLE ionic. It is proposed a RVLE molecular framework is capable of describing AGWA absorption and stripping in ample ranges of (T, P) and composition, using aqueous solvents containing, individually or blended, MEA, MDEA, DEA, and AMP. The approach is based on the Chemical Theory,38 where Complex species are formed in the liquid via chemical reactions with the Real AGWA species. Complex species are reversibly created in absorption and destroyed for solvent regeneration. This RVLE model uses a set of ChE molecular reactions. The insufficient thermodynamic knowledge to handle the high number of AGWA ions is circumvented by such hypothetic, nonionic, and nonvolatile Complexes. The nonionic assumption greatly facilitates thermodynamic calculations and is justifiable, since AGWA electrolytes are weak and created by incomplete dissociations. The RVLE model includes mass balances coupled to VLE for Real species and ChE for Complex formation, with the advantage that model tuning requires exclusively VLE AGWA data, which are more available in the Literature than their nonequilibrium counterparts. Furthermore, since AGWA absorption operates in wide pressure ranges, the properties of both phasesenthalpy and fugacityare estimated via the SRK/PR EOS with classical mixing rules,38 discarding idealities such as ideal solution, ideal gas, and Henry’s Law fugacities. And now the most important: Due to the inexorable consistency of thermodynamics, if such an RVLE model is calibrated using only ordinary VLE AGWA data, then the obscure thermal properties and heat effects must be accessible by means of standard thermodynamic constructions, which make totally unnecessary ad hoc enthalpy correlations. The rest of this work is organized in the following manner: section 2 describes the VLE AGWA database that is used for calibration of the RVLE AGWA model, whose formalism is discussed in section 3. The estimation of parameters of this model is the subject of section 4, which presents the Maximum Likelihood strategy for conducting Implicit Parameter Estimation jointly with full statistical treatment and statistical appreciation of results. The estimated ChE constants are correlated with

involved are probably only chemical conjectures or were described by means of chemical analogies28 and/or were characterized with measured properties only under diluted or very idealized conditions.28 The consideration of diluted or ideal solutions in rate-based approaches prevails everywhere, independently of the assigned method for calculation of rates and fluxes and also independently of the process conditions. As a typical example, consider a recent rate-based work17 prescribing ionic species and nonequilibrium reactions for packing tower absorption of a NG with 3% CO2 and 1% H2S with aqueous 50% w/w MDEA. The SRK EOS is used in the vapor phase whereas ElectrolyteNRTL handles the liquid phase. Interfacial Maxwell−Stefan fluxes are calculated with reaction rate enhancing factors and packed tower correlations for binary mass transfer coefficients. Despite this so heavy modeling machinery, the interfacial VLE is handled by ordinary Henry’s Law. A set of six ChE reactions, with at least two redundant equations, with some ions above, is chosen without further information. Thus, it is conceivable that they rely on simple mass action concentration-based formulas, developed elsewhere32 for dilute isothermal systems. As nonequilibrium reaction, only the liquid phase reaction CO2 + H2O + MDEA → MDEA·H+ + HCO3− is selected, but it is said17 that its rate data varies widely in the literature sources. The theory of Rinker et al.28 is chosen, which is also rate-based and was calibrated with data on diluted systems again and subatmospheric CO2. In other words, both approaches17,28 adopt reaction rate in terms of concentrations and/or partial pressures with Henry’s law in all VLE instances. The gas from the absorber17 loses less than 10% of its fed CO2. Only H2S was abated appreciably, changing the solvent temperature less than 1.5 °C, which is a very inexpressive thermal effect. In AGWA systems, the liquid phase is primarily molecular, dominated by water and amine.32 The existent ions should strongly influence thermodynamic and transport liquid properties. However, such influence, particularly conjugated with strong solvent−ion and solvent−solvent interactions (e.g., solvation and H-bonding), is not well understood. Another fact is that, despite the high concentrations, loadings, pressures, and reaction rates that prevail in industrial AGWA operations, rate-based models usually propose ideal behavior for some species in the aqueous and/or gas phases.14,16,17,19,20,27,28,32 It is also not particularly difficult to encounter high pressure applications (e.g., NG purification) with acid gas partial pressures in the chemical reaction terms.27 This text is not a detraction of rate-based AGWA models. Probably, the essence of Nature is rate-based and equilibrium would be an abstraction that only exists in a limited human perception of Nature. Nevertheless, the rate-based AGWA formalisms, even if appropriately funded, correctly installed within the unit operation, and soundly solved, are still complex and potentially unpractical. Furthermore, they require a wide range of data with lots of physical, chemical, and geometrical parameters, such as equilibrium ionic/nonionic constants, nonequilibrium kinetic constants, interfacial heat/mass transfer coefficients, thermodynamic parameters, thermal data of species, and contact and equipment data.19,20,24−27 These points make it difficult to implement rigorous rate-based models in routine process simulation and optimization. Besides, in opposition to equilibrium AGWA systems whose experimental data are easily found in the Literature, there is a serious lack of reliable nonequilibrium AGWA data for calibration purposes. 9204

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

temperature and submitted in section 5 to a process of conversion of standard states, allowing access to ideal gas properties of formation of Complexes, which are necessary objects for predicting heat effects and temperature profiles in AGWA separations.

2. EQUILIBRIUM DATABASE FOR H2S AND CO2 WITH AQUEOUS ALKANOLAMINES RVLE AGWA formalisms have parametersChE constants requiring calibration with VLE AGWA data in appropriate ranges of (T, P) and alkanolamine composition. VLE AGWA systems have been extensively studied in the last decades, allowing data to be easily found.2 Barbosa et al.39 collected VLE data40−46 up to 2004 for AGWA systems with CO2/H2S/H2O/ MEA/MDEA consolidating a database with 630 VLE runs. Each run reports temperature (°C), solvent composition, CO2/ H2S partial pressures (bar), and CO2/H2S loadings (mol/mol), with the following respective ranges: 40−180 °C, 0−5 M, 0−40 bar, and 0−1.2 mol/mol. In this work, the database of Barbosa et al.39 is extended with data from new sources,29−32,47−55 including: (i) new MEA/MDEA runs; (ii) hundreds of DEA runs; and (iii) almost 200 runs with AMP. This extended databasehenceforth referred to as the Databasecontains 1331 AGWA VLE runs with CO2/H2S/H2O/MEA/MDEA/ DEA/AMP. The majority of runs involve only one acid gas and one alkanolamine in water, but there is a high number of runs assigning one acid gas to two or three alkanolamines or both gases to one alkanolamine. About 300 runs assign both CO2 and H2S to two or three alkanolamines. There are no runs with four alkanolamines. MDEA dominates the Database in terms of number of runs (as well as in terms of %w/w or %mol solvent composition), seconded by DEA and MEA. AMP has the weakest presence, and no run uses only AMP as amine. CO2 dominates as acid gas, with more than 900 runs. The new ranges of temperature and pressure in the Database are broader than the previous ones: 24−180 °C for temperature and 0.001−40 bar for partial pressures. Data at high pressures are useful for studying new absorbers (e.g., NG purification), whereas data at high temperatures are useful for investigating new more resistant solvents under new regeneration schemes at higher pressures. Figure 1 allows a direct appreciation of the Database in terms of three images: (i) histogram of temperatures (Figure 1a) with a large presence of 40 °C, seconded by 100 °C; (ii) histogram of %w/w amine in solvents (Figure 1b), displaying MDEA densely populating the Database from 5%w/w to 50%w/w, whereas AMP is only a marginal cosolvent in a few runs from 3%w/w to 10%w/w; and (iii) plot of loadings versus partial pressures, addressing three populated temperatures (40, 100, and 120 °C) with different symbols for each acid gas (Figure 1c). Figure 1c reveals that, depending on temperature, solvent composition, and acid gas, the entire range of loadings can be swept by varying partial pressure; that is, it strongly suggests linear (parallel) correlations between loading and the logarithm of partial pressure. Parts a−d of Figure 2 depict four “tomographic” views of Figure 1c, each focusing on one amine (MEA, MDEA, DEA, and AMP, respectively), with CO2 plotted as squares and H2S as diamonds, coding runs at 40, 100, and 120 °C with red, blue, and green colors, with the remaining in black. Such images suggest interrelationships between output variable loading (mol/mol) and input variables temperature (K), logarithm of solute partial pressure (bar), and solvent composition (%w/w

Figure 1. AGWA Database: (a) histogram of temperatures (°C); (b) histogram of %w/w amine in solvents; (c) CO2/H2S loadings (mol/ mol) vs partial pressures (bar).

of amine), supporting the identification of short-cut predictors of loadings. Figure 2e depicts a histogram of temperatures in terms of amines and solutes, complemented by Figure 2f, with a histogram of temperatures in terms of chemical reactions. As 9205

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Figure 2. AGWA Database: CO2/H2S loadings (mol/mol) vs partial pressures (bar) with (a) MEA, (b) MDEA, (c) DEA, and (d) AMP; histograms of (e) temperature vs amine, (f) temperature vs chemical reaction, (g) loadings (mol/mol), and (h) partial pressures (bar).

seen in Figure 2e and f, CO2−MDEA is the most frequent pair, with runs at all temperatures, whereas MDEA−H2S runs are only present between 40 and 120 °C. The reactive pair CO2− MDEA has the best coverage (751 occurrences), seconded by H2S−MDEA (384 occurrences). The pair with poorest coverage is H2S−AMP (75 occurrences). Since the Database

contains two acid gases and four amines, eight chemical reactions are possible, but Figure 2f shows that all eight reactions are only seen in the 40 °C extract of runs. Temperatures of 70, 100, and 120 °C are the next in data size, covering seven reactions. However, data at 70, 75, and 80 °C can be regrouped as a new larger 75 °C level, covering all 9206

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Figure 3. Selected data for estimation of ChE constants: histograms of chemical reactions at (a) 40 °C, (c) 75 °C, (e) 100 °C, and (g) 120 °C; distribution of loadings vs partial pressures at (b) 40 °C, (d) 75 °C, (f) 100 °C, and (h) 120 °C.

2.1. Isothermal Data Extracts for Estimation of ChE Constants. Since a true ChE constant is a function of temperature only, ChE constants will be estimated at fixed temperatures, so that the Database should be screened to list isothermal slices of data with sufficient information for estimation. Figure 1a shows the Database with 13 temperatures,

eight reactions. Figure 2g and h display histograms of acid gas partial pressures and loadings, denoting that low pressures and loadings are more frequent for both CO2 and H2S. Probably due to its lower reactivity and wider applications, CO2 runs explore the entire range of partial pressures (0−40 bar), whereas H2S partial pressures are confined from 0 to 15 bar. 9207

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

AMP and CO2−DEA (10 occurrences each) share the worst. Figure 3f presents the distribution of loadings and partial pressures in this set, with an adequate range of loadings and a shrunken range of partial pressure toward high values. (iv) The temperature level of 120 °C, with 84 runs, is chosen in fourth place. As seen in Figure 2f, there are no data for H2S− MEA. Nevertheless, no better choice is available. Hence, all 84 runs at 120 °C are selected, creating a NLP problem for the seven qualified ChE constants with about 2400 variables to estimate and 1600 constraints to satisfy. Figure 3g shows the histogram of reactions for this selection, confirming H2S−MEA as invisible and the remaining seven reactions exposed: CO2− MDEA (57 occurrences) has the best coverage, whereas CO2− MEA (12 occurrences) has the worst. Figure 3h presents the distribution of loadings and partial pressures in this set, with a good range of loadings and a yet more shrunken range of partial pressure toward high values. Table 1 summarizes the final choices of four isothermal extracts of selected runs for parameter estimation. Apparently,

but only a few offer sufficient data to support parameter estimation. Even in the case when a large set of isothermal data is encountered, there would be no assurance that all eight chemical reactions of interest (CO2−MEA, CO2−MDEA, CO2−DEA, CO2−AMP, H2S−MEA, H2S−MDEA, H2S− DEA, H2S−AMP) would be adequately exposed. For example, in Figure 2f, several temperatures do not encompass all eight reactions; only 40 °C has sufficient data to estimate all eight ChE constants. Thus, the strategy is to seek an abridged list of temperatures, each one with a subset of isothermal runs as large as possible, in order to guarantee the following: (i) statistical coverage of all eight reactions in all isothermal extracts; and (ii) sufficient number of different isothermal extracts allowing thermodynamic processing of the estimated ChE constants for further inferences of Complex properties (e.g., ideal gas heat capacity and standard enthalpy of formation), which means that at least four estimated values of a ChE constant at four temperatures must be available for regression of threeparameter formulas such as ln K = A + B/T + C ln T. As shown in section 3, each run involves 30−37 observed variables and 20−28 nonlinear constraints to be satisfied by estimated values. For instance, when one selects 100 isothermal runs for parameter estimation, the numerical burden corresponds to more than 3000 variables and 2000 constraints. Hence, the size of the isothermal set of runs must be considered to avoid overburden. It is impossible to identify four temperatures with adequate coverage of all eight reactions in this Database. Nevertheless, reliable results are possible for at least six reactions and partial results for the remaining two reactions, according to the following ad hoc procedure: (i) The largest temperature level of 40 °C, with 536 runs (Figure 1a), is first chosen. In order to limit the numerical burden and distribute loadings, partial pressures, and reactions, a subset with 320 runs is selected at 40 °C. This gives a NLP problem for eight ChE constants with about 9000 variables to estimate and 6000 constraints to satisfy. Figure 3a shows the histogram of chemical reactions for this selection, confirming all eight reactions exposed: reaction CO2−MDEA (164 occurrences) has the best coverage, whereas H2S−AMP (26 occurrences) has the poorest. Figure 3b presents the distribution of loadings and partial pressures in this set, confirming the respective ranges as entirely swept. (ii) The temperature levels of 70, 75, and 80 °C are merged in an extended 75 °C level with 244 runs. A price is paid in terms of loss of accuracy, but this regrouping now covers all eight reactions. A subset with 171 runs is selected in this 75 °C level, invoking an NLP problem for eight ChE constants with about 5000 variables to estimate and 3000 constraints to satisfy. Figure 3c shows the histogram of reactions for this selection, confirming all eight reactions as exposed: CO2−MDEA (79 occurrences) has the best coverage, whereas CO2−AMP (22 occurrences) has the worst. Figure 3d presents the distribution of loadings and partial pressures in this set, evidencing the respective ranges as swept. (iii) The temperature level of 100 °C, with 193 runs, is chosen. As seen in Figure 2f, there are no data for the H2S− AMP reaction, but no better option is available. Hence, a subset with 114 runs is selected at 100 °C, creating an NLP problem for the seven qualified ChE constants with about 3000 variables to estimate and 2000 constraints to satisfy. Figure 3e shows the histogram of reactions for this selection, confirming H2S−AMP as invisible and the remaining seven reactions exposed: CO2− MDEA (56 occurrences) has the best coverage, whereas CO2−

Table 1. Selected VLE AGWA Runs for Estimation of ChE Constants isothermal extracts of data 40 °C selected runs

320

CO2 loading (mol/mol) H2S loading (mol/mol) CO2 partial pressure (bar) H2S partial pressure (bar) ChE constants for estimation “invisible” ChE constant

[0,1.3] [0,1.5] [0,32] [0,23] 8

75 °C 171 Ranges of [0,1.1] [0,1.2] [0,38] [0,10] 8

100 °C 114 Data Values [0,0.9] [0,0.9] [0,39] [0,3.8] 7 KH2S‑AMP

120 °C 84 [0,0.7] [0,0.8] [0,38] [0,11] 7 KH2S‑MEA

there is sufficiency of data in the sets of 40 and 75 °C for all eight ChE constants. The H2S−AMP and H2S−MEA ChE constants cannot be estimated, respectively, at 100 and 120 °C. However, this is only the tip of the iceberg.

3. EQUILIBRIUM MODEL FOR AGWA SYSTEMS The RVLE approach to AGWA systems prescribes the superposition of VLE onto ChE for a complete set of chemical reactions. The intent is to develop a framework, precise enough for engineering purposes, reproducing reversible equilibrium changes in operations, easily tunable, and installable in process simulators. It should be emphasized that not only the absorption/stripping composition changes are of interest but also the associated thermal effects, as economic and environmental evaluation of AGWA technologies relies on heat consumption. Since thermodynamics does not depend on a precise knowledge of the internal structure of phases, the approach is molecular, with fictitious species in the liquid phasethe Complexesreplacing the actual weak AGWA ions. Complexes are nonvolatile, not participating in the VLE, but influencing the liquid phase properties. This is manifested via the ChE constants, which are estimated with VLE AGWA data. The molecular nature of both phases allows them to be described via conventional SRK/PR EOS supporting applications up to high pressures and all compositions. The present approach is inherently simpler than any rate-based counterpart; nevertheless, it relies on the solid foundations of equilibrium; that is, it is asymptotically exact at limit conditions of AGWA operations and has the advantage that its calibration only 9208

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

constraints, and its variables are removed from the set of variables to be estimated. (vi) All Complex forming reactions in eq 1 are independent, since the atomic matrix of Real species has rank five in the space of atoms (H, C, O, N, S). Thus, eq 1 is a complete set for chemical changes. The stoichiometric matrix (n x nC) of Real species [CO2 H2S H2O MEA MDEA DEA AMP] is given by the following:

demands VLE AGWA data. The formalism adopts the Chemical Theory,38 which is applicable in complex phase equilibrium systems (e.g., with strong H-bonds, solvation, etc). Chemical Theory prescribes a set of independent reactions in ChE, converting Real species into Complex species. Such Complex species are not actually in the system but can mimic the observed effects via the composition changes induced by their formation. Tuning this framework to data allows the apparent thermodynamic behavior in question to be reproduced. Chemical Theory is implemented with the following premises: (i) Two kinds of species: Real and Complex. Real species include solutes CO2 and H2S, amines (MEA, MDEA, DEA, AMP), and water. Water has a central role in all reactions forming Complexes; without it, there would be no reversible affinity of acid gases for the liquid. (ii) Each Complex forming reaction has as reactants an acid gas, water, and an amine. Reactions are defined with stoichiometric coefficient 1 for all reactants. Despite that water could react with a proportion higher than 1 mol per mol of Complex, the coefficient 1 is kept because it is simpler and sufficient for all that matters. A possible issue is that unitary amine coefficients make it difficult to produce loadings well above 1 mol/mol (Figure 1c evidences loadings of 1.1−1.3 mol/mol at high pressures). Nevertheless, the unitary choice is maintained, as it does not even conflict with the aphorism that a limiting value of 1/2 should exist for CO2 loading in MEA. The reason lies on ChE: despite the unitary coefficients, ChE can mimic the loading barrier, if necessary, by blocking the reaction extension near 1/2 in order to match data. (iii) The eight Complexes CO2−H2O−MEA, CO2−H2O− MDEA, CO2−H2O−DEA, CO2−H2O−AMP, H2S−H2O− MEA, H2S−H2O−MDEA, H2S−H2O−DEA, and H2S−H2O− AMP are created by eight chemical reactions in eq 1 from four amines, water, and two acid gases. Equation 1 displays the eight ChE constants and the chosen standard states (section 3.1) as Pure Saturated Liquid (Liq) at T for Complexes and Pure Ideal Gas (g) at T and 1 bar for Real species:

⎡− 1 ⎢ ⎢0 ⎢− 1 Π = ⎢⎢−1 ⎢0 ⎢0 ⎢ ⎣0

−1 0 −1 0 −1 0 0

−1 0 −1 0 0 −1 0

−1 0 −1 0 0 0 −1

0 −1 −1 −1 0 0 0

0 −1 −1 0 −1 0 0

0 −1 −1 0 0 −1 0

0⎤ ⎥ −1⎥ −1⎥ ⎥ 0⎥ 0⎥ 0 ⎥⎥ −1⎦

(2)

(vii) Both phases are molecular without ions. The thermodynamic properties of vapor and liquid are calculated via the SRK/PR EOS. Thus, the critical constants (TC, PC, ω) and normal boiling point (TB) of Complexes are therefore necessary. These values, per se, do not have to make sense or be accurate, but they must imply negligible vapor pressures as viewed by the EOS. These constants are estimated in Table 2 Table 2. Constants of Complexes via the Joback Method56 complex

MW

TC (K)

PC (bar)

ω

TB (K)

CO2−H2O−MEA CO2−H2O−MDEA CO2−H2O−DEA CO2−H2O−AMP H2S−H2O−MEA H2S−H2O−MDEA H2S−H2O−DEA H2S−H2O−AMP

123.11 181.19 167.16 151.16 113.18 171.26 157.23 141.23

845.8 880.4 903.0 874.0 758.4 794.1 819.9 788.3

61.24 42.58 49.29 48.98 60.18 42.21 48.68 46.09

1.699 1.991 2.103 1.643 1.154 1.630 1.723 1.192

659.3 714.2 732.4 686.4 560.6 628.3 648.3 595.3

via the Joback Contribution Method56 assuming “plausible structures” for the Complexes in eq 1. Other fundamental properties of Complexes, such as thermal data of formation, are not estimated by contribution methods. A sound thermodynamic approach for obtaining them is discussed after the estimation of ChE constants. At the present point, the constants in Table 2 are the only values needed to work with the EOS. Notice the high TB’s, confirming low vapor pressures of Complexes. 3.1. Thermodynamic Modeling of AGWA VLE Experiments. Model calibration demands a thermodynamic description of AGWA VLE experiments. A “theorized” experiment in Figure 4 consists in charging numbers of moles (N) of Real species CO2, H2S, H2O, MEA, MDEA, DEA, and AMP in a closed system under a specified temperature (T). At equilibrium, the system is characterized by the following: pressure, P; mole fractions of Real species in liquid and vapor phases (X, Y); mole fractions of Complexes in the liquid phase (XC); molar content of phases (L,V); solute partial pressures (PS); and solute loadings (αS). Subscripts A, S, C and superscripts V, L refer to Amines, Solutes, Complexes, Vapor, and Liquid phases; n, nC, nA, nS are numbers of Real, Complex, Amine, and Solute species. The equations imposed on experimental variables are the experimental constraints in Table 3: eq 3 for mass balances of Real species (n eqs); eqs 4

K1

CO2 (g) + H 2O(g) + MEA(g) ↔ CO2 −H 2O−MEA(Liq) CO2 (g) + H 2O(g) + MDEA(g) K2

↔ CO2 −H 2O−MDEA(Liq) K3

CO2 (g) + H 2O(g) + DEA(g) ↔ CO2 −H 2O−DEA(Liq) K4

CO2 (g) + H 2O(g) + AMP(g) ↔ CO2 −H 2O−AMP(Liq) K5

H 2S(g) + H 2O(g) + MEA(g ) ↔ H 2S−H 2O−MEA(Liq) H 2S(g) + H 2O(g) + MDEA(g) K6

↔ H 2S−H 2O−MDEA(Liq) K7

H 2S(g) + H 2O(g) + DEA(g) ↔ H 2S−H 2O−DEA(Liq) K8

H 2S(g) + H 2O(g) + AMP(g) ↔ H 2S−H 2O−A MP(Liq) (1)

(iv) Only Real species are in the VLE; Complexes have null vapor pressures and remain liquid. (v) All species participate in the ChE expressed by the eight ChE constraints in eq 1. If a reaction is absent (due to an absent reactant), its ChE is withdrawn from the set of 9209

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

μ̲ = μ̲ °(T ) + RT ln a ̂,

μ̲ C = μ̲ C° (T ) + RT ln aĈ (11)

ΔG°(T ) = μ̲ C° (T ) + ΠT μ̲ ° (T ) = − RT ln K̲ (T ) L

(12)

L

a ̂ = f ̂ ·/ f̲ ° = f ̂ = φ̂ L ·PX̲

(13)

L aĈ = f ̂ ·/ f̲ ° (T , P̲ CSAT(T )) C

C

= X̲ C · γ̲ C · f̲ Figure 4. VLE experiment with the AGWA system.

constraints [label]

expression

N̲ + L ΠX̲ C − LX̲ − VY̲ = 0̲

∑ Xj + ∑ XCk − 1 = 0

2 Normalizations [SX]

(3)

nC

n j=1

k=1

(4a)

n

∑ Yj − 1 = 0 nS Solute Partial Pressures [PP]

P̲ S − P S S Y̲ = 0̲

(5)

T

nS Solute Loadings [LD]

( 1̲ S AN̲ ) α̲ S − S S(N̲ − VY̲ ) = 0̲

n VLE of Real Species [LV]

V L ln f ̂ − ln f ̂ = 0̲ L

(6)

(7)

ln X̲ C + ΠT ln f ̂ − ln K̲ (T ) = 0̲

nC Reactions ChE [CE]

(8)

for normalization of liquid/vapor mole fractions (2 eqs); eq 5 for solute partial pressures (nS eqs); eq 6 for solute loadings (nS eqs); eq 7 for VLE of Real species (n eqs); and eq 8 for ChE of Complex formations (nC eqs). Vectors f ̂L, f ̂V, ϕ L̂ , and ϕ V̂ represent fugacities and fugacity coefficients of Real species in liquid and vapor phases in eq 9 (· applies position-to-position algebra). Fugacity coefficients are calculated with the SRK or PR EOS. ChE constraints use the classical eq 10, where μ̲ is the vector of chemical potentials of Real species with identical values in both phases, related to activities (a )̂ and standard state potentials (μ̲°(T)) via eq 11a. Activities of Real species have identical values in both phases. Complex activities in eq 11b are restricted to the liquid phase. The vector (nC by 1) of ChE constants (K̲ (T)) is written from the vector of the standard Gibbs free energy of reactions (ΔG°(T)) as in eq 12. The simplest assignment of standard state to Real species is Pure Ideal Gas at T and P = 1 bar ((g) state) as in eq 13 with f°̲ = 1 bar. Due to their high stability as liquids, the choice of (g) state for Complexes leads to numerically ill-conditioned EOS liquid phase fugacities (≈10−23), preventing convergence of the estimation algorithm. To avoid this, standard states of Complexes are chosen as Pure Saturated Liquids at T ((Liq) state) as in eqs 14 and 1. The RHS of eq 14 is simplified with the following: (i) small molar volumes of Complexes cancel pure liquid fugacities; (ii) unitary liquid phase activity coefficients for Complexes (γC ≅ 1). With eqs 11−14, the ChE constraints in eq 10 can be put as in eq 8 in Table 3. V

f ̂ = PY̲ · φ̂V ,

μ̲ C + Πt · μ̲ = 0̲

L

f ̂ = PX̲ · φ̂ L

L C

(T , P̲ CSAT(T )) = X̲ C

(14)

Z̲ T = [T P̲ ST α̲ ST N̲ T X̲ T Y̲ T X̲ CT P V L ]

(15a)

F̲ (Z̲ , θ̲ ) = 0̲

(15b)

For the CO2, H2S, H2O, MEA, MDEA, DEA, and AMP systems, the maximum number of experimental variables is nZ = 3n + 4 + nC + 2nS = 37 and the maximum number of experiment constraints is nF = 2n + 2 + nr + 2nS = 28. Thus, at a fixed temperature level, the problem consists in estimating nθ = nC = 8 parameters using a set of runs, each one with a maximum of 37 variables and 28 nonlinear constraints to be satisfied by estimated variables. 3.2. Statistical Assumptions for Parameter Estimation in AGWA Systems. In Table 1, four extracts of isothermal runs are selected from the AGWA Database for estimating ChE constants of the RVLE model at four temperatures. Consider the isothermal extract at temperature T with nE runs. Sizes nZ and nF in section 3.1 are limiting values; that is, a run may use only a fraction of the set of Real species activating just a subset of constraints. Define, for the ith run, nZi and nFi as its numbers of Active measured variables and Active constraints (nZi ≤ nZ, nFi ≤ nF). The constraint vector for run i is F̲i = 0̲ with size nFi × 1. Experimental Variables have Observed, Correct, and Estimated Values. A symbolic discrimination among such classes of active variables for the ith run is introduced as vectors (all nZi × 1) of: (i) correct experiment variables (Z̲ i); (ii) estimated experiment variables ( Ẑ i); and (iii) observed experimental variables (Z̲ Ei ). Analogously, two vectors of model parameters are defined (all nθ × 1, nθ = nC): (i) correct model parameters (θ̲); and (ii) ̂ Correct values are supposed estimated model parameters ( θ ). to exist but are unknown; only observed and estimated values are available. At this point, praxis invokes the following assumptions: Assumption I: The model in eqs 15 is correct; Assumption II: Observed values (Z̲ Ei ) for run i are uncorrelated and normally distributed around correct values (E(Z̲ Ei ) = Z̲ i), with variance−covariance matrix in eq 17 with a known positive definite diagonal matrix W i and σ2E an unknown Fundamental Variance; Assumption III: Experiments are independent. Assumption I implies that correct and estimated values of parameters and variables, respectively, satisfy eqs 16a and 16b. Assumption II implies that σ2E can be estimated a posteriori. Equation 17a is a simple model with a single unknown (σ2E). Other variance models can be conceived57 but are beyond the scope. With Assumption II, the probability density function (PDF) of observed Z̲ Ei is given by eq 17b, giving E(Z̲ Ei ) = Z̲ i and COV (Z̲ Ei ) = σ2E W −1 i . Assumption III implies that the PDF of the

(4b)

j=1

(T , P ) · / f ̲

Experimental constraints and variables (eqs 3−8) are written in eq 15 after grouping variables as Z. Equation 15b is the set of experiment constraints. Model parameters (θ̲ = K̲ ) comprise ChE constants for all Complex reactions at the selected temperature, i.e., nθ = nC = 8 parameters per temperature level.

Table 3. Set of Constraints in AGWA VLE Experiments n Real Species Mass Balance [MB]

L C

(9) (10) 9210

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

set of nE isothermal runs is given by the product of all forms (eq 17b) as in eq 17c. F̲ i(Z̲ i , θ̲ ) = 0̲

(i = 1, ..., nE)

(16a)

Fi( Ẑ i , θ )̂ = 0̲

(i = 1, ..., nE)

(16b)

COV (Z̲ iE) = σE2 W i−1

(i = 1, ..., nE)

them. The estimation algorithm effectively establishes a hierarchy of search priorities, primordially persecuting the TOVs, secondarily the best UOVs, and, with least priority, the worst UOVs. The completion of the data set with UOVs is important because, as said in statistics, an independent coarse estimate is always better than no estimate. Albeit untrue, UOVs can be sensibly guessed: pressure (P) is guessed as 110% of the sum of partial pressures of solutes (P̲S); molar hold up of liquid (L) is guessed as 90% of the loaded solvent (NMEA + NMDEA + NDEA + NAMP + NH2O); preparation mole numbers of solutes are guessed as twice the loading times of the total of loaded amine (N̲ S = 2α̲ S(NMEA + NMDEA + NDEA + NAMP)), and so on. Nonetheless, more important than guessing UOVs is the assignment of the respective “experimental” standard deviations, which must always be conservative. The discrimination between TOVs and UOVs of run i is coded into matrix W i, whose diagonal elements are inverses of experimental variances of the observed variables (eq 17a). Small standard deviations (0.1% to 5% of the observed value, lower-bounded absolutely by 10−3) are assigned to TOVs giving large (heavy) positive diagonal positions in W i. UOVs have large standard deviations (order of 10% to 200% of the observed value, absolutely lowerbounded by 1), giving small (light) positive diagonal positions in W i. During parameter estimation, the “pressure” for adherence of predictions ( Ẑ i) onto observed data (Z̲ Ei ) is dosed by W i: its “heavier” TOV positions orient the search to adhere on TOVs much more intensely than UOVs. In other words, since parameter estimation always reduces to an optimization search, good adherence of predictions onto TOVs implies higher gains in the objective function than a similar adherence onto UOVs. For this reason, when presenting results of parameter estimation against AGWA data, the focus is primarily placed on the adherence of predictions onto TOVs such as partial pressures and loadings.

(17a)

⎡ ⎛ 1 PDF(Z̲ iE) = ⎢ | W i| exp⎜ − 2 (Z̲ iE − Z̲ i)T ⎢⎣ ⎝ 2σE ⎞⎤ × W i(Z̲ iE − Z̲ i)⎟⎥ ⎠⎥⎦

⎤ ⎡ ⎢(2πσE2)nZi /2 ⎥ ⎦⎥ ⎣⎢

(i = 1, ..., nE) (17b)

⎡ E PDF(Z̲ 1E , ..., Z̲ nE )=⎢ ⎢⎣ ⎛ −1 × exp⎜⎜ 2 ⎝ 2σE

nE

∏ | W i| i

nE

⎞⎤

i

⎠⎦

∑ (Z̲ iE − Z̲ i)T W i(Z̲ iE − Z̲ i)⎟⎟⎥⎥

⎡ ⎤ ⎢(2πσ 2)(∑inE nZi)/2 ⎥ E ⎢⎣ ⎦⎥

(17c)

3.3. True and Untrue Observed Values of Experimental Variables. A VLE experiment i typically has only a subset of Z̲ Ei accessible in the Literature, i.e., temperature (T), solute partial pressures (P̲S), solute loadings (α̲S), and solvent composition (NMEA, NMDEA, NDEA, NAMP, NH2O). Since Z̲ Ei in section 3.1 contains several other variables (mole fractions, phase mole numbers, pressure, etc.), two classes of Observed Values were devised: True Observed Values (TOV) and Untrue Observed Values (UOV). TOVs are actually measured values in the AGWA Literature, while UOVs are crude guesses for the remaining variables. There are also TOVs not actually measured, which can be inferred without doubt (e.g., zero mole fraction of Real species inexistent in a run). Since a statistically sound estimation procedure is used, the discrimination between TOVs and UOVs is passed to the algorithm through the respective standard deviations: TOVs have small experimental standard deviations ranging from 0.1% to 5%, whereas UOVs (as crude guesses produced with TOVs) have larger “experimental” standard deviations ranging from 10% to 200% (or more). Guessing UOVs is not particularly difficult. For example, if an experiment at 40 °C reports partial pressures of H2S and CO2 with standard deviations of 5%, the pressure can be guessed as the sum of partial pressures increased by 1% to 10% to account for water and MEA; that is, this UOV uncertainty implies an equivalent “experimental” standard deviation of, say, 20%. In this work, the methodology of parameter estimation requires the entire vector Z̲ iE as numerically available; that is, its TOV portion should be complemented with UOVs for the rest of the variables. The insertion of UOVs is done exclusively to allow implicit parameter estimation based on the implementation of the maximum likelihood criterion (MLC) known as the method of estimated deviations,57 which requires that all variables to be estimated must have observed counterparts. The insertion of UOVs does not hinder or distort the statistical estimation because it is also true that large uncertainties are attached to

4. MLC PARAMETER ESTIMATION OF CHE CONSTANTS Consider the extract of nE isothermal runs from the AGWA database in section 3.2. The goal is to estimate 8 ChE constants for this set. However, AGWA experiments are complex with RVLE of highly nonideal liquids and nonideal gases. Also, due to variable assignments of species in AGWA experiments, runs may vary widely in terms of number of active variables and constraints, and in general, a given run is affected by only part of the set of ChE constants under estimation. The set of constraints of run i is shown in eq 16b, and the set of experiment variables is shown in eq 15a. The number of variables (≈37) is always greater than the number of constraints (≈28), all supposed to be qualified. Due to the highly nonlinear constraints (VLE and ChE), it is impossible to explicit part of the variables in terms of the rest; that is, it is unfeasible to separate variables into the traditional dichotomy of input and output variables, where input variables are not subjected to experimental errors, in opposition to output variables that bear the incidence of all errors (model + experiment). It turns out that the present parameter estimation is nonlinear and implicit, but it can be well-posed, depending on the selected data (section 2.1) and on the appropriateness of the model (section 3.1). Additionally, the problem is large in size, especially in the 40 °C extract, with 320 selected runs. To estimate the 8 ChE constants at 40 °C, about 9000 variables have to be estimated, and 6000 nonlinear constraints have to be 9211

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

satisfied. At least two traditional NLP solvers were tried in this 40 °C case, and both have failed. The approach herein proposed is an implementation of the MLC for implicit parameter estimation known as the method of estimated deviations57 (MED). The MED generates MLC estimates for model parameters ( θ )̂ and all variables ( Ẑ i (i = 1, ..., nE)), whereas the method of observed deviations (MOD), also known as Deming’s method,57 is a simpler strategy that does not estimate variables but can be used to generate initial parameter estimates to MED. The MOD was not used in this work. Instead, MED was implemented in stages, as a succession of progressively harder MED applications, such that good initial estimates for a MED in this succession are furnished by the converged values in the previous one. This implementation was successful and able to handle AGWA NLPs with huge numbers of variables and constraints. 4.1. MLC for Implicit Parameter Estimation. For the extract of nE isothermal runs, Assumptions II and III give the joint PDF of all measured variables in eq 17c. However, since the correct values Z̲ i are not accessible, a likelihood function, LF, is built in eq 18, which mimics the true PDF in eq 17c, replacing correct values by estimated ones: E , Z1̂ , ..., Ẑ nE) LF(Z̲ 1E , ..., Z̲ nE ⎡ nE ⎛ 1 = ⎢ ∏ | W i| exp⎜⎜ − 2 ⎢⎣ ⎝ 2σE i

⎞⎤ × W i(Z̲ iE − Ẑ i)⎟⎟⎥ ⎠⎦⎥

W i( Ẑ i − Z̲ iE) + H iT λ̲ i = 0̲ nE

∑ J iT λ̲ i = 0̲ i

F̲ i( Ẑ i , θ )̂ = 0̲

i

(21)

(i = 1, ..., nE) (i = 1, ..., nE)

(m)

(m)

F̲ i( Ẑ i , θ )̂ ≅ F̲ i(m) + J ( θ ̂ − θ ̂ ) + H i( Ẑ i − Ẑ i )

nE

i

E

∑ (Z̲ i − Ẑi)T

(i = 1, ..., nE)

i

⎡ ⎤ ⎢(2πσ 2)(∑inE nZi)/2 ⎥ E ⎢⎣ ⎥⎦

(23)

(̂ m) where F̲(m) = F̲i( Ẑ (m) ). Vector F̲(m) and matrices H i and J i i i , θ i (m) (̂ m) ̂ are calculated at ( Z i , θ ). Equating the RHS of eq 23 to zero and substituting it into eq 21c allows eqs 21 to be solved for Ẑ i, θ ,̂ λ̲i (i = 1, ..., nE) giving the following:

(18)

(m)

θ̂ = θ̂

+ S̲ θ(m)

(24a)

(m) Ẑ i = Ẑ i + S̲ Zi(m)

(i = 1, ..., nE)

(24b)

(m) λ̲ i = M i(F̲ i(m) + H i(Z̲ iE − Ẑ i ) + J S̲ θ(m)) i

(i = 1, ..., nE)

(24c)

where search directions

{ θ ̂, Z1̂ , ..., Ẑ nE}

S̲(m) θ

and

S̲(m) Zi

are given by eqs 25:

nE

s. t . (i = 1, ..., nE)

(m)

S̲ θ(m) = −[Ξ]−1 {∑ J T M i(F̲ i(m) + H i(Z̲ iE − Ẑ i ))} i

(19)

i

(25a)

This problem is solved with the associate Lagrangian:

S̲ Zi(m) = { I − W i−1 H iT M i H i}(Z̲ iE −

nE

(m) Ẑ i )

M i{F̲ i(m) + J S̲ θ(m)}

Λ = (1/2) ∑ ( Ẑ i − Z̲ iE)T W i( Ẑ i − Z̲ iE)

i

i

− W i−1 H iT (25b)

which incorporate new matrices (sizes, respectively, nFi × nFi and nθ × nθ) defined in eqs 26:

nE

λ̲ iT F̲ i( Ẑ i ,

(22)

It is assumed that all active constraints of run i (Fi) are qualified. Hence, H i and J i have maximum ranks: rank( H i) = nFi, rank( J i) = nθ (i = 1, ..., nE). 4.2. Gauss−Newton Method (GNM) for MLC Implicit Estimation. To solve eqs 21, local searches via the iterative Gauss−Newton method57 (GNM) are invoked. The present implementation is based on the MED strategy, generating simultaneous estimates for all variables ( Ẑ i) and parameters ̂ GNM relaxes constraints (eq 21c) around the estimates ( θ ). (m) ̂ ( Z i , θ (̂ m)) at iteration (m) via eq 23:

i

i

(at Ẑ i , θ )̂

J = (∇̲ θ F̲ iT )T

min Ψ = (1/2) ∑ ( Ẑ i − Z̲ iE)T W i( Ẑ i − Z̲ iE)



(at Ẑ i , θ )̂

H i = (∇̲ Z F̲ iT )T

nE

+

(i = 1, ..., nE)

where jacobian matrices H i, J i (sizes nFi × nZi and nFi × nθ) for the ith run with respect to estimated variables ( Ẑ i) and model parameters ( θ )̂ are defined in eq 22:

The MLC implicit estimation problem57 for parameters ( θ )̂ and variables ( Ẑ i) consists in maximizing LF in eq 18 subjected to the constraints in eq 16b. This entails minimizing the sum of squares in eq 18, subjected to the constraints (eq 16b), as in eq 19:

F̲ i( Ẑ i , θ )̂ = 0̲

(i = 1, ..., nE)

θ )̂ (20)

M i = [ H i W i−1 H iT ]−1

where λi is the vector (nFi × 1) of Lagrange multipliers for run i. The condition of stationary point for eq 20, imposed in terms of Ẑ i , θ ̂, λ̲ i (i = 1, ..., nE), gives eqs 21:

(i = 1, ..., nE)

(26a)

nE

Ξ=

∑ J iT M i J i i

9212

(26b)

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Since rank( H i) = nFi and rank( J i) = nθ (i = 1, ..., nE), the matrix inversion in eq 26a is feasible. Similarly, rank( M i) = nFi and rank(Ξ i) = nθ; that is, M i (nFi × nFi) and Ξ (nθ × nθ) are implies that also invertible. The definite positiveness of W −1 i −1 T Ξ is also H H M −1 M Thus, = and are positive definite. W i i i i i positive definite, as a sum of positive definite matrices. Equations 24−26 are used iteratively until convergence, where the norms of search directions become small. However, this is not globally guaranteed due to the local characteristics of GNM; that is, it can only find a solution θ ,̂ Ẑ i (i = 1, ..., nE) near a good initial estimate ( θ (̂ 0), Ẑ (0) i (i = 1, ..., nE)). For this reason, a damping factor t (t ∈ [0, 1]) multiplies the last term in eqs 24a and 24b as a function of the norms of search directions; that is, large norms entail small t. When GNM approaches convergence, t is released to 1, allowing O(2). Despite these shortcomings, GNM is suited to the present application, because: (i) if adequately programmed (with sparse block handling and sparse linear solvers), GNM can deal with thousands of nonlinear constraints and thousands of variables, within an average PC; and (ii) convergence is fast for wellposed problems (i.e., jacobians sustain maximum rank; ChE constants are identifiable with the data; UOVs and their standard deviations are cautiously guessed for data completion; and initial estimates are good). Appendix A (see Supporting Information) describes the GNM algorithm for solving the implicit parameter estimation in AGWA systems, which was executed successfully to estimate 24 out of 32 ChE constants for 8 chemical reactions at 4 isothermal sets (section 4.5). The GNM in Appendix A can also execute estimation for a subset of parameters, with the remaining parameters kept fixed. This is useful when processing AGWA data with insufficient information to spot some reaction(s) in the original set of 8 reactions. In such cases, GNM can only converge if the deficient ChE constant leaves the target vector and is kept fixed. 4.3. Statistical Analysis of MLC Implicit Parameter Estimation. After convergence of GNM to θ ,̂ Ẑ i (i = 1, ..., nE) for an isothermal extract Z̲ Ei , W i (i = 1, ..., nE) with nE runs, a statistical analysis is performed. Despite the highly nonlinear constraints in Table 3, statistical analysis is based on linear extrapolators, requiring the assumption that the GNM solution θ ,̂ Ẑ i (i = 1, ..., nE) is close to the respective (unknown) correct values; that is: θ ̂ ≅ θ̲ ,

Ẑ i ≅ Z̲ i

(i = 1, ..., nE)

With eqs 16a and 16b, the linearized set of constraints for run i reads as follows: J ( θ ̂ − θ̲ ) + H i( Ẑ i − Z̲ i) = 0̲ i

J = (∇̲ θ F̲ iT )T i

(at Z̲ i , θ̲ ) (at Z̲ i , θ̲ )

W i( Ẑ i − Z̲ iE) + H iT λ̲ i = 0̲ nE i

J ( θ ̂ − θ̲ ) + H i( Ẑ i − Z̲ i) = 0̲

(i = 1, ..., nE)

i

(31)

Equations 31 translate the MLC estimation of θ ,̂ Ẑ i (i = 1, ..., nE) asymptotically near the correct values θ ,̂ Ẑ i (i = 1, ..., nE). These equations are linear with maximum rank in all jacobians, and are solved for Ẑ i, θ ,̂ λ̲i (i = 1, ..., nE) in terms of Z̲ Ei − Z̲ i, giving the following: nE

θ ̂ = θ̲ − [Ξ]−1 {∑ J T M i H i(Z̲ iE − Z̲ i)} i

i

(32)

nE

λ̲ i = − M i J [Ξ]−1 {∑ J T M j H j(Z̲ jE − Z̲ j)} i

j

j

+ M i H i(Z̲ iE − Z̲ i)

(i = 1, ..., nE)

(33)

Ẑ i = Z̲ iE − W i−1 H iT M i H i(Z̲ iE − Z̲ i) nE

+ W i−1 H iT M i J [Ξ]−1 {∑ J T M j H j(Z̲ jE − Z̲ j)} i

j

j

(34)

where square matrices M i (nFi × nFi) and Ξ (nθ × nθ) sustain rank and are invertible (eqs 35): M i = [ H i W i−1 H iT ]−1

(i = 1, ..., nE)

(35a)

nE

Ξ=

(27)

∑ J iT M i J i

(35b)

i

W −1 i ,

M −1 i

T H i W −1 i Hi

From the definite positiveness of = is also positive definite and so is M i . This implies that Ξ is also positive definite. Matrices H i, J i, W i, M i, and Ξ are all functions of correct values; that is, they are true constants. The random variables on the RHS of eqs 32−34 are only the observed Z̲ Ei (i = 1, ..., nE), obeying normal behavior by Assumption II and eq 17b. Equations 32−34 are linear formulas, with constant matrices, in terms of normal random variables Z̲ Ei (i = 1, ..., nE). Thus, the estimates on the LHS of eqs 32−34 are normal variables also. Assumption II and eq 17b give the following:

(28)

Then, with eq 27, one can write the asymptotically correct expression in eq 29:

E(Z̲ iE) = Z̲ i ,

F̲ i( θ ̂, Ẑ i) ≅ F̲ i( θ̲ , Z̲ i) + J ( θ ̂ − θ̲ ) + H i( Ẑ i − Z̲ i)

COV (Z̲ iE) = σE2 W −1

(36)

Additionally, Assumption III implies that

i

(i = 1, ..., nE)

(i = 1, ..., nE)

∑ J iT λ̲ i = 0̲

(i = 1, ..., nE) (i = 1, ..., nE)

(30)

Equation 30 replaces the nonlinear constraints of run i, allowing linear propagation for statistical analysis. Returning to the Lagrangian in eq 20, replacing the model constraints by eq 30, and taking the stationary condition, a set similar to eqs 21 is reached as shown in eqs 31:

By Assumption I, correct values and estimated values must both satisfy model constraints in eqs 16a and 16b. Now let H i and J i represent the jacobian matrices of the ith run, respectively, for the correct experimental variables (Z̲ i) and for the correct parameters (θ̲), both with maximum rank (rank( H i) = nFi, rank( J i) = nθ (i = 1, ..., nE)), and at (Z̲ i, θ̲): H i = (∇̲ Z F̲ iT )T

(i = 1, ..., nE)

E((Z̲ iE − Z̲ i)(Z̲ jE − Z̲ j)T ) = 0

(29) 9213

( i ≠ j)

(37)

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

nonlinear implicit models are involved. More specifically, it seems that the Literature is somewhat confused (or at least, not clear) when an estimator for σ2E (i.e. σ̂2E) is sought in implicit regressions. This statistic is critical for variance−covariance matrices (eqs 41 and 42) and for confidence intervals and domains of correct values. In other words, if σ̂2E is wrongly estimated, several estimators and statistics for validation of results can be compromised and report inappropriate values. In this context, consider the works of Britt and Luecke57 and of Anderson et al.,58 both landmarks in nonlinear MLC parameter estimation that have been cited along decades by hundreds of subsequent authors. The form chosen to explain σ̂2E in both works may have led others to wrong results. The critical item affecting σ̂2E is its number of degrees of freedom (DF). Two kinds of problems can be detected in posterior papers citing Britt and Luecke57 and/or Anderson et al.:58 (i) DF was wrongly reported because one of the two original sources used a wrong DF and propagated such misuse; or (ii) DF in the original source was correct but was wrongly adapted in the citing paper. Britt and Luecke57 considered MLC implicit parameter estimation with a total of q measured variables, k constraints, and n parameters, such that (ibid.) q ≥ k > n. The DF to estimate σ̂2E was chosen, without demonstration, as DF = k − n, justified by only a brief citation of a previous book by Deming.59 Anderson et al.58 considered MLC implicit estimation with N experiments, M measured variables per experiment, 2 constraints per experiment, and L parameters. The DF to estimate σ̂2E was chosen as DF = N − L, but again no demonstration was given and the same brief citation of Deming59 was made. These two DF’s are obviously in conflict, because if one applies in the Anderson et al. case the formula of Britt and Luecke, one would have DF = 2N − L, where 2N is the total number of constraints. Since the correct estimate of the fundamental variance is a decisive issue in the present work, it is clarified in Appendix B (see Supporting Information) via a formal demonstration of the correct formulas for the number of degrees of freedom (DF) and the estimator of fundamental variance σ̂2E. This demonstration leads to eqs 46 and 47 below:

Then, it can be shown that eqs 32−34 produce the following expected values: E( Ẑ i) = Z̲ i ,

E( θ )̂ = θ̲

(38)

In other words, the expectation of MLC estimates in eqs 32−34 is the respective correct values; that is, MLC estimates in the implicit parameter estimation are asymptotically unbiased statistics with normal behavior. 4.3.1. Variance−Covariance Matrices. With eqs 36−38, variance−covariance matrices of estimates of parameters ( θ )̂ and variables ( Ẑ i) are obtained from eqs 39 and 40: COV ( θ )̂ = E(( θ ̂ − θ̲ )( θ ̂ − θ̲ )T ) COV ( Ẑ i) = E(( Ẑ i − Z̲ i)( Ẑ i − Z̲ i)T )

(39)

(i = 1, ..., nE) (40)

With eqs 32, 34, and 36−40, it can be shown that: COV( θ )̂ = σE2[Ξ]−1

(41)

COV ( Ẑ i) = σE2 Ψi

(42)

Ψi = W i−1{ W i − H iT M i H i + H iT M i J [Ξ]−1 J T M i H i} W i−1 i

(43)

i

If necessary, the covariance matrix between variable estimates of different experiments (E(( Ẑ i − Z̲ i)(Ẑ j − Z̲ j)T)) can also be produced (it is not zero). Equations 41 and 35b establish that the MLC estimates of parameters ( θ )̂ are also coherent, because COV ( θ )̂ goes to zero as nE increases. This can be proved from the fact that all matrices M i are positive definite, implying, when nE increases to infinity, that the limiting Ξ (eq 35b) is an infinite sum of positive definite matrices, becoming infinite with a vanishing inverse: nE

lim Ξ = lim

nE →∞

nE →∞

∑ J iT M i J i ≈ [∞] ⇒ i

lim COV ( θ )̂

nE →∞

= lim σE2[Ξ]−1 ≈ 0 nE →∞

nE

DF = (44)

∑ nFi − nθ i

This does not occur with COV ( Ẑ i) since eqs 42 and 43 attain a different limit when nE increases:

(46)

⎛ ⎞ nE 1 ⎟ ∑ (Z̲ iE − Ẑ i)T W (Z̲ iE − Ẑ i) σ̂E2 = SR2 ≡ ⎜⎜ nE ⎟ i ⎝ ∑i nFi − nθ ⎠ i

lim COV ( Ẑ i) = σE2 W i−1 − σE2 W i−1 H iT M i H i W i−1

nE →∞

(45)

(47)

The limiting COV ( Ẑ i), due to the definite positiveness of M i, is smaller in magnitude comparatively to the variance− covariance matrix for run i in eq 36. 4.3.2. Degrees of Freedom and Fundamental Variance Estimator in Implicit Estimation. The literature on parameter estimation and variance analysis is complete, and probably exhausted, when the focus is on multivariate regression with (linear or nonlinear) explicit models.57 In explicit models the separation of variables into input and output variables is feasible; that is, each output variable has a predictor dedicated to it, expressed only in terms of input variables and model parameters. A key assumption is that input variables do not bear experimental errors, but even if they do, errors are transferred to the output variables, which is the same thing. However, the scenario is somewhat nebulous when multivariate

In words: In the problem of implicit parameter estimation with several measured variables and several constraints per experiment, the best estimator for the fundamental variance is S2R in eq 47. Under assumptions I, II, and III, S2R is proportional to a Chi-Square with number of degrees of freedom DF in eq 46. This DF equals the total number of qualif ied constraints subtracted by the number of model parameters. Equation 46 agrees with the Britt and Luecke’s formula, despite the lack of explanation about their DF choice. On the other hand, the DF of Anderson et al.58 does not agree with eq 46 and is supposed to be wrong. 4.3.3. Variance Analysis for MLC Implicit Parameter Estimation. With eq 47 and calculating all matrices at the MLC estimates θ ,̂ Ẑ i (i = 1, ..., nE), it is possible to access estimators for variance−covariance matrices of estimated 9214

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

due to eclipsing by other more present reactions. Symptoms of these situations were ill-conditioned jacobians and lack of GNM convergence unless the problematic ChE constant is inactivated and locked. These 8 problematic or impossible estimations were as follows: (i) KH2S−MEA and KCO2−MEA at 120 °C (Table ) due to lack of MEA runs; (ii) simultaneously KH2S−DEA and KH2S−MDEA at 75 °C (i.e., estimation was only possible for one of them because DEA and MDEA actuated mixed in most runs); hence, KH2S−DEA was inactivated at 75 °C; (iii) KH2S−AMP was inactivated at 75 and 100 °C due to lack of data, and values 10.3785 and 4.1477 in Table 5 were guessed with healthy vicinal KH2S−AMP; (iv) KCO2−AMP was inactivated at 40, 75, and 100 °C due to similar issues as in the KH2S−AMP case. Therefore, it is advised that AMP ChE constants may exhibit bias due to lack of data to discriminate it from other amines, since AMP actuates in only a few runs and always mixed, in low concentration, with one or two other amines. In summary, the ChE constants for these 8 deficient cases were guessed and kept Inactive in the GNM searches for the remaining 24 active ChE constants. Table 4 reports the assignment of Active and

parameters and estimated variables, whose corrected values are given in eqs 41 and 42: ̂ ( θ )̂ = SR2[Ξ]−1 COV

(48)

̂ ( Ẑ i) = SR2 Ψ COV i

(49)

where Ψ i is given by eq 43 and is calculated at θ ,̂ Ẑ i (i = 1, ..., nE). Estimators for standard deviations of estimated values of the kth parameter (σ̂ θ k̂ ) and the jth variable in experiment i (σ̂ Ẑ ij) are given by the square root of the respective diagonal positions in COV̂ ( θ )̂ and COV̂ ( Z î ) as shown in eq 50: σ ̂ θk̂ =

SR2 [Ξ−1]kk ,

σ ̂ Zjî =

SR2 [Ψi]jj

(50)

Since the estimates θ̂k and Ẑ ji are approximately normal and uncorrelated with S2R, the t-Student theorem60 can be used to write confidence intervals for correct θk and Zji, with probability (1 − α)·100%, as shown in eqs 51a and 51b. Limits of ± 3 standard deviations for θk follow from eq 50, giving eq 52. θk̂ − t1 − α /2 SR2 [Ξ−1]kk ≤ θk ≤ θk̂ + t1 − α /2

SR2

−1

[Ξ ]kk

Table 4. Active (A) or Inactive (I) ChE Constants in Parameter Estimation

(51a)

Zjî − t1 − α /2 SR2 [Ψi]kk ≤ Zji ≤ Zjî + t1 − α /2 SR2 [Ψi]kk (51b)

θk̂ − 3 SR2 [Ξ−1]kk ≤ θk ≤ θk̂ + 3 SR2 [Ξ−1]kk

(52)

where t1−α/2 is the t-Student abscissa with probability (1 − α/ 2)·100% and ∑nE i nFi − nθ degrees of freedom (α = 0.1, 0.01, 0.001, etc). In Appendix B (see Supporting Information), it is shown that the Chi-Squares in eqs B18 and B24 are uncorrelated; hence, the F theorem60 establishes that their ratio, times the inverse ratio of the respective degrees of freedom, follows the F distribution with nθ and ∑nE i nFi − nθ degrees of freedom as written in eq 53: ( θ ̂ − θ̲ )T [Ξ]( θ ̂ − θ̲ ) nθ ·SR2

∑ nFi − nθ) i

(53)

Representing by ϕ1−α the F abscissa with probability (1 − α)·100% and (nθ, ∑nE i nFi − nθ) degrees of freedom (α = 0.1, 0.01, 0.001, etc), the Confidence Domain for correct parameters θ̲ with (1 − α)·100% of probability is as follows: ( θ̲ − θ )̂ T [Ξ]( θ̲ − θ )̂ ≤ nθ ·φ1 − α ·SR2

no.

40 °C

75 °C

100 °C

120 °C

KCO2−MEA KCO2−MDEA KCO2−DEA KCO2−AMP KH2S−MEA KH2S−MDEA KH2S−DEA KH2S−AMP

1 2 3 4 5 6 7 8

A A A I A A A A

A A A I A A I I

A A A I A A A I

I A A A I A A A

Inactive ChE constant showing that 3 ChE constants were successfully estimated at all 4 temperatures (KCO2−MDEA, KCO2−DEA, KH2S−MDEA), whereas other 3 were estimated at only 3 temperatures (KCO2−MEA, KH2S−MEA, KH2S−DEA). The worst cases, KH2S−AMP and KCO2−AMP, were statistically estimated at only 2 and 1 temperatures, respectively. In each one of these 4 isothermal estimations, all selected runs had experiment variables ( Ẑ i) estimated simultaneously with parameters ( θ )̂ guaranteeing all active constraints satisfied. Table 5 consolidates all results of estimations. It is seen that all estimated ChE constants decrease with T, as expected, since absorption reactions are exothermic. In order to concisely present illustrations generated by the estimation of ChE constants, analogous graphical objects for 40, 75, 100, and 120 °C are grouped. Figure 5 presents grouped histograms of logarithms of absolute residues of all constraints (active + inactive) after GNM convergence, confirming that the great majority of constraints (labels in Table 3) have attained negligible absolute residues below e−15, with only a few reaching e−10 (these few “tougher” residues are VLE and ChE residues). Figure 5a for 40 °C, with 8960 (active + inactive) constraints, reports its 4319 active constraints (Table 5), as satisfied by the 6172 estimated variables. Figure 5b for 75 °C, with 4788 (active + inactive) constraints, portraits its 2352 active constraints (Table 5) as satisfied by the 3364 estimated variables. Figure 5c for 100 °C, with 3192 (active + inactive) constraints, reports its 1456 active constraints (Table 5) as satisfied by the 2086 estimated variables. Figure 5d for 120 °C, with 2352 (active +

nE

→ F(ν1 = nθ , ν2 =

ChE constant

(54)

Formula 54 defines the interior of an 8D hyperellipsoid with (1 − α)·100% of probability for correct parameters. It can be reduced to 3D (ellipsoid) or 2D projections in the space of parameters by putting the nonprojected (correct) parameters at the estimated values (θk = θ̂k). 4.4. Results of Parameter Estimation for AGWA Model. The 8 ChE constants in Figure 2f and eq 1 were submitted to MLC estimation for 4 isothermal extracts at 40, 75, 100, and 120 °C, totaling 32 parameters to be estimated via 4 GNM searches. As shown in Table 1, due to inexistence of information in the 100 and 120 °C extracts, ChE constants KH2S‑AMP and KH2S−MEA cannot be estimated at 100 and 120 °C, respectively. Additionally, there are also another 6 ChE constants in certain isothermal extracts, whose identification was precarious, due to lack of data or, albeit with some data, 9215

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Table 5. Summary of Parameter Estimation for ChE in Eq 1 data temp

40 °C

320 experiments (nE) variables (37nE) 11840 active variables 6172 constraints (28nE) 8960 active constraints 4319 active ChE const 7 DF 4312 Estimated ChE Constants KCO2−MEA 367.61 9.74 × KCO2−MDEA KCO2−DEA 15.6 × KCO2−AMP 7.21 × KH2S−MEA 2.16 × 3.89 × KH2S−MDEA KH2S−DEA 9.80 × KH2S−AMP 195.64 Variance Estimator (Eq 47) S2R 11.146

75 °C 171 6327 3364 4788 2352 5 2347

105 106 103 103 106 108

0.845 4.80 × 66.6 × 30.31 11.17 19.4 × 7.04 × 10.38 25.343

103 103

103 106

100 °C

120 °C

114 4218 2086 3192 1456 6 1450

84 3108 1832 2352 1290 6 1284

0.098 9.31 55.43 10.75 1.74 27.37 15.5 × 103 4.15

0.031 0.604 1.545 5.72 0.027 5.41 0.84 2.26

10.011

3.191

inactive) constraints, portraits its 1290 active constraints (Table 5) as satisfied by the 1832 estimated variables. Figure 6 depicts grouped logarithmic plots of predicted versus observed loadings and partial pressures, variance− covariance matrices of estimated parameters, and also the estimated ChE constants at 40 and 75 °C with 99% probability confidence limits and ±3 standard deviation limits. Figure 7 depicts grouped predicted and observed partial pressures versus run number for CO2 and H2S. Figure 8 is analogous for predicted and observed loadings versus run number for CO2 and H2S. Figure 9 shows 99% confidence hyperellipsoids 3D projected for two triads of ChE constants per temperature level. Figure 10 conveys three singular images to be addressed later. 4.4.1. Estimation of ChE Constants at T = 40 °C. This isothermal set has the largest size with 320 runs, but the adherence of predicted onto observed values is good with a S2R of 11.146. Figure 6a expresses this via a logarithmic plot of predicted versus observed partial pressures and loadings of CO2 and H2S. The coalescence on the diagonal is convincing for values above 10−1, with some unavoidable scattering below it, thanks to the logarithmic scale. Figure 6e sketches absolute magnitudes in the 8 × 8 variance−covariance matrix of estimated ChE constants, numbered as in Table 4. Variances (and standard deviations) are small; the larger ones belong to KCO2−AMP, KH2S−DEA, and KH2S−AMP (parameters 4, 7, and 8). Figure 6f is a logarithmic plot of estimated ChE constants against parameter number, including 99% confidence limits (eq 51a) and ±3 standard deviation limits (eq 52). Limits are hard to see due to their relative tightness to the magnitude of ChE constants. Parts a and b of Figure 7 draw comparisons of predicted and observed partial pressures (bar) of CO2 and H2S against run number, where confidence limits for correct values were not plotted. Parts a and b of Figure 8 show similar plots for loadings (mol/mol) of CO2 and H2S. Concordance of predicted and observed values is evident for partial pressures in Figure 7a and b, whereas for loadings some discrepancies appear, mainly above 1.1 mol/mol, since the model can hardly produce loadings greater than 1 mol/mol. Parts a and b of Figure 9 depict 99% probability 3D projected confidence ellipsoids of correct values for two triads of ChE constants rendered by eq 54. The lengths of the principal axes of

Figure 5. Histograms of converged constraint residues after estimation of ChE constants at (a) 40 °C, (b) 75 °C, (c) 100 °C, and (d) 120 °C [MB: mass balance, SX: normalization, LV:VLE, CE:ChE, PP: partial pressures, LD: loadings].

ellipsoids traduce, with 99% probability, intervals of correct values of ChE constants. Thus, it is seen in Figure 9a that correct KCO2−MEA and KH2S−MEA lie in the intervals [367, 368] and [2163, 2165]. Figure 9b addresses a similar 99% 3D 9216

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Figure 6. Predicted vs observed partial pressures and loadings at (a) 40 °C, (b) 75 °C, (c) 100 °C, and (d) 120 °C; parameter 8 × 8 variance− covariance matrices at (e) 40 °C and (g) 75 °C; estimated ChE constants with 99% confidence limits and ±3 standard deviation limits at (f) 40 °C and (h) 75 °C.

4.4.2. Estimation of ChE Constants at T = 75 °C. This set comprises 171 runs. Estimation gives a S2R of 25.343, the largest in Table 5, perhaps influenced by the “forced” aggregation of runs from vicinal temperatures (70, 75, and 80 °C). Figure 6b shows the logarithmic plot of predicted versus observed partial pressures and loadings of CO2 and H2S with a convincing

confidence ellipsoid for KCO2−MDEA and KH2S−MDEA with values 1000−2000 times those of the MEA counterparts. Figure 10a shows, in addendum, the 99% 3D confidence ellipsoid of correct KH2S−AMP, KCO2−MEA, and KH2S−MEA, showing, with 99% probability, that KH2S−AMP lies in the interval [194, 198]. 9217

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Figure 7. Predicted and observed partial pressures vs run number for (a) CO2 at 40 °C, (b) H2S at 40 °C, (c) CO2 at 75 °C, (d) H2S at 75 °C, (e) CO2 at 100 °C, (f) H2S at 100 °C, (g) CO2 at 120 °C, and (h) H2S at 120 °C.

coalescence on the diagonal above 10−1. Visual scattering appears below 10−1, enhanced by the logarithmic scale. Figure 6g sketches magnitudes in the 8 × 8 variance−covariance matrix of estimated ChE constants. The larger variances belong to KCO2−AMP, KH2S−DEA, and KH2S−AMP (parameters 4, 7, and 8 in Table 4, all inactive in this search) with a maximum standard deviation of ≈1 for KH2S−DEA. Figure 6h is a logarithmic plot of

estimated ChE constants against parameter no., including 99% confidence limits (eq 51a) and ±3 standard deviation limits (eq 52), barely seen due to their relative tightness to the magnitudes of ChE constants. Parts c and d of Figure 7 represent predicted and observed partial pressures (bar) of CO2 and H2S against run number, where confidence limits for corrected variables (eq 51b) are also very tight and were not 9218

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Figure 8. Predicted and observed loadings vs run number for (a) CO2 at 40 °C, (b) H2S at 40 °C, (c) CO2 at 75 °C, (d) H2S at 75 °C, (e) CO2 at 100 °C, (f) H2S at 100 °C, (g) CO2 at 120 °C, and (h) H2S at 120 °C.

by eq 54, showing that correct KCO2−MEA, KH2S−MEA, and

plotted. Parts c and d of Figure 8 show analogous plots for loadings (mol/mol) of CO2 and H2S. The concordance of predicted and observed values is evident for partial pressures. Parts c and d of Figure 9 depict 99% probability 3D projected confidence ellipsoids for two triads of ChE constants rendered

KCO2−MDEA lie, respectively, in the intervals [0.6, 1.1], [10.9, 11.5], and [4802, 4803]. As before, MDEA ChE constants have values about 1000 times those of the MEA counterparts. 9219

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

Figure 9. 3D projected 99% probability confidence ellipsoids for triads of correct ChE constants: (a) 40 °C KH2SMEA, KCO2MDEA, KCO2MEA, (b) 40 °C KH2SMDEA, KH2SMEA, KCO2MDEA, (c) 75 °C KH2SMEA, KCO2MDEA, KCO2MEA, (d) 75 °C KH2SMDEA, KH2SMEA, KCO2MDEA, (e) 100 °C KH2SMEA, KCO2MDEA, KCO2MEA, (f) 100 °C KH2SMDEA, KH2SMEA, KCO2MEA, (g) 120 °C KH2SAMP, KH2SMDEA, KCO2MDEA, (h) 120 °C KH2SMDEA, KCO2DEA, KCO2MDEA.

4.4.3. Estimation of ChE Constants at T = 100 °C. This set comprises 114 runs. The adherence of predicted onto observed values is good with a S2R of 10.011. Figure 6c shows a logarithmic plot of predicted versus observed partial pressures and loadings of CO2 and H2S, with impressive coalescence on the diagonal above 10−1. Parts e and f of Figure 7 depict predicted and observed partial pressures (bar) of CO2 and H2S versus run number. Parts e and f of Figure 8 show analogous

plots for loadings of CO2 and H2S. To demonstrate the relative tightness of the 99% confidence limits for correct variables, the distributions of H2S loadings in Figure 8f were replotted in Figure 10b with superposition of the 99% probability confidence limits (eq 51b) and ±3 standard deviation limits (eq 52). The ∞ diverging lines in Figure 10b are associated with the confidence limits for runs without H2S. Parts e and f of Figure 9 depict the 99% probability 3D projected confidence 9220

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

to the same reason, adherence of predicted onto observed values is very good with a S2R of 3.191. Figure 6d shows a logarithmic plot of predicted versus observed partial pressures and loadings of CO2 and H2S, with a very good coalescence on the diagonal. Parts g and h of Figure 7 depict predicted and observed partial pressures (bar) of CO2 and H2S. Parts g and h of Figure 8 show similar plots for loadings of CO2 and H2S. The distribution of CO2 loadings in Figure 8g is replotted in Figure 10c with superposition of 99% probability confidence limits (eq 51b) and ±3 standard deviation limits (eq 52). The ∞ diverging lines are associated with confidence limits for runs without CO2. Parts g and h of Figure 9 depict 99% probability 3D projected confidence ellipsoids for two triads of ChE constants rendered by eq 54, showing that correct KH2S−AMP, KH2S−MDEA, and KCO2−MDEA lie, respectively, in the intervals [1.4, 3.2], [5.2, 5.7], and [0.5, 0.7]. 4.5. Correlation of Estimated ChE Constants with Absolute Temperature. The four values of estimated ChE constants for each absorption reaction (ln KEST (k = 1, ..., 8)) in k Table 5 were fitted by three-parameter formulas, eq 55, in terms of absolute temperature: B ln Kk = Ak + k + Ck ln(T ) (k = 1, ..., 8) (55) T Equation 55 allows differentiation of ChE constants, allowing access to thermal properties of formation of Complex species. Since it is desirable that all forms in eq 55 exhibit an everywhere convex behavior, monotonously decreasing with temperature, it is wise to take measures to avoid concavity inversions, maxima, minima, etc. Such anomalies may arise in some cases if a direct regression of eq 55 is attempted with only four data points. Equation 55 was then fitted, for each reaction of Complex formation, via a two step procedure: (i) First a twoparameter formula ln KBASE = Ak + (Bk/T) was adjusted with k the four estimated values of a ChE constant, giving always negative Ak and large positive Bk; (ii) Second, the residue RESk = ln KEST − ln KBASE was regressed against the logarithm of k k temperature with RESk = Ck ln(T), giving always a small Ck. Equation 55 was then built adding the terms of Ak, Bk, and Ck. In both (i) and (ii) stages, linear regressions weighted active ChE constants (Table 4) by a factor of 10 relative to inactive ones. The resulting coefficients are shown in Table 6. Since the Table 6. Fitted Coefficients of ln Kk(Liq) = Ak(Liq) + Bk(Liq)/T + Ck(Liq) ln T

Figure 10. 99% probability confidence domains for (a) 40 °C correct KH2SAMP, KH2SMEA, and KCO2MEA, (b) 100 °C correct predictions of H2S loadings vs run number, and (c) 120 °C correct predictions of CO2 loadings vs run number.

ellipsoids for two triads of ChE constants rendered by eq 54, showing that the correct KCO2−MEA, KH2S−MEA, and KCO2−MDEA lie, respectively, in the intervals [0, 0.3], [1.5, 2.0], and [9.2, 9.5]. 4.4.4. Estimation of ChE Constants at T = 120 °C. This set comprises 84 runs. A problem inherent to this case is its insufficiency of data to estimate MEA ChE constants. But, due

ChE constant

k

Ak(Liq)

Bk(Liq)

Ck(Liq)

KCO2−MEA KCO2−MDEA KCO2−DEA KCO2−AMP KH2S−MEA KH2S−MDEA KH2S−DEA KH2S−AMP

1 2 3 4 5 6 7 8

−45.3921 −57.5658 −63.1761 −23.4209 −39.0465 −53.9065 −72.0679 −16.7682

15975.8 22522.2 25230.0 9847.11 14592.7 21761.7 29316.0 6890.58

2.97 × 10−5 −4.50 × 10−5 −7.11 × 10−5 2.62 × 10−5 6.77 × 10−6 −2.36 × 10−5 −15.36 × 10−5 7.46 × 10−6

formulas adjusted via eq 55 refer to ChE constants with the choice of Pure Saturated Liquid (Liq) as standard state at T for Complex species, Table 6 writes these regressors as K(Liq) (T). k Figure C-1 in Appendix C (See Supporting Information) plots each regressed ChE constant in eq 55 against temperature with the respective estimated values in Table 5. Figures C-1a, C-1b, C-1c, and C-1d present such plots for the ChE constants of 9221

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

CO2 reactions with MEA, MDEA, DEA, and AMP. Figures C1e, C-1f, C-1g, and C-1h are analogous for the H2S ChE constants. The well-known exothermic absorption of CO2/H2S in aqueous alkanolamines is coherently reproduced by the proposed RVLE AGWA model, whose ultimate results for modeling AGWA separations are the formulas in Table 6. As shown in Figure C-1, the estimated (Table 5) and regressed (Table 6) ChE constants monotonously decrease with temperature, implying negative standard enthalpies of reaction (ΔH°), since (d ln K/dT) = (ΔH°/RT2) < 0. This behavior, exhibited by all chemical reactions in eq 1, is a prerequisite for the physical consistency of the proposed model.

ln Pksat(T ) = Ak* +

Bk* + Ck* ln T T

(61a)

ln Kk(Liq) = Ak(Liq) +

Bk(Liq) + Ck(Liq) ln T T

(61b)

ln Kk(g )(T ) = Ak(g ) +

Bk(g ) + Ck(g ) ln T T

(62)

Ck(g ) = Ck(Liq) + Ck*

(63)

Equation 61a is regressed over the EOS saturation curve of pure Complex k (ln PSat k vs T) via a set of subcritical T’s, where, for each T, ln ϕVk − ln ϕLk = 0 is solved for Psat k with liquid and vapor EOS fugacity coefficients of Complex k. Table 6 displays , B(Liq) , and C(Liq) in eq 61b, and Table D-1 in Appendix D A(Liq) k k k (see Supporting Information) shows the regressed parameters (g) (g) of eq 61a (A*k , B*k , C*k ) and coefficients A(g) k , Bk , and Ck for ChE constants with Complexes in state (g) using the SRK EOS. (Liq) Figure D-1, also in Appendix D, depicts ln Psat , and K(g) k , Kk k versus T for reaction of Complex k (k = 1, ..., 8) with coefficients in Tables 6 and D-1. Figures D-1a to D-1d display reactions of CO2 with MEA, MDEA, DEA, and AMP, while Figures D-1e to D-1h are analogous for H2S. Figure D-1 shows that with the (Liq) state eqs 1 are exothermic, but with the state (g) they become endothermic or weakly exothermic due to the large transition enthalpy between the (Liq) and (g) states. 5.1. Ideal Gas Properties of Formation of Complex Species. With the typical temperature dependence of the ideal gas heat capacities of species, the thermal properties of formation of Complex k in the (g) state such as H̅ °k (g)(T0), 3 −2 (g) C̅ (g) pk (T) = αk + βk·T + γk·T + ωkT , and G̅ k° (T0) can be accessed with eq 62. This was accomplished in Appendix E (see Supporting Information) with numerical results in Table E-1.

6. CONCLUDING REMARKS This work presents an equilibrium model founded in the ideas of the Chemical Theory for RVLE in AGWA (Acid Gas− Water−Amine) systems. This AGWA model was designed to reproduce absorption/stripping of CO2/H2S in/from aqueous phases with alkanolamines MEA, MDEA, DEA, and AMP. The formalism models AGWA absorption/stripping via multireactive VLE, with ChE reactions forming hypothetical nonvolatile and nonionic Complexes in liquid phase. The Pure Saturated Liquid standard state (Liq) was chosen for Complex species in the ChE constraints. There are 8 ChE reactions to model the VLE of CO2/H2S with water and the above 4 amines. The formalism was devised to accept wide ranges of temperature, pressure, liquid composition, and acid gas partial pressures. The set of VLE and ChE constraints is modeled with a traditional EOS (SRK & PR) applied to both

(57)

With reactants in state (g) and defining ΔG°k ,(Liq) and ΔG°k ,(g) as standard ΔG of reaction of eq 1k with Complex k, respectively, in standard states (Liq) and (g), one can write the following: °, (g) ΔGk°, (Liq)(T ) = G̅ k°, (Liq)(T ) − GAC.GAS (T ) − G̅ H°,2(g) ̅ O (T )

(58a)

°, (g) ΔGk°, (g)(T ) = G̅ k°, (g)(T ) − GAC.GAS (T ) − G̅ H°,2(g) ̅ O (T ) °, (g) (T ) − GAMINE ̅

(60)

Ak(g ) = Ak(Liq) + Ak*, Bk(g ) = Bk(Liq) + Bk*,

Equation 56 allows the conversion of the standard state for Complex k by eq 57:

°, (g) (T ) − GAMINE ̅

ln Kk(g ) = ln Kk(Liq) + ln Pksat(T )

where

(56)

G̅ k°, (g)(T0) = G̅ k°, (Liq)(T0) − RT0 ln Pksat(T0)

(59)

where K(Liq) and K(g) k k consider Complex k respectively in the standard states (Liq) and (g). With the T dependences of ln PSat k and ln K(Liq) in eq 61, eqs 62 and 63 are written for ln K(g) k k :

5. STANDARD STATE CONVERSION OF CHE CONSTANTS The ChE of Complex formations (eq 1) adopt the standard state of Pure Saturated Liquid (Liq) at T for Complexes (f °i (T) sat = Psat i (T)), P = Pi (T)) and Pure Ideal Gas (g) at T and 1 bar for Real species ( f i°(T) = P0 = 1 bar). The estimation of ChE constants in section 4.4 was based on isothermal VLE AGWA data only, without knowledge of thermal properties such as ideal gas heat capacity (C̅ (g) P (T)) and standard enthalpy of formation (H̅ °,(g)(T0)). However, the modeling of AGWA columns requires enthalpies calculated with ideal gas enthalpies added to residual enthalpies from the EOS. The ideal gas properties of Real species are obtained from literature tables, but for Complexes this is impossible. Formation properties H̅ °,(g)(T0) and C̅ (g) P (T) of complexes in the (g) state may be obtained from the temperature dependence of ChE constants rewritten for the (g) state as K(g) k (T), which can be obtained (T) in eq 55. The transformation K(Liq) (T) → from K(Liq) k k K(g) k (T) (k = 1, ..., 8) defines the standard state conversion of Complexes from Pure Saturated Liquid at T (Liq) to Pure Ideal Gas at T and 1 bar (g). The two standard states of a pure Complex k at T0 = 298.15 K can be linked via a transition from (g) state to (Liq) state in the following way: (i) Pure Ideal Gas k at (T0, P0 = 1 bar) with μ°k ,(g)(T0); (ii) Pure Ideal Gas k at (T0, ,(g) (T0) + RT0 ln PSat PSat k (T0)) with μk = μ° k k (T0); (iii) Since Sat Pk (T0) is small, this state is saturated vapor k at (T0, PSat k (T0)) Sat with μk = μ0,(g) k (T0) + RT0 ln Pk (T0), which is also the chemical potential of Pure Saturated Liquid k at (T0, PSat k (T0)); that is: μk°,(Liq)(T0 , PkSat(T0)) = μk°,(g )(T0) + RT0 ln PkSat(T0)

−ΔGk°, (g)(T ) −ΔGk°, (Liq)(T ) = + ln Pksat(T ) RT RT

(58b)

Combining eqs (58) and 57: 9222

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

ChE constants, an indication of a well-posed parameter estimation problem. In spite of the simplicity of this model and its small number of parameters (one per loading phenomenon), the adherence of predicted responsesi.e. partial pressures and loadings of acid gasesover a so extensive set of data is undoubtedly promising. Obviously with a more balanced or segregated set of data, some localized lacks of adherence of predictions can be improved. This will be the object of future work. The apparent limitations of loadings near 1 mol/mol exhibited by this model must be attributed to the characteristics of the data in question and not to the model itself. As any phenomenon-based predictor of loadings, this model can also predict high values of loadings if the concentration of alkanolamines is low enough and the partial pressure of acid gas is high enough. Another result accomplished in this work was the indirect estimation of ideal gas thermal properties of formation for the eight Complexes in this AGWA context. Such properties are essential for calculating enthalpies of reactive liquid phases with a conventional EOS, which, by their turn, are necessary to solve energy balances and heat effects in AGWA columns. In this case, a procedure based on standard state conversion of the 8 Complexes was implemented in Appendices D and E via Tables D-1 and E-1 (see Supporting Information). The intent was to rewrite the ChE constants using the Pure Ideal Gas standard state (g), such that ideal gas heat capacity, ideal gas enthalpy, and ideal gas Gibbs free energy of formation could be obtained for Complex species. Obviously, to model RVLE AGWA separations, the (Liq) choice of standard state of Complexes is still required in eq 8 with the formulas in Table 6. The (g) state is only relevant via Table E-1 to obtain ideal gas properties of formation for Complexes allowing calculation of liquid phase enthalpies and entropies with an EOS.

phases with classical mixing rules. An interesting point is that one does not need to establish bounds often encountered in commercial process simulators, such as, for example, a maximum %w/w of amines in the solvent. The reason is that the proposed model is ChE based with three reactants for each reaction, entailing that absorption of acid gas naturally does not evolve if water is not sufficiently present as a liquid. Hundreds of VLE experiments of AGWA systems were gathered from the Literature29−32,40−55 and used for estimating the 8 ChE constants of Complex forming reactions. Parameter estimation was applied in sequence to 4 isothermal extracts of VLE−AGWA experiments, namely 40, 75, 100, and 120 °C, respectively, with 320, 171, 114, and 84 VLE runs. This implied that a total of 32 ChE constants had to be estimated with the observed data; 8 ChE constants per isothermal extract. Each isothermal estimation of ChE constants used a MLC approach for parameter estimation, with nonlinear implicit models working with floating number of constraints and variables per experiment; that is, experiments are not required to have the same species and same AGWA reactions. Following Britt and Luecke,57 parameter estimation was developed based on a GNM implementation with full management of sparse objects. This algorithm was able to deal with large algebraic systems, sparse in blocks, within ordinary desktop computers. Despite its profile as a local optimizer, a remarkable point in this GNM implementation was that all four isothermal parameter estimations conducted here were successful despite the huge numbers of highly nonlinear constraints to satisfy and of variables to be estimated. Namely, 4319, 2352, 1456, and 1290 constraints and 6172, 3364, 2086, and 1832 variables, were simultaneously handled respectively for the isothermal extracts of 40, 75, 100, and 120 °C. Unfortunately, due to insufficient volume of nonredundant information capable to expose the individual role of some ChE reactions (mainly AMP reactions), 8 ChE constants from the original target of 32 ChE constants could not be statistically estimated. To negotiate with this obstacle, these 8 ChE constants were manually guessed and kept inactive during the estimation of the remaining active ChE constantsobviously, a crucial step was to identify and isolate the problematic ChE constants. As shown in Table 4, the inactive ChE constants were 5 (out of 8) AMP ChE constants, 1 (out of 8) DEA ChE constant, and 2 (out of 8) MEA ChE constants. The MDEA ChE constants were estimated with maximum comfort and confidence, a consequence of the massive presence of MDEA in the AGWA VLE Database. Statistical analysis of the estimation results was carried out for each isothermal extract in terms of fundamental variance estimator, variance−covariance matrices of estimates, confidence limits for correct individual parameters, and confidence ellipsoidal domains for correct triads of ChE constants. In view of the existence of certain confusion in the Literature concerning the correct determination of the Number of Degrees of Freedom (DF) and the Fundamental Variance Estimator (S2R) in the context of large nonlinear implicit parameter estimation problems, this work also presented in Appendix B (see Supporting Information) an analytical demonstration of the correct formulas for such entities. In this paper, the correct DF was relevant for conducting sound statistical analysis, which is based on confidence intervals, variance−covariance matrices, and ellipsoidal domains of (correct) ChE constants, all demanding a reliable estimator of the Fundamental Variance. It was observed that confidence limits are narrow relative to the absolute estimated values of



ASSOCIATED CONTENT

S Supporting Information *

GNM algorithm for implicit parameter estimation; correct number of degrees of freedom and the estimator of fundamental variance in multivariate parameter estimation with nonlinear implicit models with eqs B1−B27; regressions ln Kk(Liq) = Ak(Liq) + Bk(Liq) /T + Ck(Liq) ln T for complexes k (k = 1, ..., 8); formulas ln Kk(g) = Ak(g) + Bk(g)/T + Ck(g) ln T for complexes k (k = 1, ..., 8); formation properties of complexes in the ideal gas state (g); and selected illustrations at higher resolution. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*Tel.: +55-21-2562-7535. Fax: +55-21-2562-7637. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS J.L.d.M. and O.Q.F.A. gratefully acknowledge CNPq-Brazil for financial support. NOMENCLATURE

Abbreviations

AGWA = acid gas, water, and alkanolamine ChE = chemical equilibrium 9223

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

X̲ , Y̲ = vectors (n × 1) of mole fractions of Real species in liquid and vapor phases X̲ C = vector (nc × 1) of mole fractions of Complex species in liquid phase V = number of moles of vapor phase Z̲ , Z̲ E, Ẑ = vectors of correct, observed, and estimated experiment variables W i = diagonal weight matrix of experiment i

EOS = equation of state GNM = Gauss−Newton method LHS = left-hand side MED = method of estimated deviations MLC = maximum likelihood criterion MOD = method of observed deviations NG = natural gas NLP = nonlinear programming PR = Peng−Robinson RHS = right-hand side RVLE = reactive vapor−liquid equilibrium SRK = Soave−Redlich−Kwong TOV = true observed value UOV = untrue observed value VLE = vapor−liquid equilibrium

Greek Symbols

α = probability of a confidence interval be unable to include a correct value αj, βj, γj, δj, ωj = coefficients of molar heat capacity of species j as Pure Ideal Gas (g) Δαk, Δβk, ΔγkΔδk, Δωk = coefficients of heat capacity change of reaction of Complex k in state (g) α̲S = vector (ns × 1) of solute loadings (mol/mol amine) in liquid phase 0,(g) = heat capacity and enthalpy changes of ΔC(g) Pk , ΔHk reaction of Complex k in state (g) ΔG0,(Liq) = standard Gibbs free energy of reaction of k Complex k in standard state (Liq) ΔG0,(g) = standard Gibbs free energy of reaction of Complex k k in standard state (g) ϕ̂ k, ϕ1−α = fugacity coefficient of k and F abscissa with probability 1 − α γC̲ = vector (nC × 1) of activity coefficients of Complex species λ̲i = vector (nFi × 1) of Lagrange multipliers for experiment i Λ = lagrangian function μ̲, μ̲C = chemical potentials of Real (n × 1) and Complex species (nC × 1) (kJ/mol) Ξ = matrix in eq 26 Π̲ = Real species stoichiometric matrix (n × nc) for reactions in eq 1 Ψ i = matrix associated with experiment i in eq 37c σ2E = fundamental variance in the estimation process θ̲, θ ̂ = vector of correct and estimated model parameters (nθ × 1) χ2ν = chi-square with ν degrees of freedom

Math Symbols

a ,̂ a Ĉ = vectors of activities of Real species (n × 1) and Complex species (nC × 1) Ak, Bk, Ck = coefficients in ChE constant and vapor pressure in eqs 55, 61, 62 COV , COV (·) = variance−covariance matrix and variance− covariance operator C̅ (g) Pj (T) = molar heat capacity of species j in the standard state of Pure Ideal Gas (g) DF = number of statistical degrees of freedom in the estimation problem E(.), VAR(.) = expectancy and variance operators f i° = standard state fugacity of species i (bar) f ̂L, f ̂V = vectors (n × 1) of Real species fugacities in liquid and vapor phases (bar) F̲(Z̲ ,θ̲) = 0̲ = vector of experiment constraints with correct variables and parameters ΔG°(T) = vector (nC × 1) of standard Gibbs free energy of Complex reaction (kJ/mol) H i, J i = jacobians of constraints of experiment i in eqs 22a and 22b Kk, K̲ = ChE constant of reaction k; vector (nC × 1) of ChE constants L = liquid phase mole number LF = likelihood function M i = matrix in eq 26 associated with experiment i n, nC , nA, nS = no. of Real species, Complex species, amine species and solute species nθ , nZ , nF, nE = no. of parameters, variables, constraints, and AGWA experiments N̲ = vector (n × 1) of Real species mole numbers in experiment P, P̲SAT C (T) = pressure and vector of vapor pressures of Complex species (nC × 1) (bar) PDF = probability density function P̲S = vector (ns × 1) of partial pressures of solutes (bar) S̲A, S̲S = selection matrices of amines (nA x n) and solutes (nS x n) (m) S̲(m) θ , S̲Zi = search directions in parameter and variable spaces at iteration m of GNM S2R = weighted sum of squares of residues of variables T, TB = temperature and temperature of normal boiling point (K) t, tν = GNM damping factor and t-Student variable with ν degrees of freedom TC, PC, ω = critical temperature, pressure, and acentric factor



REFERENCES

(1) Sheilan, M. H.; Spooner, B. H.; Hoorn, E. Amine Treating and Sour Water Stripping, 3rd ed.; 2007. (2) Kundu, M.; Bandyopadhyay, S. S. Modelling VapourLiquid Equilibrium of CO2 in Aqueous N-Methyldiethanolamine through the Simulated Annealing Algorithm. Can. J. Chem. Eng. 2005, 83, 344. (3) Jamal, A.; Meisen, A.; Lim, C. J. Kinetics of Carbon Dioxide Absorption and Desorption in Aqueous Alkanolamine Solutions Using a Novel Hemispherical Contactor-II. Experimental Results and Parameter Estimation. Chem. Eng. Sci. 2006, 61, 6590. (4) Kvamsdal, H. M.; Rochelle, G. T. Effects of the Temperature Bulge in CO2 Absorption from Flue Gas by Aqueous Monoethanolamine. Ind. Eng. Chem. Res. 2008, 47, 867. (5) Lepaumier, H.; Picq, D.; Carrette, P. L. New Amines for CO2 Capture. I. Mechanisms of Amine Degradation in the Presence of CO2. Ind. Eng. Chem. Res. 2009, 48, 9061. (6) Sartori, G.; Ho, W. S.; Savage, D. W.; Chludzinski, G. R.; Wlechert, S. Sterically-Hindered Amines for Acid-Gas Absorption. Sep. Purif. Rev. 1987, 16 (2), 171. (7) Hook, R. J. An Investigation of Some Sterically Hindered Amines as Potential Carbon Dioxide Scrubbing Compounds. Ind. Eng. Chem. Res. 1997, 36 (5), 1779. 9224

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

(8) Samanta, A.; Bandyopadhyay, S. S. Absorption of Carbon Dioxide into Aqueous Solutions of Piperazine Activated 2-amino-2-methyl-1propanol. Chem. Eng. Sci. 2009, 64, 1185. (9) Mandal, B.; Bandyopadhyay, S. S. Simultaneous Absorption of CO2 and H2S into Aqueous Blends of N-Methyl-diethanolamine and Diethanolamine. Environ. Sci. Technol. 2006, 40, 6076. (10) Li, M. H.; Shen, K. P. Solubility of Hydrogen Sulfide in Aqueous Mixtures of Monoethanolamine with N-Methyl-diethanolamine. J. Chem. Eng. Data 1993, 38, 105. (11) Rebolledo-Libreros, M. E.; Trejo, A. Gas solubility of CO2 in aqueous solutions of methyldiethanolamine and diethanolamine with 2-amino-2-methyl-1-propanol. Fluid Phase Equilib. 2004, 218, 261. (12) Rebolledo-Libreros, M. E.; Trejo, A. Gas solubility of H2S in aqueous solutions of methyl-diethanolamine and diethanolamine with 2-amino-2-methyl-1-propanol at 313, 343, and 393K in the range 2.5− 1036 kPa. Fluid Phase Equilib. 2004, 224, 83. (13) Glasscock, D. A.; Critchfield, J. E.; Rochelle, G. T. CO2 Absorption/Desorption in Mixtures of Methyl-diethanolamine with Monoethanolamine or Diethanolamine. Chem. Eng. Sci. 1991, 11, 2828. (14) Alatiqi, I.; Sabri, M. F.; Bouhamra, W.; Alper, E. Steady-State Rate-Based Modelling for CO2/Amine Absorption-Desorption Systems. Gas Sep. Purif. 1994, 8, 3. (15) Ko, J. J.; Li, M. H. Kinetics of Absorption of CO2 in Solutions of MDEA+Water. Chem. Eng. Sci. 2000, 55, 4139. (16) Al-Baghli, N. A.; Pruess, S. A.; Yesavage, V. F.; Selim, M. S. A Rate-Based Model for the Design of Gas Absorbers for the Removal of CO2 and H2S Using Aqueous Solutions of MEA and DEA. Fluid Phase Equilib. 2001, 185, 31. (17) Bolhàr-Nordenkampf, M.; Friedl, A.; Koss, U.; Tork, T. Modeling Selective H2S Absorption and Desorption in an Aqueous MDEA-Solution Using a Rate-Based non-Equilibrium Approach. Chem. Eng. Proc. 2003, 43, 701. (18) Kabouche, A.; Meniai, A. H.; Bencheikh-Lehocine, M. Modeling the Absorption of Acidic Gases by Alkanolamine Solutions Incorporating Several Reactions. Chem. Eng. Technol. 2005, 28, 67. (19) Tobiesen, F. A.; Svendsen, H. F. Experimental Validation of a Rigorous Absorber Model for CO2 Post-combustion Capture. AIChE J. 2007, 53 (4), 846. (20) Aroonwilas, A.; Veawab, A. Characterization and Comparison of the CO2 Absorption Performance into Single and Blended Alkanolamines in a Packed Column. Ind. Eng. Chem. Res. 2004, 43, 2228. (21) Astarita, G.; Savage, D. W. Gas Absorption and Desorption with Reversible Instantaneous Chemical Reaction. Chem. Eng. Sci. 1980, 35, 1755. (22) Danckwerts, P. V. Absorption by Simultaneous Diffusion and Chemical Reaction. Trans. Faraday Soc. 1950, 46, 300. (23) Danckwerts, P. V. Gas−Liquid Reactions; McGraw-Hill: New York, 1979. (24) Billet, R.; Schultes, M. Predicting Mass Transfer in Packed Columns. Chem. Eng. Technol. 1993, 16, 1. (25) Bravo, J. L.; Rocha, J. A.; Fair, J. R. Comprehensive Model for the Performance of Columns Containing Structured Packings. Inst. Chem. Eng. Symp. Ser. 1992, 1, A489. (26) Brito, M. H.; von Stockar, U.; Bangerter, A. M.; Bomio, P. Effective Mass-Transfer Area in a Pilot Plant Column Equipped with Structured Packings and with Ceramic Rings. Ind. Eng. Chem. Res. 1994, 33, 647. (27) Leye, L.; Froment, G. F. Rigorous Simulation and Design of Columns for Gas Absorption and Chemical ReactionI. Packed Columns. Comput. Chem. Eng. 1986, 10, 493. (28) Rinker, E. B.; Ashour, S. S.; Sandall, O. C. Kinetics and Modelling of Carbon Dioxide Absorption into Aqueous Solutions of N-Methyl-diethanolamine. Chem. Eng. Sci. 1995, 50, 755. (29) Gabrielsen, J.; Michelsen, M. L.; Stenby, E. H.; Kontogeorgis, G. M. A Model for Estimating CO2 Solubility in Aqueous Alkanolamines. Ind. Eng. Chem. Res. 2005, 44, 3348.

(30) Gabrielsen, J.; Michelsen, M. L.; Stenby, E. H.; Kontogeorgis, G. M. Modeling of CO2 Absorber Using an AMP Solution. AIChE J. 2006, 52, 3443. (31) Gabrielsen, J.; Svendsen, H. F.; Michelsen, M. L.; Stenby, E. H.; Kontogeorgis, G. M. Experimental validation of a rate-based model for CO2 capture using an AMP solution. Chem. Eng. Sci. 2007, 62, 2397. (32) Chunxi, L.; Fürst, W. Representation of CO2 and H2S Solubility in Aqueous MDEA Solutions Using an Electrolyte Equation of State. Chem. Eng. Sci. 2000, 55, 2975. (33) Deshmukh, R.; Mather, A. A Mathematical Model for Equilibrium Solubility of Hydrogen Sulfide and Carbon Dioxide in Aqueous Alkanolamine Solutions. Chem. Eng. Sci. 1981, 36, 355. (34) Austgen, D. M.; Rochelle, G. T.; Peng, X.; Chen, C. C. Model of vapour-liquid equilibria for aqueous gas-alkanolamine systems using the electrolyte-NRTL equation. Ind. Eng. Chem. Res. 1989, 28, 1060. (35) Austgen, D. M.; Rochelle, G. T.; Chen, C. C. Model of vaporliquid equilibria for aqueous gas-alkanolamine systems. 2. Representation of hydrogen sulfide and carbon dioxide solubility in aqueous MDEA and carbon dioxide solubility in aqueous mixtures of MDEA with MEA or DEA. Ind. Eng. Chem. Res. 1991, 30, 543. (36) Posey, M. L.; Rochelle, G. T. A Thermodynamic Model of Methyldiethanolamine−CO2−H2S−Water. Ind. Eng. Chem. Res. 1997, 36, 3944. (37) Vallée, G.; Mougin, P.; Jullian, S.; Fürst, W. Representation of CO2 and H2S Absorption by Aqueous Solutions of Diethanolamine Using an Electrolyte Equation of State. Ind. Eng. Chem. Res. 1999, 38, 3473. (38) Prausnitz, J. M.; Lichtenthaler, R. N.; Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 2nd ed.; Prentice-Hall: Englewood Cliffs, 1986. (39) Barbosa, L. C.; Medeiros, J. L.; Araújo, O. Q. F. An Equilibrium Description for the Absorption of CO2/H2S with Aqueous Solutions of Ethanolamines. In ENPROMER 20054th Mercosur Congress on Process Systems Engineering, Rio de Janeiro, E-papers, 2005. (40) Shen, K. P.; Li, M. H. Solubility of CO2 in Aqueous Mixtures of MEA with MDEA. J. Chem. Eng. Data 1992, 37, 96. (41) Li, M. H.; Shen, K. P. Solubility of H2S in Aqueous Mixtures of MEA with MDEA. J. Chem. Eng. Data 1993, 38, 105. (42) Jou, F. Y.; Carroll, J. J.; Mather, A. E.; Otto, F. D. Solubility of Mixtures of Hydrogen Sulfide and Carbon Dioxide in Aqueous NMethyl-diethanolamine Solutions. J. Chem. Eng. Data 1993, 38, 75. (43) Dawodu, O. F.; Meisen, A. Solubility of CO2 in Aqueous Mixtures of Alkanolamines. J. Chem. Eng. Data 1994, 39, 548. (44) Jou, F. Y.; Mather, A.; Otto, F. Solubility of H2S and CO2 in Aqueous Methyl-diethanolamine Solutions. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 539. (45) Kent, R. L.; Eisenberg, B. Better Data for Amine Treating. Hydrocarbon Process. 1976, 55, 87. (46) Kuranov, G.; Rumpf, B.; Maurer, G.; Smirnova, N. VLE modelling for aqueous systems containing methyl-diethanolamine, carbon dioxide and hydrogen sulphide. Fluid Phase Equilib. 1997, 136, 147. (47) Kundu, M.; Bandyopadhyay, S. S. Solubility of CO2 in water + diethanolamine + N-methyldiethanolamine. Fluid Phase Equilib. 2006, 248, 158. (48) Mandal, B.; Bandyopadhyay, S. S. Absorption of Carbon Dioxide into Aqueous Blends of 2-amino-2-methyl-1-propanol and Monoethanolamine. Chem. Eng. Sci. 2006, 61, 5440. (49) Mokhatab, S.; Poe, W.; Speight, J. Handbook of Natural Gas Transmission and Processing, 1st ed.; 2006. (50) Murrieta-Guevara, F.; Rebolledo-Libreros, M. E.; RomeroMartínez, A.; Trejo, A. Solubility of CO2 in aqueous mixtures of diethanolamine with methyl-diethanolamine and with 2-amino-2methyl-1-propanol. Fluid Phase Equilib. 1998, 150−151, 721. (51) Sidi-Boumedine, R.; Horstmann, S.; Fischer, K.; Provost, E.; Fürst, W.; Gmehling, J. Experimental determination of carbon dioxide solubility data in aqueous alkanolamine solutions. Fluid Phase Equilib. 2004, 218, 85. 9225

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226

Industrial & Engineering Chemistry Research

Article

(52) Sidi-Boumedine, R.; Horstmann, S.; Fischer, K.; Provost, E.; Fürst, W.; Gmehling, J. Experimental determination of hydrogen sulfide solubility data in aqueous alkanolamine solutions. Fluid Phase Equilib. 2004, 218, 149. (53) Tomcej, R.; Otto, F. Absorption of CO2 and N2O into Aqueous Solutions of Methyl-diethanolamine. AIChE J. 1989, 35, 861. (54) Vrachnos, A.; Voutsas, E.; Magoulas, K.; Lygeros, A. Thermodynamics of Acid Gas-MDEA-Water Systems. Ind. Eng. Chem. Res. 2004, 43, 2798. (55) Benamor, A.; Aroua, M. K. Modeling of CO2 solubility and carbamate concentration in DEA, MDEA and their mixtures using the Deshmukh−Mather model. Fluid Phase Equilib. 2005, 231, 150. (56) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids, 4th ed.; McGraw-Hill: New York, 1987. (57) Britt, H. I.; Luecke, R. H. The Estimation of Parameters in NonLinear Implicit Models. Technometrics 1973, 15, 233. (58) Anderson, T. F.; Abrams, D. S.; Greens, E. A. Evaluation of Parameters for Nonlinear Thermodynamic Models. AIChE J. 1978, 24 (1), 20. (59) Deming, W. E. Statistical Adjustment of Data; Wiley: New York, 1943. (60) Himmelblau, D. M. Process Analysis by Statistical Methods; John Wiley & Sons: 1970.

9226

dx.doi.org/10.1021/ie302558b | Ind. Eng. Chem. Res. 2013, 52, 9203−9226