Modeling of gas–solid equilibrium of binary systems aided by Frankenstein PSO

Juan A. Lazzús

doi:10.1016/j.crci.2016.01.001

Résumé

Binary systems containing supercritical CO₂ + hydrocarbons were used for modeling gas–solid equilibrium by a combined method including a thermodynamic model and a meta-heuristic algorithm. The Peng–Robinson (PR) equation of state was used in the classical solubility equation. In addition, Wong–Sandler (WS) mixing rules were used, and the van Laar model (VL) was included in order to evaluate the excess Gibbs free energy that appears in these mixing rules. Then, a variant of particle swarm optimization (PSO), called Frankenstein PSO (FPSO) was implemented for minimizing the difference between calculated and experimental solubility values. The results showed that the FPSO algorithm is a very powerful tool for parameter estimation on the PR-WS-VL model with good performance and accuracy, and considerably low deviations. Therefore, values calculated by the combined method (PR-WS-VL + FPSO) are considered accurate enough for physical and engineering calculations, among other uses.

1 Introduction

The parameter estimation is highly significant in various scientific areas for the development of mathematical models that are based on parameters obtained from experimental data, which applies particularly to thermodynamic models [1]. The accurate prediction of physical properties on phase equilibrium is one of the most important applications for thermodynamic processes and it is where parameter estimation is required [2]. Particularly, an accurate prediction of phase equilibria of gas–solid mixtures is required for separation processes involving the so called supercritical fluid (SCF) [3].

Therefore, the best way to proceed is to adjust experimental data for a thermodynamic model that describe the gas–solid phase and to use the model obtained by the adjusted parameters for predicting the properties of these phases. However, the existing methods for solving the phase equilibrium only obtain local solutions. In several cases, optimum values for binary interaction parameters have been proven to be dependent on search interval and initial values [4]. Hence, it is necessary to develop new parameter optimization methods for creating global solutions and for determining the most suitable solution for a problem under a given set of restrictions associated with an objective function.

The current optimization problem lies, then, in an intelligent search issue where agents for determining an optimum value in a restricted search space can be used. This has led to modern optimization techniques, also known as meta-heuristic algorithms, that have caught great interest within the scientific community due to the great capacity of solving this kind of problem, and many other highly nonlinear problems [1].

In this work, twenty-two gas–solid binary systems containing supercritical CO₂ + hydrocarbons were used for modeling the gas–solid equilibrium using a combined method including a thermodynamic model and a meta-heuristic algorithm. The Peng–Robinson (PR) equation of state [5] was used in the classical solubility equation. Additionally, the Wong–Sandler (WS) mixing rules [6] were used, and the van Laar model (VL) [7] was applied for calculating the excess Gibbs free energy in the mixing rules. Then, a variant of particle swarm optimization (PSO), called Frankenstein PSO [8] was implemented in order to minimize the difference between calculated and experimental solubility values.

2 Thermodynamics of GSE

When using the fundamental equation of phase equilibria at a given pressure and temperature, the fugacity of the component in the gas phase is equal to the fugacity of the same component in the solid phase. Considering a subscript “2” for the solid component, the following equilibrium condition is obtained ${\bar{F}}_{2}^{s} = {\bar{F}}_{2}^{g}$ , and $P_{2}^{s} ϕ_{2}^{s} = P_{2}^{g} ϕ_{2}^{g}$ is obtained when considering a pure substance in the solid phase. Here, it can be assumed that $ϕ_{2}^{s} \approx 1$ , and the solid volume is considered as pressure independent. Then, solid solute solubility (y₂) at equilibrium in a supercritical phase can be solved by the following equilibrium condition:

y_{2} = \frac{P_{2}^{s} ϕ^{s}}{P ϕ^{SCF}} exp [\frac{V_{2}^{s} (P - P_{2}^{s})}{R T}]

(1)

This equation can be obtained from the equifugacity between the solid component and the fluid phase; where $P_{2}^{s}$ is the solid vapor pressure of the pure substance, $V_{2}^{s}$ is the solid molar volume, T is the temperature, and $ϕ^{s}$ is the fugacity coefficient of solid at the pressure P. In this case, the fugacity coefficient is calculated from standard thermodynamic relations, such as [7]:

