Phase Equilibria for Carbon Capture and Storage

1. Introduction

The worldwide anthropogenic carbon dioxide (CO₂) releases in atmosphere, accounting for over 80% of greenhouse gases (GHGs: > 81% CO₂, 10% methane, 7% nitrous oxide, ∼3% fluorinated gases [1, 2]), are emitted through numerous processes such as burning solid waste, trees, other biological materials, and especially fossil fuels (coal, oil, and natural gas) in power plants and energy production facilities, land use change, worldwide bushfires, and also as a result of certain chemical reactions in different industries (e.g., cement and steel factories) [3, 4, 5].

In addition, the Mauna Loa Observatory, Hawaii, USA forecast that the annual CO2 concentration will continue to rise in 2020, exceeding for the first time 415 ppm monthly [6], after the historical negative record of 415.39 ppm in May 2019. The higher and higher concentrations of CO2 together with methane, nitrous oxide, and fluorinated gases into atmosphere trap heat which causes the greenhouse effect producing the global warming and the climate change, triggering an escalation in the number and strength of natural disasters such tornados, floods, droughts, melting of ice from both poles, melting of mountain glaciers, and an increase in average sea levels [7, 8].

Among the many ways for carbon mitigation [9, 10], such as wind, solar, or hydro (waves) based electricity generation [4], waste management [11], chemical conversion of carbon [5], etc., one of the most mature technologies is the carbon capture and storage (CCS), which could contribute approximately 20% to CO2 emission reductions by 2050, as recommended by International Energy Agency (IEA) [12].

A key obstacle in the utilization on a large scale of CCS remain the operating and capital costs despite its great potential in many industries and petroleum refineries given their large CO₂ emissions [13, 14]. In the cases of industrial processes that generate pure or rich CO₂ gas streams the costs of CCS could be reduced but they are still an impediment to its deployment [13, 14].

In the last decade, many researchers focused almost exclusively on the development of improved sorbents with increased CO₂ capacity and/or reduced heat of regeneration [14]. However, recent results [3, 5] show that whilst equilibrium CO₂ capacity is a key determinant of process performance, phase equilibria, transport properties (e.g., viscosity, diffusion coefficients, etc.) and other thermophysical properties (e.g., heat capacity, density, etc.) have a significant effect on the capital cost, and thus on the price of the carbon captured.

In this context we recently focused on the phase behavior research of physical solvents for carbon dioxide capture. Phase equilibria (PE) at high-pressures (HP) of carbon dioxide with different classes of organic substances, such as alkanes [15], cycloalkanes [16, 17, 18], ethers [19, 20], alcohols [21, 22, 23, 24], esters, ketones, were investigated to illustrate the functional group effect on the solvent ability to dissolve CO₂. As the experiments are usually expensive and very time consuming, equations of state (EoS) models are the most common approach for the correlation and/or prediction of phase equilibria and properties of the mixtures [25]. In this chapter we will also show modeling results with the General Equation of State (GEOS) [26, 27], Peng–Robinson (PR) [28] and Soave–Redlich–Kwong (SRK) [29] equations of state (EoS) coupled with both classical van der Waals (two-parameter conventional mixing rule, 2PCMR). Our approach is rather predictive, and we are interested in the global phase behavior of the systems, not only on pressure-composition diagrams, for instance. Therefore, phase diagrams construction is also explained.

2. Experimental method and setup for measuring PE at HP

One of our laboratory experimental setup for measuring phase equilibria at high pressure consists of three modules: equilibrium cell, sampling and analyzing systems and control systems (Figure 1). The main component of the experimental apparatus is the high-pressure windowed cell with variable volume [30, 31, 32], coupled with the sampling and analyzing system [33, 34]. The sampling system consists of two high-pressure electromechanical sampling valves, namely the rapid on-line sampler injector (ROLSI™, MINES ParisTech/CEP-TEP – Centre énergétique et procédés, Fontainbleau, France [35]). The ROLSI valves are connected to the equilibrium cell and to a gas chromatograph (GC) through capillaries. The expansion chamber of the sampler injector is heated with a heating resistance, so the liquid samples are rapidly vaporized. A linear resistor coupled to an Armines/CEP/TEP regulator is used to heat the transferring lines between ROLSI and the GC. The GC (Perichrom) is equipped with a thermal conductivity detector, TCD, and a HP-Plot/Q column 30 m long and 0.530 mm diameter. Helium is the GC carrier gas at a flow rate of 30 mL/min. The setup is completed with a syringe pump Teledyne ISCO model 500D.

