Scalable simulation, design, and optimization of the CO₂ capture processes start with modeling of the thermodynamic properties, specifically vapor-liquid equilibrium (VLE) and chemical reaction equilibrium, as well as calorimetric properties. Accurate modeling of thermodynamic properties requires availability of reliable experimental data and a good model and accurate model to simulate VLE data. A thermodynamic property model capable of accurate representation of the vapor-liquid equilibrium (VLE) of the aqueous MEA-CO₂ system is essential for a successful computer simulation of the process. In the last few decades, considerable progress has been made in modeling VLE of the acid gas (CO₂ and H₂S) in aqueous alkanolamine systems, including the aqueous MEA-CO₂ system.

Many models have been applied to predict the adsorption of CO₂ by MEA, some of them have been successful and others can be applied in a specific range of parameters such as temperature and concentration. The extended UNIQUAC model sticks out to be one of the best and they have been fewer applications in this respect despite its simplicity and the small number of its parameters. In this paper, we will be using this model using a different approach of regressing its parameters, and not based on the traditional regression methods but on multiple objective function using differential evolution algorithm.

Thermodynamic framework

Vapor liquid equilibrium: The processes discussed involve both chemical equilibria and multi-component phase equilibria. The liquid phase comprises both molecular species and ionic species, which makes the modelling non-trivial. The chemical reactions taking place in the liquid phase for MDEA-CO₂-H₂O can be expressed as:

Water ionization

2H₂O↔ H₃O⁺+ OH⁻ (1)

Dissociation of carbon dioxide

CO₂ + 2H₂O↔ H₃O⁺ + HCO₃⁻ (2)

Dissociation of bicarbonate

HCO₃⁻ + H₂O↔ H₃O⁺ + CO₃²⁻ (3)

Dissociation of protonated amine

MEAH⁺ + H₂O↔ H₃O⁺+ + MEA (4)

MEACOO⁻ (carbamate) reversion to bicarbamate

C₃H₆NO₃⁻ (MEACOO⁻ ) + H₂O↔C₂H₇NO(MEA) + HCO₃⁻ (5)

The model

Extended UNIQUAC: The extended UNIQUAC model expresses the excess Gibbs energy as the combination of three terms contributing to the total excess Gibbs energy: the entropic and enthalpy terms of the original UNIQUAC equation to consider the non-electrostatic interactions and, the electrostatic term (Debye Hückel):

(6)

The first term i.e. the entropic term, is to describe the deviation from ideality due to the shapes and sizes of individual species in the solution and it is determined by the composition. This term is given as

(7)

Where z is the coordination number arbitrarily set to 10. Volume fraction Ø_i and area fraction θ_i are

x_i is the mole fraction of component i. The volume parameter r_i and the surface area parameter q_i are treated as adjustable parameters in this work.

The residual or the enthalpy term, is the other short range term of the UNIQUAC equation which is meant to take into account the energetic interactions between like and unlike species,

(8)

Where ψ_ji is

The adjustable interaction parameters (u_ji) of the UNIQUAC enthalpy term are assumed to be temperature dependent and are fitted to the following function of temperature

(9)

The Debye-Hückel expression used is the simplification of the original term given by Debye and Hückel for the electrostatic contribution to the excess Gibbs energy:

(10)

x_wand M_w are the mole fraction and molar mass of water respectively. The parameter b is considered to be a constant equal to 1.50 (kg/mol)^1/2. A is the Debye-Hückel constant given by

Where F (Cmol^-1) is the Faraday’s constant, NA (mol^-1) is Avogadro’s number, ε0 (C2 J^-1m^-1) is the vacuum permittivity, R (Jmol^-1 K^-1) the gas constant and T (K) is the temperature. d (kgm^-3) and εr are the density and the relative permittivity of the solution respectively.

The Debye-Hückel parameter is considered temperature dependent which for the temperature range of 273.15 to 383.15 K can be written as

I (mol/kg H₂O) is the ionic strength expressed as:

(11)

Where m_i (mol/kgH₂O) and z_i are the molality and the charge number of ionic species i respectively.

There are no adjustable parameters in the Debye-Hückel term and the Debye-Hückel parameter is based on the density and dielectric constant of pure water. This means that the effect of MEA and other solute species on the dielectric constant of the solution is not considered and also by considering water as the single solvent, the density of pure water rather than that of the solution can be used. This strategy simplifies the calculations while the reasonable precision of the model is sustained. The vapor phase fugacities are calculated using the Soave-Redlich-Kwong (SRK) equation of state, no parameters to be adjusted for the SRK-EOS (Equation of State). Now that all the terms contributing to the total excess Gibbs energy are defined, the activity coefficients of the molecular and ionic species are obtained by partial molar differentiation of the Gibbs energy expressions

(12)

The symmetrical activity coefficient for water and the asymmetrical activity coefficients for the ions can be calculated as:

(13)

(14)

The Gibbs-Helmholtz equation is used for calculation of the symmetrical excess enthalpy of MEA-water solutions. The equation defines the temperature dependence of the excess Gibbs energy and therefore also of the activity coefficients:

(15)

The experimental values and those calculated by the model are indicated by “exp” and “calc” respectively. P (bar) is the bubble point pressure, HE (Jmol^-1) the excess enthalpy and ΔG° (J mol^-1) is the change in the chemical potential by transferring one mole from liquid to solid phase and R (mol/kg H₂O) is the gas constant.

Differential evolution

Storn and Price proposed Differential Evolution (DE) as a family of algorithms to solve real-parameter optimization problems [1]. The variants of DE are differentiated from each other by varying the mutation and/or the recombination operation within a common framework. However, in this study we used the DE/rand/1/exp variant of DE. Therefore, we describe DE with this variant.