R T ln ϕ_{i} = \sum_{V}^{\infty} [{(\frac{\partial P}{\partial n_{i}})}_{T, V, n_{j}} - \frac{R T}{V}] d V - R T ln Z

(2)

Note that, the fugacity coefficient of the solid component can have considerably low values for these calculations (10⁻⁴ to 10⁻⁵) [9].

When the fugacity coefficient is used in both phases, the method for the solution of the phase equilibrium problem is known as “the equation of state method”. In order to express the fugacity coefficient as a function of temperature-pressure-concentration an equation of state (EoS) with a set of mixing rules is needed [7]. In this context, cubic equations derived from van der Waals EoS are the most common EoS used for correlating phase equilibrium in mixtures at high and low pressures [9]; among these EoS, the Peng–Robinson equation of state (PR EoS) has proven to be capable of combining simplicity and accuracy required for the prediction and correlation of volumetric and thermodynamic properties of several compounds [3]. Commonly in GSE, van der Waals (vdW) mixing rules for both the energy and co-volume parameters are employed with one or two interaction parameters. However, an alternative to vdW mixing rules is the use of EoS/G^E mixing rules, which combine the advantages of successful cubic EoS and G^E models [9]. A good example of these alternative mixing rules are the Wong–Sandler mixing rules [6]. Here, the Peng–Robinson EoS [5] and Wong–Sandler mixing rules [6], are used as a thermodynamic model for calculating the fugacity coefficient ϕ.

The PR EoS can be formulated as follows [5]:

P = \frac{R T}{V - b} + \frac{a}{V (V + b) + b (V - b)}

(3)

with

a = 0.457235 \frac{R^{2} T_{c}^{2}}{P_{c}} α (T_{r})

(4)

b = 0.077796 \frac{R T_{c}}{P_{c}}

(5)

α (T_{r}) = {[1 + κ (1 - \sqrt{T_{r}})]}^{2}

(6)

κ = 0.37646 + 1.54226 ω - 0.26992 ω^{2}

(7)

where T_r = T/T_c is the reduced temperature. Note that, Peng–Robinson EoS is a two-parameter EoS (a and b), and it is completely predictive (in this form) once the critical temperature T_c, critical pressure P_c, and acentric factor ω, are given [5].

For a mixture:

P = \frac{R T}{V - b_{m}} + \frac{a_{m}}{V (V + b_{m}) + b_{m} (V - b_{m})}

(8)

where a_m and b_m are the EoS constants to be calculated using defined mixing rules. To this end, Wong–Sandler (WS) mixing rules are used here and can be expressed as follows: [6]:

b_{m} = \frac{\sum_{i}^{N} \sum_{j}^{N} y_{i} y_{j} {(b - \frac{a}{R T})}_{i j}}{1 - \sum_{i}^{N} \frac{y_{i} a_{i}}{b_{i} R T} + \frac{A_{\infty}^{E} (y)}{Ω R T}}

(9)

{(b - \frac{a}{R T})}_{i j} = \frac{1}{2} (b_{i} + b_{j}) - \frac{\sqrt{a_{i} a_{j}}}{R T} (1 - l_{i j})

(10)

a_{m} = b_{m} (\underset{i}{\sum^{N}} \frac{y_{i} a_{i}}{b_{i}} + \frac{A_{\infty}^{E} (y)}{Ω})

(11)

where a_m and b_m are the equation of state constants, l_ij is an adjustable parameter, Ω is equal to 0.34657, and

A_{\infty}^{E} (y)

is calculated assuming that

A_{\infty}^{E} (y) \approx A_{0}^{E} (y) \approx G_{0}^{E} (y)

For a binary mixture:

b_{m} = \frac{y_{1}^{2} {(b - \frac{a}{R T})}_{1} + 2 y_{1} y_{2} {(b - \frac{a}{R T})}_{12} + y_{2}^{2} {(b - \frac{a}{R T})}_{2}}{1 - \frac{y_{1} a_{1}}{b_{1} R T} - \frac{y_{2} a_{2}}{b_{2} R T} + \frac{G_{0}^{E} (y)}{Ω R T}}

(12)

{(b - \frac{a}{R T})}_{12} = \frac{1}{2} (b_{1} + b_{2}) - \frac{\sqrt{a_{1} a_{2}}}{R T} (1 - l_{12})

(13)