The working procedure to determine isothermal phase equilibria is briefly described. Firstly, the entire internal loop of the apparatus including the equilibrium cell is flushed several times with carbon dioxide. Any trace of carbon dioxide is then evacuated from the cell and lines with a vacuum pump. The next step is to load the organic substance in the equilibrium cell, which is previously degassed by using a vacuum pump and vigorously stirring. The lighter component (in this case CO2) is also filled with the syringe pump into equilibrium cell and the pressure is set to the desired value. The experiment continues by heating the cell to the required temperature. The mixture in the cell is stirred for a few hours to facilitate the approach to an equilibrium state. Then the stirrer is switched off for about one hour allowing complete separation of the coexisting phases. Samples from both liquid and vapor phases are withdrawn by ROLSI valves and analyzed with the GC. To check the repeatability, at least six samples of the liquid phase are normally analyzed at the equilibrium temperature and pressure. The sample sizes being very small, the equilibrium pressure in the cell remains constant.

The calibration of the TCD for CO2 and organic substance is done by injecting known amounts of each component using gas chromatographic syringes. Calibration data are fitted to quadratic polynomials to obtain the mole number of the component versus chromatographic area. Then correlation coefficients of the GC calibration curves are calculated.

The uncertainties all variables and properties are estimated [21]. For the phase equilibrium compositions, the relative uncertainty of the mole fraction in the liquid and vapor phases are calculated using the procedure given by Scheidgen [37] and Chirico et al. [38].

A different working procedure is used when measuring critical points of systems at high-pressures. Thus, we start with a homogenous phase with unknown composition at random temperature and pressure. The pressure is modified by varying the volume of the cell with the manual pump to determine if we obtain a bubble or a dew point. If we obtain a bubble point, the temperature is slowly increased until the first dew point is observed, then the pressure is increased to a homogeneous phase and the composition is determined by sampling. Then the pressure is very slowly decreased until the first drops of liquid are observed. At this point, the temperature is slowly decreased simultaneously with reducing the volume, so the system is at the limit between homogeneous (single phase)-heterogeneous (two phases). The decreasing of temperature continues until the first gas bubbles are observed. The procedure is then repeated by introducing new amounts of CO2 and slowly cooling.

3. Modeling

Equations of state (EoS) models are the most common approach for the correlation and prediction of phase equilibria and properties of the mixtures. Cubic and generalized van der Waals equations of state were actively studied, since van der Waals proposed his famous equation in 1873 [39]. Although the cubic EoSs have their known limitations [39, 40, 41], they are frequently used for practical applications. They offer the best balance between accuracy, simplicity, reliability, and speed of computation [39, 40, 41]. Therefore, cubic equations of state remain important and easy tools to calculate the phase behavior in many systems [40], even for complex mixtures like petroleum fluids [25].

The general cubic equation of state (GEOS) has the form:

The four parameters a, b, c, d for a pure component are expressed by:

The temperature function used is:

with the reduced temperature Tr=T/Tc. The expressions of the parameters Ωa, Ωb, Ωc, and Ωd are:

with

where αcis the Riedel’s criterion.

It can be easily noticed that a, b, c, d coefficients of the cubic GEOS equation are in fact functions of the critical data (T_c, P_c, and V_c), m, and α_c parameters.

The cubic GEOS equation is a general form for all the cubic equations of state with two, three, and four parameters [42]. This is the meaning of the statement cubic “general equation of state” used for GEOS. To obtain the parameters of the Peng–Robinson (PR) or Soave–Redlich–Kwong (SRK) equations of state from the Eqs. (2–5) we set the following restrictions [42]: Ωc=−2Ωb2and Ωd=−Ωbor Ωc=−Ωb22and Ωd=−Ωb2respectively. It results:

or

giving BPR=0.2296and ZcPR=0.3074or BSRK=0.2467.

The equations for B can be solved iteratively, starting with an initial approximation of B in the right hand term. The corresponding value for Ωa, Ωb, Ωc, Ωd are given by Eq. (4). The Ωa, Ωb, Ωc, Ωd parameters are compound dependent for GEOS, while for SRK and PR are universal. It must be also pointed out that in GEOS the value of Zc is the experimental one.

