## 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_{2} + 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 ${\overline{F}}_{2}^{\mathrm{s}}={\overline{F}}_{2}^{\mathrm{g}}$, and ${P}_{2}^{\mathrm{s}}{\varphi}_{2}^{\mathrm{s}}={P}_{2}^{\mathrm{g}}{\varphi}_{2}^{\mathrm{g}}$ is obtained when considering a pure substance in the solid phase. Here, it can be assumed that ${\varphi}_{2}^{\mathrm{s}}\approx 1$, and the solid volume is considered as pressure independent. Then, solid solute solubility (y_{2}) at equilibrium in a supercritical phase can be solved by the following equilibrium condition:

$${\mathit{y}}_{2}=\frac{{P}_{2}^{\mathrm{s}}{\varphi}^{\mathrm{s}}}{P{\varphi}^{\text{SCF}}}\text{exp}\left[\frac{{V}_{2}^{\mathrm{s}}\left(P-{P}_{2}^{\mathrm{s}}\right)}{RT}\right]$$ | (1) |

This equation can be obtained from the equifugacity between the solid component and the fluid phase; where ${P}_{2}^{\mathrm{s}}$ is the solid vapor pressure of the pure substance, ${V}_{2}^{\mathrm{s}}$ is the solid molar volume, T is the temperature, and ${\varphi}^{\mathrm{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]:

$$RT\phantom{\rule{0.25em}{0ex}}\text{ln}\phantom{\rule{0.25em}{0ex}}{\varphi}_{i}=\sum _{V}^{\infty}\left[{\left(\frac{\partial P}{\partial {n}_{i}}\right)}_{T,V,{n}_{j}}-\frac{RT}{V}\right]\mathrm{d}V-RT\phantom{\rule{0.25em}{0ex}}\text{ln}\phantom{\rule{0.25em}{0ex}}Z$$ | (2) |

Note that, the fugacity coefficient of the solid component can have considerably low values for these calculations (10^{−4} to 10^{−5}) [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{RT}{V-b}+\frac{a}{V\left(V+b\right)+b\left(V-b\right)}$$ | (3) |

$$a=0.457235\frac{{R}^{2}{T}_{\mathrm{c}}^{2}}{{P}_{\mathrm{c}}}\alpha \left({T}_{\mathrm{r}}\right)$$ | (4) |

$$b=0.077796\frac{R{T}_{\mathrm{c}}}{{P}_{\mathrm{c}}}$$ | (5) |

$$\alpha \left({T}_{\mathrm{r}}\right)={\left[1+\kappa \left(1-\sqrt{{T}_{\mathrm{r}}}\right)\right]}^{2}$$ | (6) |

$$\kappa =0.37646+1.54226\omega -0.26992{\omega}^{2}$$ | (7) |

_{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{RT}{V-{b}_{\mathrm{m}}}+\frac{{a}_{\mathrm{m}}}{V\left(V+{b}_{\mathrm{m}}\right)+{b}_{\mathrm{m}}\left(V-{b}_{\mathrm{m}}\right)}$$ | (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}_{\mathrm{m}}=\frac{\sum _{i}^{N}\sum _{j}^{N}{y}_{i}{y}_{j}{\left(b-\frac{a}{RT}\right)}_{ij}}{1-\sum _{i}^{N}\frac{{y}_{i}{a}_{i}}{{b}_{i}RT}+\frac{{A}_{\infty}^{E}\left(y\right)}{\text{\Omega}RT}}$$ | (9) |

$${\left(b-\frac{a}{RT}\right)}_{ij}=\frac{1}{2}\left({b}_{i}+{b}_{j}\right)-\frac{\sqrt{{a}_{i}{a}_{j}}}{RT}\left(1-{l}_{ij}\right)$$ | (10) |

$${a}_{\mathrm{m}}={b}_{\mathrm{m}}\left(\underset{i}{\sum ^{N}}\frac{{y}_{i}{a}_{i}}{{b}_{i}}+\frac{{A}_{\infty}^{E}\left(y\right)}{\text{\Omega}}\right)$$ | (11) |

_{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}\left(y\right)$ is calculated assuming that ${A}_{\infty}^{E}\left(y\right)\approx {A}_{0}^{E}\left(y\right)\approx {G}_{0}^{E}\left(y\right)$.

For a binary mixture:

$${b}_{\mathrm{m}}=\frac{{y}_{1}^{2}{\left(b-\frac{a}{RT}\right)}_{1}+2{y}_{1}{y}_{2}{\left(b-\frac{a}{RT}\right)}_{12}+{y}_{2}^{2}{\left(b-\frac{a}{RT}\right)}_{2}}{1-\frac{{y}_{1}{a}_{1}}{{b}_{1}RT}-\frac{{y}_{2}{a}_{2}}{{b}_{2}RT}+\frac{{G}_{0}^{E}\left(y\right)}{\text{\Omega}RT}}$$ | (12) |

$${\left(b-\frac{a}{RT}\right)}_{12}=\frac{1}{2}\left({b}_{1}+{b}_{2}\right)-\frac{\sqrt{{a}_{1}{a}_{2}}}{RT}\left(1-{l}_{12}\right)$$ | (13) |

$${a}_{\mathrm{m}}={b}_{\mathrm{m}}\left(\frac{{y}_{1}{a}_{1}}{{b}_{1}RT}+\frac{{y}_{2}{a}_{2}}{{b}_{2}RT}+\frac{{A}_{0}^{E}\left(y\right)}{\text{\Omega}}\right)$$ | (14) |

$$\frac{{G}_{0}^{E}}{RT}=\frac{\left(\frac{{A}_{12}}{RT}\right){y}_{1}{y}_{2}}{{y}_{1}\left(\frac{{A}_{12}}{{A}_{21}}\right)+{y}_{2}}$$ | (15) |

_{12}and A

_{21}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

_{12}, A

_{21}, and parameter l

_{12}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:

$$\text{min}\phantom{\rule{0.25em}{0ex}}f=\frac{100}{{N}_{D}}\sum _{i=1}^{{N}_{D}}{\left|\frac{{y}_{2}^{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}-{y}_{2}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{{y}_{2}^{\mathrm{e}\mathrm{x}\mathrm{p}}}\right|}_{i}$$ | (16) |

_{D}is the number of points in the experimental dataset and y

_{2}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:{\mathrm{\mathbb{R}}}^{d}\to \mathrm{\mathbb{R}}$, a population of particles P = {p_{1},⋯,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 ${\mathcal{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\left(l{b}_{i}^{t}\right)\le f\left(p{b}_{j}^{t}\right)\forall {p}_{i}\in {\mathcal{N}}_{i}$ [8].

For a given particle p_{i}, φ can be decomposed into ${\phi}_{k}=\phi /\left|{\mathcal{N}}_{i}\right|,\phantom{\rule{0.5em}{0ex}}\forall {p}_{k}\in {\mathcal{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 {\mathcal{N}}_{i}}{\phi}_{k}{U}_{i}^{t}\left(p{b}_{k}^{t}-{s}_{i}^{t}\right)$$ | (17) |

_{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) |

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}_{\mathrm{m}\mathrm{a}\mathrm{x}}-t}{w{t}_{\mathrm{m}\mathrm{a}\mathrm{x}}}\left({w}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{w}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+{w}_{\mathrm{m}\mathrm{i}\mathrm{n}}$$ | (19) |

_{max}marks the time at which ${w}^{t}={w}_{\mathrm{m}\mathrm{i}\mathrm{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.

**Table 1**

Scheme of the FPSO algorithm (taken from [8]).

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].

**Table 2**

Benchmark functions used.

Function | Mathematical formulation | Search space | min f |

Rastrigin | ${f}_{1}\left(x\right)=\sum _{i=1}^{D}\left[{x}_{i}^{2}-10\phantom{\rule{0.25em}{0ex}}\mathrm{cos}\left(2\mathrm{\pi}{x}_{i}\right)+10\right]$ | [−5.12;5.12]^{D} | 0 |

Griewank | ${f}_{2}\left(x\right)=\frac{1}{400}\sum _{i=1}^{D}{x}_{i}^{2}-\prod _{i=1}^{D}\text{cos}\left(\frac{{x}_{i}}{\sqrt{i}}\right)+1$ | [−600;600]^{D} | 0 |

Ackley | ${f}_{3}\left(x\right)=-20\text{exp}\left(-0.2\sqrt{\frac{1}{D}\sum _{i=1}^{n}{x}_{i}^{2}}\right)-\text{exp}\left(\frac{1}{D}\sum _{i=1}^{D}\text{cos}2\mathrm{\pi}{x}_{i}^{2}\right)+20+\mathrm{e}$ | [−32;32]^{D} | 0 |

**Table 3**

Configurations for all algorithms used on benchmark functions.

Algorithm | Parameters settings |

GA | Crossover probability = 0.8, mutation probability = 0.1 |

PSO | w_{max} = 0.9, w_{min} = 0.4, φ_{1} = φ_{2} = 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.

**Table 4**

Average and standard deviation of the best-of-run solution obtained for 200 runs for FPSO, PSO, and GA.

f | D | P | Gen | Average (standard deviation) | ||

PSO-ACO | GA | PSO | ||||

f_{1} | 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_{2} | 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_{3} | 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
_{12}[–0.1; 0.4], SS2 for A_{12}[0; 10], SS3 to A_{21}[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_{12}, A_{21}, and l_{12}) 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_{12}, A_{12}, A_{21}, 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
_{12}, A_{12}, A_{21}, 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
_{12}, A_{12}, A_{21}, P^{s}] is created.

## 4 Database used

Seven binary gas–solid phase mixtures of supercritical CO_{2} + 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.

**Table 5**

Details of GSE data used in this study.

Mix. | CO_{2}+ | Syst. | N_{D} | T(K) | ΔP (MPa) | $\text{\Delta}{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 |

**Table 6**

Thermodynamic properties of the substances involved in this study [23].

Substance | Formula | T_{c} (K) | P_{c} (MPa) | ω | V^{s} (m^{3}/kmol) |

Naphthalene | C_{10}H_{8} | 748.4 | 4.05 | 0.3020 | 0.1100 |

Biphenyl | C_{10}H_{12} | 773.0 | 3.38 | 0.4029 | 0.1312 |

2,3-Dimethylnaphthalene | C_{12}H_{12} | 785.0 | 3.22 | 0.4240 | 0.1547 |

2,6-Dimethylnaphthalene | C_{12}H_{12} | 777.0 | 3.17 | 0.4177 | 0.1547 |

Anthracene | C_{14}H_{10} | 873.0 | 2.90 | 0.4857 | 0.1430 |

Phenanthrene | C_{14}H_{10} | 869.0 | 2.90 | 0.4707 | 0.1528 |

Pyrene | C_{16}H_{10} | 936.0 | 2.61 | 0.5074 | 0.1591 |

Carbon dioxide | CO_{2} | 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_{2}.

## 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}^{\mathrm{s}}$, volume ${V}_{2}^{\mathrm{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_{12}, A_{12}, A_{21}, P^{s}) in this method.

Thus, an accurate modeling of the solubility is very important in order to obtain optimal values of l_{12}, A_{12}, A_{21}, and accurate value of P^{s}. The following deviations were considered to evaluate accuracy:

$$\left|\text{\%}\Delta {y}_{2}\right|=\frac{100}{{N}_{D}}\sum _{i=1}^{{N}_{D}}{\left|\frac{{y}_{2}^{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}-{y}_{2}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{{y}_{2}^{\mathrm{e}\mathrm{x}\mathrm{p}}}\right|}_{i}$$ | (20) |

$$\left|\text{\%}\Delta {P}^{\mathrm{s}}\right|=\frac{100}{{N}_{D}}{\sum _{i=1}^{{N}_{D}}\left|\frac{{P}^{\mathrm{s},\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}}-{P}^{\mathrm{s},\mathrm{e}\mathrm{x}\mathrm{p}}}{{P}^{\mathrm{s},\mathrm{e}\mathrm{x}\mathrm{p}}}\right|}_{i}$$ | (21) |

In these equations, N_{D} is the number of points in the experimental dataset, y_{2} 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_{12}, A_{12}, A_{21} using the combined method (PR-WS-VL + FPSO) for the systems used in this work. In addition, this table presents the deviations for y_{2} and P^{s} for each solid hydrocarbons calculated by Eqs. (20) and (21).

**Table 7**

Binary interaction parameters (l_{12}, A_{12}, A_{21}), and deviations for y_{2} and P^{s} obtained by the PR-WS-VL model + FPSO algorithm.

Syst. | l_{12} | A_{12} | A_{21} | $\left|\text{\%\Delta}{y}_{2}\right|$ | P^{s,calc}(Pa) | $\left|\text{\%\Delta}{P}^{\mathrm{s}}\right|$ |

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_{2} and P^{s}. Clearly, for y_{2} the deviations $\left|\text{\%\Delta}{y}_{2}\right|$ were lower than 6%, and deviations $\left|\text{\%\Delta}{P}^{\mathrm{s}}\right|$ were lower than 7.5% for the calculations of P^{s}.

The range for the interaction parameters (A_{12} and A_{21}) 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_{12} parameter of ${\left(b-\frac{a}{RT}\right)}_{12}$ (Eq. (13)) was defined as [–0.1; 0.4] with theoretical bases [6]. Fig. 1b shows the best values for the l_{12} parameters of the PR-WS-VL model.

Fig. 2 shows an example of the combined method performance (PR-WS-VL model + FPSO algorithm) on a binary mixture of supercritical CO_{2} + hydrocarbon. Fig. 2a shows the minimization of the objective function (Eq. (16)) as a function of the number of iterations for CO_{2} + Naphthalene (mixture 1, see Table 7). Fig. 2b shows a comparison between experimental values and calculated values of y_{2} 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).

A comparison was made between the results in the correlation of y_{2} 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% $\left(\left|\text{\%\Delta}{y}_{2}\right|>5\right)$ for the three models reported in the literature [25], and $\left|\text{\%\Delta}{y}_{2}\right|<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.

**Table 8**

Solubility deviations $\left|\text{\%\Delta}{y}_{2}\right|$ between three thermodynamic models reported in the literature [25] versus the PR-WS-VL model implemented in this study.

MixtureCO_{2}+ | 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 $\left|\text{\%\Delta}{P}^{\mathrm{s}}\right|$ 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].

**Table 9**

Deviations obtained for the estimation of sublimation pressure using: Coutsikos's GCM [26], and the combined method (PR-WS-VL + FPSO).

Compound | T (K) | P^{s,exp} (Pa)Ref. [23] | $\left|\text{\%\Delta}{P}^{\mathrm{s}}\right|$ Coutsikos's GCM [26] | $\left|\text{\%\Delta}{P}^{\mathrm{s}}\right|$ 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_{2} + 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_{2}; (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.