a_{m} = b_{m} (\frac{y_{1} a_{1}}{b_{1} R T} + \frac{y_{2} a_{2}}{b_{2} R T} + \frac{A_{0}^{E} (y)}{Ω})

(14)

where

G_{0}^{E} (y)

is calculated using the van Laar model (VL). Note that, the VL model has proven to be capable of good performance in high pressure phase equilibrium calculations [3]. In the case of a binary mixture, the VL model for

G_{0}^{E} (y)

is as follows [7]:

\frac{G_{0}^{E}}{R T} = \frac{(\frac{A_{12}}{R T}) y_{1} y_{2}}{y_{1} (\frac{A_{12}}{A_{21}}) + y_{2}}

(15)

where A₁₂ and A₂₁ are the binary interaction parameters into the van Laar G^E-model. At this point, the problem is reduced in order to determine the interaction parameters A₁₂, A₂₁, and parameter l₁₂ into the combined thermodynamic model (PR-WS-VL), from T–P–y data of GSE at a high pressure. Using an optimization algorithm, these optimal interaction parameters were calculated for minimizing the following objective function:

min f = \frac{100}{N_{D}} \sum_{i = 1}^{N_{D}} {| \frac{y_{2}^{c a l c} - y_{2}^{e x p}}{y_{2}^{e x p}} |}_{i}

(16)

where N_D is the number of points in the experimental dataset and y₂ is the mole fraction of the solid solute in the gas phase. The superscript denotes the experimental (exp) data point and the calculated (calc) values.

In this work, a variant of particle swarm optimization, called “Frankenstein PSO” is used for the minimization the objective function. Details of this meta-heuristic algorithm are present in the following section.

3 Frankenstein PSO

Frankenstein's particle swarm optimization (FPSO) was proposed by Montes de Oca and co-authors [8]. This variant of particle swarm optimization [10] consists of three main algorithmic components, such as: 1) a population topology reducing connectivity over time; 2) a fully informed particle swarm (FIPS) as a mechanism for updating the particle velocity, and; 3) a linearly decreasing inertia weight. These components are taken from adaptive hierarchical PSO [11], fully informed PSO [12], and the lineally decreasing inertia weight PSO [13], respectively.

In this meta-heuristic algorithm, each particle can move through a multidimensional search space, adjusting its position according to its own experience and the experience of neighboring in order to find an optimum solution [14]. Thus, in a d-dimensional objective function $f : ℝ^{d} \to ℝ$ , a population of particles P = {p₁,⋯,p_n} is randomly initialized in the solution space, and the performance of each particle is evaluated using a fitness function that defines the optimization problem [3, 15]. At any time step t, a particle p_i is associated with a position vector $s_{i}^{t}$ and a velocity vector $v_{i}^{t}$ , and a personal best vector $p b_{i}^{t}$ saves the best position that the particle has visited. If the particle p_i has a topological neighborhood $N_{i} \subseteq P$ of particles; then, the best personal best vector in a particle's neighborhood, called local best is a vector $l b_{i}^{t}$ such that $f (l b_{i}^{t}) \leq f (p b_{j}^{t}) \forall p_{i} \in N_{i}$ [8].

For a given particle p_i, φ can be decomposed into $φ_{k} = φ / | N_{i} |, \forall p_{k} \in N_{i}$ [12], and the velocity of each particle can be expressed as follows [8]:

v_{i}^{t + 1} = w^{t} v_{i}^{t} + \sum_{p_{k} \in N_{i}} φ_{k} U_{i}^{t} (p b_{k}^{t} - s_{i}^{t})

(17)

where w_t is the time-dependent inertia weight [13];

U_{i}^{t}

are d × d diagonal matrices with indiagonal elements distributed in the interval [0; 1) uniformly at random and generated at every iteration; φ_k are the acceleration coefficients [10]; φ (i.e., the sum of the acceleration coefficients) is equally distributed among all the neighbors of a particle [12];

s_{i}^{t}

is the current position of the particle, and

v_{i}^{t + 1}

is the new velocity at time t + 1;

p_{t}^{i}

is one of the the best solutions that this particle has reached [14].

The new particle position is computed by adding the velocity vector to the current position. Thus, the position of each particle at each generation is updated according to the following equation:

s_{i}^{t + 1} = s_{i}^{t} + v_{i}^{t + 1}