The coefficients a, b, c, d were obtained for mixtures using the classical van der Waals two parameters conventional mixing rules (2PCMR for a, b) extended correspondingly for c and d:

Generally, negative values are common for the c parameter of pure components.

The GEOS parameters m and α_c of each component were estimated by constraining the EoS to reproduce the experimental vapor pressure and liquid volume on the saturation curve between the triple point and the critical point.

The phase equilibrium is represented using the Gibbs energy minimization criterion, which can be equivalently written as iso-fugacity criterion [43]:

where f is the fugacity of the component i in the liquid (L) or vapor (V) phase.

Using the fugacity model for both phases, this equation becomes:

where ϕiLand ϕiVare the fugacity coefficients of component i in the liquid and vapor phase respectively. Using the equilibrium constant, Ki=yixi, Eq. (10) becomes:

The fugacity coefficient of component i in the mixture is calculated from any equation of state.

The vapor–liquid equilibrium calculation for bubble p or dew p type problems involves solving the following system of nonlinear equations:

or

respectively.

The calculations were made using the software package PHEQ (Phase Equilibria Database and Applications), developed in our laboratory [44], and the GPEC (Global Phase Equilibrium Calculations) [45]. In our in-house software (CRIMIX), the critical curves are calculated using the method proposed by Heidemann and Khalil [46] with the numerical derivatives given by Stockfleth and Dohrn [47].

4. Phase diagrams

The phase diagram shows the domains occupied by the different phases of a system, the boundaries that separate these regions and the special points of the system, as a function of two independent variables. A practical choice of these variables is that of pressure (P) and temperature (T), which can be measured experimentally. The phase law states that for a one-component system, the coexistence curves, whether liquid-vapor, solid-vapor or solid–liquid, are monovariants. The curves separate in the phase diagram the domains of existence of each phase (vapor, liquid, solid).

The most known classification of phase diagram types was proposed by Scott and van Konyneburg [48]. They applied the van der Waals equation combined with van der Waals-type mixing rules in binary systems and quantitatively predicted almost all types of phase equilibria of fluids known from experiments. The results were presented in a global phase diagram and its pressure–temperature (P–T) projections. Scott and van Konyneburg [48] classified the different types of phase diagrams taking into account the nature of P–T projections and, in particular, the presence or absence of three-phase lines and azeotropic lines, as well as how critical lines connect with these. According to Scott and van Konyneburg’s classification, six main types of phase diagrams can be distinguished. The first five types of fluid behavior were calculated by [48] with the above equation, the sixth type being calculated with other equations of state [32]. More recently, Privat and Jaubert [49], presented an updated version of classification of the phase diagrams in binary systems together with the transitions between the various types of systems. The six main types of phase diagrams are as follows (Figure 2 [32]):

• Type I – presents a critical liquid-vapor curve which connects the critical points of the two pure components;

• Type II – presents two critical curves. It has a critical curve as in the case of type I and, in addition, a liquid–liquid critical curve that starts at an upper critical endpoint (UCEP) and evolves rapidly at high pressures and a three-phase liquid–liquid-vapor (LLV) equilibrium curve ending in UCEP;

• Type III – presents three critical curves. A critical liquid-vapor curve that starts at the critical point of one of the pure components and ends in a UCEP where the curve of the three liquid–liquid-vapor phases ends; another liquid-vapor critical curve starting from the critical point of the other pure component and connecting with a liquid–liquid critical curve at very high pressures;

• Type IV – presents two distinct liquid-vapor critical curves starting from each of the critical points of the pure components and a liquid–liquid–vapor equilibrium curve of the three phases bordering the two liquid–vapor critical curves, the two critical points where it intersects representing a UCEP and, respectively, a lower critical endpoint (LCEP). Also, in this type of diagram, a liquid–liquid critical curve is found at low temperatures, ending in another UCEP, in which an equilibrium curve of the three LLV phases ends;

• Type V – presents a critical curve similar to that of type IV, but without the liquid–liquid critical curve and without the equilibrium curve of the three-phases at low temperatures;

• Type VI – presents a critical curve as in the case of type I, but also a critical liquid–liquid curve that stretches between an LCEP and a UCEP, which also border the equilibrium curve of the three phases.

Subsequently, the classification of phase diagrams was enriched with several new types of critical curves [32]:

• Type VII - this phase diagram shows the same type of liquid–liquid critical curve as type VI, but the liquid-vapor critical curve is composed of two sections ending in a UCEP, respectively in an LCEP, which in turn border a liquid–liquid-vapor curve of the three phases;

• Type VIII -has a critical curve as in the case of type III which additionally presents a liquid–liquid critical line that starts at a lower critical limit point, LCEP and then goes to high pressures [30];

• Type Vm - the critical curve is interrupted by a curve of the three phases and has a maximum and a minimum in pressure.

In the case of binary systems, especially, the sections P-X at T = constant or T-X at P = constant are widely used. Such sections are constructed as in Figures 3 and 4. In P–T projections, the axis of the composition is hidden by the pressure–temperature plane. The vapor pressure curve for component 1 is denoted 1 and for component 2, 2. In Figure 3, a P–T projection of the type I diagram is presented, but, in addition, for a composition, X, we represented the bubble pressure curves (BPC) and dew pressure curves (DPC) that intersects at a point on the critical liquid-vapor curve, where the liquid and vapor phases become identical.

Each critical point on the critical liquid-vapor curve corresponds to a certain composition. The points near the critical point of component 1 are points in the mixture where the concentration of component 1 is very high, while the points located near the component 2 are critical points of the mixture in which its concentration is very high.

To construct a P-X diagram at T = constant (Figure 3), for example T₁, a line parallel to the pressure axis is drawn. This line, corresponding to the temperature T₁, firstly meets the vapor pressure curve of component 2, then the dew pressure curves (DPC), respectively the boiling pressure curves (BPC) and, finally, the vapor pressure curve of component 1.

All mixtures of components 1 and 2 are represented at low pressure by a single vapor phase (single-phase domain) and at high pressure by a single liquid phase (single-phase domain). The biphasic domain, delimited by the liquid curve, respectively by that of the vapors, corresponds to the equilibrium state between phases (liquid + vapors).

To construct a T-X diagram at P = constant (Figure 4), one can proceed similarly, but draw a line parallel to the temperature axis.

The complexity of P–T and P-X projections increases as the phase behavior is more sophisticated. As examples, the P–T and P-X projections for a type III phase diagram are shown in Figures 5 and 6 for two constant temperatures: one lower than the critical point of the pure component 1 and the other located between the critical point 1 and the UCEP.

5. Predicting the phase behavior of systems containing carbon dioxide of interest for CCS

Binary systems containing carbon dioxide and organic substances are of interest for CCS both from fundamental point of view and for the ability of different substances to dissolve CO₂. Different classes of organic substances are studied to explore the functional group effect on the solvent ability to dissolve CO₂. Several examples are discussed in the next sections.

5.1 Phase diagram predictions for the carbon dioxide +1-dodecanol binary system

As explained in previous section, the phase diagrams can be calculated using equations of state. The experimental data can be correlated at each temperature and sets (k_ij, l_ij) of optimized binary interaction parameters (BIPs) are obtained. Each set then can be used to calculate the phase diagram of that system generating as many diagrams as the temperatures used.

One example is shown in Figure 7 where for the carbon dioxide +1-dodecanol binary system [22] are plotted the critical curves calculated with the PR/2PCMR model using the optimized binary interaction parameters obtained at 293.15 K (k₁₂ = 0.1277; l₁₂ = −0.0026), 303.15 K (k₁₂ = 0.1266; l₁₂ = −0.0068), 313.15 K (k₁₂ = 0.1234; l₁₂ = −0.0055), 333.15 K (k₁₂ = 0.1535; l₁₂ = −0.0193), and 353.15 K (k₁₂ = 0.1380; l₁₂ = −0.0100).

Another possibility is to use one constant set of BIP to calculate the phase diagram. In the same figure, the critical curve (dark red continuous line) of the system is also calculated with a constant set of BIPS obtained by the averaging the optimized BIPs for the five temperatures mentioned (k₁₂ = 0.1278; l₁₂ = −0.0006).

The calculation with the average value of BIPs is not representing well the experimental critical data (blue full circles), therefore a modified set of BIPs (k₁₂ = 0.0900; l₁₂ = −0.0100) was obtained by a trial and error procedure [22]. This set of parameters is predicting reasonably well the type phase behavior of the carbon dioxide +1-dodecanol binary system.

5.2 Phase diagram predictions for carbon dioxide + different classes of substances at high-pressures