DE works with a population of individuals xⁱ _G, i = 1, 2, . . . , p each representing a solution to the problem. DE individuals are encoded as real vectors of size which is the dimension of the problem. The number of individuals in a population is called population size and is denoted by P and the generation number is denoted by G. The initial population, P₁ is created by randomly creating the vectors in appropriate search ranges. Then the fitness score of each individual is calculated through evaluation.

DE practices random parent selection regardless of their fitness values. In every generation, each individual xi G gets a chance to become the principal parent and to breed its own offspring mating with other randomly chosen auxiliary parents. Formally, for every principal parent xⁱ _G, i= 1, 2, . . . , p, three other auxiliary parents x_G^r1 , x_G^r2, x_G^r3 are selected randomly such that r1, r2, r3 ∈ {1, 2, . . . , p} and i ≠ r1 ≠ r2 ≠ r3. Then these three auxiliary parents participate in differential mutation operation to create a mutated individual x_G^must as follows:

Where F is the amplification factor, a real-valued control parameter chosen from [0.1, 1.0] [2]. Subsequently, the mutated vector, x_G^must, participates in exponential crossover operation with the principal parent x_Gⁱto generate the trial individual or offspring x_G^child. Exponential crossover is actually a cyclic two-point crossover in which CR, another control parameter of DE, determines how many consecutive genes of the mutated vector, x_G^must, on average are copied to the offspring x_G^child.

The selection scheme used in DE is also known as knockout competition. As the name suggests, DE plays a one-to-one competition between the principal parent, xⁱ_G, and its offspring x_G^child to select the survivor for the next generation. The DE selection scheme can be described as follows:

(16)

Repeating the above mentioned mutation and crossover operations on each individual of the current generation, DE creates a new generation of population which replaces the current generation. And this generation alternation process is iterated until the termination criteria is satisfied. The control parameters of DE (F, and P) are chosen beforehand and are kept constant throughout the search in this canonical version of the algorithm. The pseudo-code description of canonical DE is presented.

DE Algorithm

The estimation of the parameters of extended UNIQUAC was performed by minimizing the following objective function :

(17)

Where I, L, and N are the total number of experimentally determined, activity coefficient, VLE and excess enthalpy. In this investigation the data weighting factor, W_j=1, was used in all the parameter estimations. If, say for example, no VLE data was used in the parameter optimisation, the data weighting factor W_p was set to zero. The minimum objective function was minimized using differential evolution and Levenberg-Marquardt algorithms, where both results were compared.

During the parameters estimation, the quality of the fit was determined by the average absolute difference of the fit. The average absolute percentage difference were calculated from

(18)

(19)

(20)

Before the ability of the UNIQUAC model was tested to simultaneously correlate all the available aqueous MEA aqueous mixture data, the three individual types of data: activity coefficient, total pressure, and excess enthalpy were regressed individually to ascertain the scatter in the data. In all cases, any data sets that significantly deviated from the regression results were rejected. In this work a rather liberal tolerance of 20% deviation from the regression value was used to reject out-lying data points.

A complete list of the data sets used in this investigation is presented in the following tables (Tables 1 to 6).

Table 1: MEA-H₂O system data.

Table 2: Densities data for MEA-H₂O system.

Table 3: Ternary data for MEA-CO₂-H₂O system.

Table 4: Other sources of data used.

Table 5: UNIQUAC volume parameters ® and surface area parameters (q).

Table 6: UNIQUAC binary interaction parameters u_ij⁰ and u_ij^T .

The volume and surface area parameters, r and q respectively, for MEA, MEAH⁺ and MEA carbamate are determined by fitting to experimental data. The binary interaction parameters that are chosen to be adjusted are regressed to all types of experimental data in order to get a well-rounded model as well as to ensure accurate representation of different properties.

Figure 1 shows the effect of CO₂ concentration on the total pressure of the mixture. The model fits well with the experimental and especially for both cases studied, two extreme, 25^oC and 120^oC, temperatures and concentrations of carbon dioxide.

Figure 1: Model predictions and the experimental data (Lee et al. [38]) of the total pressure for CO₂-MEA-H₂O system at different MEA concentrations and temperatures.

Figure 2 presents a comparison fit between the DE algorithm and LM. The discrepancy between the two is quite noticeable at higher concentration of MEA. It is obvious when comparing the ADD% (γ_water) in the cases of DE and LM. ADD% (γ_water) for DE and LM are equal to 0.19 and 0.32, respectively.

Figure 2: Comparison of the predicted using UNIQUAC with DE and LM and the experimental water activity coefficients in (MEA+H₂O) mixture.

Figure 3 shows a good representation of the model by predicting the total pressure of the binary system (MEA+H₂O) at three different temperatures.

Figure 3: Prediction of the total pressure of the mixture (MEA+H₂O) at different temperatures. Nath & Bender [3] data (T=60, 78 and 91.7°C).

Presented in Figure 4, are the excess Gibbs energy and excess entropy calculated by the model (line) and LM approach (dots). The plots give noticeable discrepancies between the two approaches.

Figure 4: Predicted Excess Gibbs energy and excess Entropy of (MEA+H₂O) using UNIQUAC with DE and LM methods.

The extended UQIQUAC model has been successfully applied with Differential Evolution algorithm to calculate the interaction parameters and to correlate the experimental data on thermodynamic properties of MEA-H₂O-CO₂ system. The model has validated a lot of experimental data, but more research should be carried out to compare DE with other genetic algorithms in Multiple Objective Optimization. The model can be used to support process modeling and simulation of the CO₂ capture process with MEA.