(18)

where

s_{i}^{t + 1}

is the new particle position at time t + 1 [10].

Variable w (inertia weight) in Eq. (17) is in charge of adjusting the equilibrium between local and global searches [13]. Note that, a high value of inertia weight implies a global search, while a low value leads to a local search [14, 11]. The following weighting function is used for Eq. (17):

w^{t} = \frac{w t_{m a x} - t}{w t_{m a x}} (w_{m a x} - w_{m i n}) + w_{m i n}

(19)

where wt_max marks the time at which

w^{t} = w_{m i n}

; w_min and w_max are the minimum and maximum values the inertia weight can take, respectively. Normally, wt_max is consistent with the maximum time allocated for the optimization process [8]. Note that, a decreasing w would counterbalance the exploratory behavior that the chosen topology change scheme could induce [15]. In FPSO, the time-varying topology starts as a fully connected one and decreases its connectivity until it ends up being a ring topology [16]. This topology transformation is performed in n−3 elimination steps [8] (please see [8, 15] for more details of the topology update). The pseudo-code for the FPSO algorithm is listed in Table 1.

The performance of the FPSO algorithm for global optimization is tested on 3 benchmark multimodal functions. The description of the functions used is listed in Table 2 (please see [17] for more details on benchmark functions). These functions are multimodal with a global minimum at the origin. Note that, finding the minimum for these functions is a very difficult problem due to its large number of local minima and its large search space. The search space values of x_i for these functions are listed in Table 2. Considering this difficulty, simulations were carried out to obtain a comparative performance analysis between the FPSO and other two meta-heuristic algorithms such as: basic particle swarm optimization (PSO) [10], and standard genetic algorithm (GA) [18]. From this analysis, the parameter settings of FPSO were carried out using this test. Table 3 shows the parameter configurations for all algorithms used on benchmark functions. These configurations are based on a trial-and-error procedure for FPSO, and based on common values used in the literature for PSO and GA [19].

Function	Mathematical formulation	Search space	min f
Rastrigin	$f_{1} (x) = \sum_{i = 1}^{D} [x_{i}^{2} - 10 \cos (2 π x_{i}) + 10]$	[−5.12;5.12]^D	0

Griewank	$f_{2} (x) = \frac{1}{400} \sum_{i = 1}^{D} x_{i}^{2} - \prod_{i = 1}^{D} cos (\frac{x_{i}}{\sqrt{i}}) + 1$	[−600;600]^D	0

Ackley	$f_{3} (x) = - 20 exp (- 0.2 \sqrt{\frac{1}{D} \sum_{i = 1}^{n} x_{i}^{2}}) - exp (\frac{1}{D} \sum_{i = 1}^{D} cos 2 π x_{i}^{2}) + 20 + e$	[−32;32]^D	0

Algorithm	Parameters settings
GA	Crossover probability = 0.8, mutation probability = 0.1
PSO	w_max = 0.9, w_min = 0.4, φ₁ = φ₂ = 2.0, V_max = ±S_max
FPSO	w_max = 0.9, w_min = 0.4, φ = 4.0, V_max = ±S_max

For benchmark functions, the mean squared error (MSE) between estimated and real responses for a number of given samples are considered as the fitness of estimated model parameters (min f). For a function of D-dimension, 200 separate runs were executed on each algorithm, and the standard deviation and the average best-of-run value were obtained. Several maximum generations Gen were used according to the complexity of the problem. In addition, three population sizes P were applied on each function, and the stopping criterion was set as reaching a fitness of 0.001. The number of generations run for each function was set to 1000, 1500 and 2000 corresponding to the dimensions 10, 20 and 30, and to populations 50, 100 and 150, respectively. Table 4 summarizes the results from this comparison. In addition this table shows the standard deviation and the mean deviation of the best-of-run solution. From the results present in Table 4, it is clear that FPSO performs better than PSO and GA. The significant improvement achieved by FPSO can be attributed to combining the best component of PSO variants.