A more severe test for an equation of state is to use the BIPs obtained for a mixture to calculate the phase diagram for other ones. Figure 8 illustrates the predictions of phase diagrams for five binary systems containing carbon dioxide and different classes of organic substances containing four C atoms: one normal alkane (n-butane, n-C₄H₁₀), one primary (1-butanol, C₄H₁₀O) and one secondary alcohol (2-butanol, C₄H₁₀O), one ester (ethyl acetate, C₄H₈O₂), and one ether (1,2-dimethoxyethane (1,2-DME), C₄H₁₀O₂).

The calculations were made with the SRK/2PCMR EoS. In all cases, the model predicts type II phase behavior, even if only for carbon dioxide +1-butanol and carbon dioxide +2-butanol there is experimental evidence that they exhibit liquid–liquid immiscibility [50]. One unique set of binary interaction parameters tailored for the carbon dioxide +2-butanol system was used to predict the phase behavior of the other four binary systems. Thus, SRK EoS was firstly used in a semi-predictive approach to obtain a set of binary parameters yielding good results for the binary system carbon dioxide +2-butanol (including VLE in the entire temperature range, critical curve). The set of binary parameters was calculated using the k₁₂-l₁₂ method [50] to obtain the experimental value of the vapor-liquid critical pressure maximum (CPM) simultaneously with the temperature of UCEP. The experimental temperature of UCEP and CPM have been traced by paths in the k₁₂-l₁₂ diagram and their intersection has given the values of the parameters [50]. The values of the interaction binary parameters fulfilling these requirements are k₁₂ = 0.0200 and l₁₂ = −0.1110. This set of BIPs was used then to predict the phase behavior for the carbon dioxide + n-butane (black curve), + ethyl acetate (green curve), + 1,2-DME (pink curve), and + 1-butanol (olive curve). Except for the carbon dioxide + n-butane binary system, the critical curves are remarkably well predicted for all mixtures. It should be noted that for the carbon dioxide + ethyl acetate system the experimental critical data are scattered (orange and brown full circles), while for the carbon dioxide +1-butanol system the model predicts the CPM at a higher temperature than the experimental one. The predicted UCEPs are increasing in temperature in the order: T_CO2 + _n-butane < T_{CO2 + ethyl acetate} < T_CO2 + _1,2-DME < T_CO2 + _2-butanol < T_CO2 + _1-butanol.

The critical lines from Figure 8 indicate that the binary of carbon dioxide +1-butanol displays the larger range of immiscibility, followed by carbon dioxide +2-butanol, carbon dioxide +1,2-DME, and carbon dioxide + ethyl acetate, compared with the binary of carbon dioxide + n-butane, the corresponding n-alkane with the same carbon number.

5.3 Comparisons of phase behavior with different models

The carbon dioxide +2-butanol binary system can be used as a model system as experimental data are available for the entire range of critical temperatures and pressures, including the UCEP [50]. Figure 9 shows the predicted critical lines for the carbon dioxide +2-butanol system with three models: GEOS, SRK, and PR using a single set of BIPs for each model [50]. The BIPs were obtained using the k₁₂-l₁₂ method explained in previous sub-section. It can be observed that all models predict very well the liquid-vapor critical curve. In addition, GEOS equation reproduces very well the CPM and the corresponding temperature, while SRK and PR predict a similar value for the CPM but shifted to a higher temperature. The models show a comparable behavior for the liquid–liquid critical curves, the only difference being the predicted UCEPs. However, the models show bigger difference in the critical pressure-compositions and critical temperature compositions projections of the phase diagrams, as can be seen in Figure 10. The prediction results by GEOS are better than PR and SRK respectively. The differences in the predictions will than be reflected in the prediction of the bubble-point and dew-point curves [50].

6. Conclusions

The phase behavior of selected hydrocarbons and functional groups substances (alcohols, ethers, esters) is investigated to assess their capability as solvents/cosolvents for CO2 capture. Simple models as cubic equations of state have the capability to accurately predict the phase behavior of carbon dioxide + different organic substances binary systems. Phase diagrams are important tools in assessing the global behavior. The type of phase behavior is providing essential information for the operating domain of binary system.

Acknowledgments

This work was supported by a grant of Ministry of Research and Innovation, CNCS - UEFISCDI, project number PN-III-P4-ID-PCE-2016-0629, within PNCDI III.

Conflict of interest

The authors declare no conflict of interest.