f	D	P	Gen	Average (standard deviation)
PSO-ACO	GA	PSO
f₁	10	50	1000	0.0166 (0.0496)	4.3688 (3.1574)	2.3681 (1.9553)
20	100	1500	0.0309 (0.0281)	15.5721 (5.8548)	12.7331 (5.6188)
30	150	2000	0.0448 (0.3057)	36.5219 (10.0374)	34.1266 (4.3796)
f₂	10	50	1000	0.0298 (0.0321)	1.1416 (1.0157)	0.1012 (0.0253)
20	100	1500	0.0103 (0.0576)	1.0328 (1.2335)	0.2025 (0.2141)
30	150	2000	0.0087 (0.0145)	1.0934 (1.2786)	0.1679 (0.2264)
f₃	10	50	1000	0.0089 (0.1391)	1.3746 (2.0271)	0.3268 (1.0528)
20	100	1500	0.0144 (0.0242)	1.6542 (3.9597)	0.3999 (2.0988)
30	150	2000	0.0084 (0.0111)	2.9388 (3.2579)	1.0662 (2.1461)

The test results have demonstrated that the FPSO can be applied for parameter estimation with good accuracy, and the next step is to use the FPSO for estimating parameters on the thermodynamic model for GSE. Thus, the FPSO algorithm was implemented to optimize the PR-WS-VL model discussed above (see, Section 2).

The step-by-step to calculate the GSE using the FPSO is described below:

1) To start, the FPSO parameters must be defined (see Table 3). And the search space for the PR-WS-VL model is defined as SS1 to l₁₂ [–0.1; 0.4], SS2 for A₁₂ [0; 10], SS3 to A₂₁ [0; 10], and to SS4 −logP^s [1; 5]. Note that, in order to provide a substantial margin of safety, the range for the interaction parameters (A₁₂, A₂₁, and l₁₂) was defined by thermodynamic considerations [3].
2) The initial population is generated with P = 50 particles associated with the positions s_i. Where SVect and VVect save a particle position and its corresponding velocity in the search space, respectively. These vectors contain the possible values of the thermodynamic parameters: l₁₂, A₁₂, A₂₁, and P^s.
3) The vector PVect is created. This vector contains the best current position of particle.
4) The objective function is evaluated for each parameter set (l₁₂, A₁₂, A₂₁, P^s). The evaluation is carried out by calculating the solubility by using the PR-WS-VL model for each value of the parameter set. Thus the vector containing the deviations between experimental and calculated values of solubility ObjFun is generated.
5) If the evaluation of the objective function does not meet the suggested criteria (Eq. (16)), then the vectors SVect and VVect are updated. Now the new position of the particles is calculated by adding the velocity to the current position (Eqs. (17) and (18)). This process is repeated until the criterion condition is satisfied (Eq. (16)), or until the maximum number of iterations is exceeded t_max = 100.
6) To complete the process, the vector containing the optimal values of the thermodynamic parameters Sol = [l₁₂, A₁₂, A₂₁, P^s] is created.

4 Database used

Seven binary gas–solid phase mixtures of supercritical CO₂ + hydrocarbon containing: naphthalene [20], biphenyl [20], 2,3-dimethylnaphthalene [21], 2,6-dimethylnaphthalene [21], anthracene [22], phenanthrene [22], and pyrene [22] are considered in this work. The experimental GSE data of 22 systems taken from the literature [20–22] are presented in Table 5. For these systems, pressure and temperature cover wide ranges from 8 to 45 MPa and from 308 to 338 K, respectively. Table 6 shows the thermodynamic properties (V^s, T_c, P_c, and ω) of the substances used [23]. Note that, these substances and their mixtures are of interest for several applications of chemical industries, and other uses.

Mix.	CO₂+	Syst.	N_D	T(K)	ΔP (MPa)	$Δ y_{2} \times 10^{4}$
1	Naphthalene [20]	1	8	308.15	9–26	75–192
2	16	328.15	8–29	13–538
3	19	333.55	10–29	52–980
4	7	338.05	15–23	247–790
2	Biphenyl [20]	5	8	308.95	10–44	104–154
6	8	318.55	15–45	156–272
7	8	322.65	15–48	178–361
3	2,3-Dimethylnaphthalene [21]	8	5	308.15	10–28	22–64
9	5	318.15	10–28	12–72
10	5	328.15	10–28	3–90
4	2,6-Dimethylnaphthalene [21]	11	5	308.15	10–28	19–45
12	5	318.15	10–28	7–68
13	5	328.15	10–28	3–92
5	Anthracene [22]	14	30	313.15	11–34	0.2–1
15	30	323.15	12–35	0.2–1.5
16	30	333.15	13–36	0.2–2
6	Phenanthrene [22]	17	26	313.15	10–35	4–24
18	30	323.15	11–34	2–31
19	30	333.15	11–35	2–41
7	Pyrene [22]	20	32	313.15	10–35	0.5–3.6
21	32	323.15	11–35	0.4–4.5
22	32	333.15	11–35	0.2–5.5

Substance	Formula	T_c (K)	P_c (MPa)	ω	V^s (m³/kmol)
Naphthalene	C₁₀H₈	748.4	4.05	0.3020	0.1100
Biphenyl	C₁₀H₁₂	773.0	3.38	0.4029	0.1312
2,3-Dimethylnaphthalene	C₁₂H₁₂	785.0	3.22	0.4240	0.1547
2,6-Dimethylnaphthalene	C₁₂H₁₂	777.0	3.17	0.4177	0.1547
Anthracene	C₁₄H₁₀	873.0	2.90	0.4857	0.1430
Phenanthrene	C₁₄H₁₀	869.0	2.90	0.4707	0.1528
Pyrene	C₁₆H₁₀	936.0	2.61	0.5074	0.1591
Carbon dioxide	CO₂	304.2	7.38	0.2236	–

In addition, the binary system data were selected only if they met two conditions: (i) GSE data must correspond to “true” gas–solid data. All data obtained from theoretical methods, correlations, or extrapolations of any kind were not considered; and (ii) GSE data must be accurate enough and pass a test of thermodynamic consistency. This method was based on the Gibbs–Duhem equation in combination with the PR-WS-VL model [24].

On these systems, the combined method (PR-WS-VL model + FPSO algorithm) was used, for considering the deviations between experimental and calculated values of solid hydrocarbon solubility in the high-pressure supercritical CO₂.

5 Results and discussion

The fugacity coefficient is calculated by the PR-WS-VL model (see, Eq. (16)). The other properties of the solid component (sublimation pressure $P_{2}^{s}$ , volume $V_{2}^{s}$ ) can be obtained from third party information or independent data (see, Eq. (1)). Note that, the sublimation pressure values are usually small, and the experimental techniques cannot be used to obtain accurate values in many cases [1]. Considering the valuable information that may be derived either directly or indirectly from sublimation data, it is rather surprising that there is so little quantitative information available in the literature on the sublimation process [1,3]. For this reason, P^s was considered as a fourth parameter to estimate in addition to binary interaction parameters (l₁₂, A₁₂, A₂₁, P^s) in this method.

Thus, an accurate modeling of the solubility is very important in order to obtain optimal values of l₁₂, A₁₂, A₂₁, and accurate value of P^s. The following deviations were considered to evaluate accuracy:

| % Δ y_{2} | = \frac{100}{N_{D}} \sum_{i = 1}^{N_{D}} {| \frac{y_{2}^{c a l c} - y_{2}^{e x p}}{y_{2}^{e x p}} |}_{i}

(20)

| % Δ P^{s} | = \frac{100}{N_{D}} {\sum_{i = 1}^{N_{D}} | \frac{P^{s, c a l c} - P^{s, e x p}}{P^{s, e x p}} |}_{i}

(21)

In these equations, N_D is the number of points in the experimental dataset, y₂ is the mole fraction of the solid solute in the gas phase, P^s is the sublimation pressure, superscripts exp and calc denote the experimental data point and calculated values, respectively.

Table 7 shows the calculated parameters l₁₂, A₁₂, A₂₁ using the combined method (PR-WS-VL + FPSO) for the systems used in this work. In addition, this table presents the deviations for y₂ and P^s for each solid hydrocarbons calculated by Eqs. (20) and (21).

Syst.	l₁₂	A₁₂	A₂₁	$\| %Δ y_{2} \|$	P^s,calc(Pa)	$\| %Δ P^{s} \|$
1	0.0716	2.7873	4.0978	3.03	28.399	0.58
2	0.0368	3.1752	0.6968	5.70	156.610	0.02
3	0.0923	2.3283	0.7834	5.80	240.120	0.15
4	0.0825	2.5742	0.8335	3.78	337.000	0.33
5	0.0651	3.1119	3.9742	5.03	4.175	0.17
6	0.0872	2.5289	2.3889	0.60	11.000	0.42
7	0.0943	2.3595	1.6574	1.62	16.500	0.41
8	0.1104	3.4514	2.3094	2.11	1.250	1.34
9	0.0665	4.1733	2.7115	5.02	3.477	0.06
10	0.0650	4.0938	1.7899	4.97	8.941	0.51
11	0.1001	4.4171	2.8696	4.98	1.213	5.57
12	0.0803	4.3035	2.1979	2.30	3.458	5.46
13	0.0770	4.2435	3.1872	5.36	9.150	3.94
14	0.2059	5.5907	0.1270	2.97	0.009	0.00
15	0.1852	5.2431	0.5796	1.56	0.026	0.00
16	0.1575	5.5070	1.9611	4.92	0.075	0.00
17	0.1585	5.0025	0.2708	2.26	0.110	0.00
18	0.1435	4.6240	0.4022	4.77	0.297	7.48
19	0.1373	5.1451	0.2493	5.53	0.944	0.11
20	0.2473	6.9176	2.1123	2.01	0.010	0.00
21	0.2090	6.3852	1.7488	3.48	0.020	0.00
22	0.1798	6.1121	1.4004	3.20	0.048	2.04

The results found using the combined method (PR-WS-VL model + FPSO algorithm) show low deviations between experimental and calculated values for y₂ and P^s. Clearly, for y₂ the deviations $| %Δ y_{2} |$ were lower than 6%, and deviations $| %Δ P^{s} |$ were lower than 7.5% for the calculations of P^s.

The range for the interaction parameters (A₁₂ and A₂₁) was defined as [0; 10] in order to provide a substantial margin of safety. This range was based on thermodynamic considerations [3], it will extremely likely contain the optimum parameter values. Fig. 1a shows the binary interaction parameters calculated using the FPSO algorithm and based on the minimization of Eq. (16). Also, the range for the l₁₂ parameter of ${(b - \frac{a}{R T})}_{12}$ (Eq. (13)) was defined as [–0.1; 0.4] with theoretical bases [6]. Fig. 1b shows the best values for the l₁₂ parameters of the PR-WS-VL model.

Fig. 1
(a) Calculated values of binary interaction parameters (A₁₂ and A₂₁) by minimization of the objective function using the FPSO algorithm. (b) Calculated values for the k₁₂ parameter by FPSO algorithm.

Fig. 2 shows an example of the combined method performance (PR-WS-VL model + FPSO algorithm) on a binary mixture of supercritical CO₂ + hydrocarbon. Fig. 2a shows the minimization of the objective function (Eq. (16)) as a function of the number of iterations for CO₂ + Naphthalene (mixture 1, see Table 7). Fig. 2b shows a comparison between experimental values and calculated values of y₂ for these systems. Fig. 2c shows the comparison between experimental values and calculated values of P^s. This figure ratifies the capabilities of the algorithm presented above (see, Table 7).

Fig. 2
(a) Example of convergence by minimization of the objective function for CO₂ + naphthalene. (b) Experimental (dots) and calculated values (dashed line) of solubility y₂ as a function of pressure for CO₂ + naphthalene (systems 1 to 4 [20]). (c) Experimental (dots) and calculated values (dashed line) for the P^s as a function of temperature for naphthalene.

A comparison was made between the results in the correlation of y₂ by the PR-WS-VL model and the results obtained from three other thermodynamic models reported in the literature [25] (for clarity, FPSO is not used on these other models). This comparison involved the following models: Peng–Robinson equation of state and van der Waals mixing rules (PR-vdW); Peng–Robinson equation of state combined with Wong–Sandler mixing rules and Universal-Quasi-Chemical theory G^E-model (PR-WS-UNIQUAC); Peng–Robinson equation of state combined with Wong–Sandler mixing rules and Non-Randon-Two-Liquid G^E-model (PR-WS-NRTL); and the PR-WS-VL model used in this study. Table 8 shows the results of this comparison applied on four mixtures with common datasets [20,21,25]. This table shows solubility deviations higher than 5% $(| %Δ y_{2} | > 5)$ for the three models reported in the literature [25], and $| %Δ y_{2} | < 5$ for the PR-WS-VL model used in this study. Particularly, the average deviations for the four mixtures were: 19.33% for the PR-vdW model, 7.60% for PR-WS-UNIQUAC, 7.69% for PR-WS-NRTL, and 4.25% for the PR-WS-VL model (for interested readers, all binary interaction parameter values of these models can be obtained from the original study [25]). From these results, the accuracy of the results obtained using the PR-WS-VL model is better than the accuracy obtained from the other models with binary interaction parameters. It has been shown that the combined method (PR-WS-VL + FPSO) is able to give a good representation of the solubility of a solid component in a supercritical fluid.

MixtureCO₂+	PR-vdW Ref. [25]	PR-WS-UNIQUAC Ref. [25]	PR-WS-NRTL Ref. [25]	PR-WS-VL This work
Naphthalene	34.6	8.2	9.02	4.58
Phenanthrene	20.3	7.4	8.06	4.19
2,3-Dimethylnaphthalene	13.5	7.2	6.42	4.03
2,6-Dimethylnaphthalene	8.9	7.6	7.25	4.21

In addition, a comparison was made between the results obtained in the estimation of P^s by the combined method (PR-WS-VL + FPSO) and the results obtained by a group contribution method (GCM) to predict the sublimation pressure [26]. Table 9 shows this comparison and shows that the deviation for the sublimation pressures $| %Δ P^{s} |$ obtained from Coutsikos's GCM [26] were greater than that obtained from the combined method (PR-WS-VL + FPSO). These results show a tremendous increase in accuracy for estimating P^s. Note that, the P^s estimation may be considered as a difficult problem in different physico-chemical processes, since nearly no experimental determinations available under 0.1 (Pa) [26]. Also, consider that of all properties involved in the solubility calculation at high pressure, P^s has received low attention in the literature although the latter considered to be directly related to solubility [3].

Compound	T (K)	P^s,exp (Pa) Ref. [23]	$\| %Δ P^{s} \|$ Coutsikos's GCM [26]	$\| %Δ P^{s} \|$ This work
Naphthalene	308.15	28.236	17.23	0.58
328.15	156.637	20.12	0.02
333.55	239.762	21.35	0.15
338.05	338.122	23.63	0.33
Biphenyl	308.95	4.168	13.87	0.17
318.55	11.046	14.37	0.42
322.65	16.433	15.18	0.41
2,3-Dimethylnaphthalene	308.15	1.267	21.12	1.34
318.15	3.479	18.39	0.06
328.15	8.987	16.61	0.51
2,6-Dimethylnaphthalene	308.15	1.149	25.50	5.57
318.15	3.279	24.63	5.46
328.15	8.803	23.00	3.94
Anthracene	313.15	0.009	6.34	0.00
323.15	0.026	7.15	0.00
333.15	0.075	7.57	0.00
Phenanthrene	313.15	0.110	22.61	0.00
323.15	0.321	21.88	7.48
333.15	0.945	21.00	0.11
Pyrene	313.15	0.010	10.71	0.00
323.15	0.020	11.12	0.00
333.15	0.049	11.89	2.04

All these results showed that the FPSO algorithm is a very powerful tool for parameter estimation on the PR-WS-VL model and has a good performance and accuracy, and considerably low deviations. Therefore, values calculated using the combined method PR-WS-VL model + FPSO algorithm are considered to be accurate enough for physico-chemical calculations, among other uses. In addition, the FPSO algorithm can be applied to optimize other thermodynamic models with good accuracy and performance.

6 Conclusions

Based on the Results and Discussion Section (see, Section 5), the following main conclusions were drawn: (i) the thermodynamic model PR-WS-VL is appropriate to describe the gas–solid equilibrium of binary systems (supercritical CO₂ + hydrocarbon) at high pressure; (ii) Frankenstein PSO algorithm is a good tool to minimize the difference between calculated and experimental solubility of a solid solute y₂; (iii) the results showed that the Frankenstein PSO algorithm is a very powerful tool for estimating binary interaction parameters with low deviations; (iv) the values calculated using the complete method (PR-WS-VL model + FPSO) are considered to be accurate enough for physico-chemical calculations, and other uses; and (v) Frankenstein PSO algorithm can be applied in order to optimize other thermodynamic models with good accuracy and performance.

Acknowledgments

The author thanks the Direction of Research of the University of La Serena (DIULS) and the Department of Physics of the University of La Serena (DFULS) for the special support that made possible the preparation of this paper.