A Practical Fitting Method Involving a Trade-Off Decision in the Parametrization Procedure of a Thermodynamic Model and Its Repercussion on Distillation Processes

Abstract

The design of processes containing information on phase equilibria must be carried out through a series of steps, experimentation ↔ verification ↔ modeling ↔ simulation. Each of these steps should be rigorously performed to guarantee a good representation of the behavior of the system under study, whose adequate modeling could be used to simulate the corresponding process. To carry out the different previous tasks, two representative systems, extracted from known database, are used. The quality checking of experimental data series is certified through several thermodynamic consistency methods. The modeling is done by applying a multi-objective optimization procedure, which allows to define a solution front (Pareto front) for different sub-models that are established in this work. The fitness of trade-off solutions, obtained from the efficient front, on the design of distillation processes is analyzed through a simulation.

1. Introduction

Nowadays, process simulation [1] plays an important role in the chemical engineering field as an indispensable tool to gain precise knowledge about the process units. The use of powerful mathematical-computational tools allows to accomplish an optimal final design of chemical processes. An important matter is to ensure that the mathematical model correctly represents the quantities that reflect the state of the studied system. Then, what is a model? A model is a mathematical relationship that links the state variables of a system, such as temperature, pressure, or compositions, influencing the process performance. Therefore, the accuracy of the selected model is essential and greatly affects the final results of the simulation and design processes. Two milestones should be considered.

The first one is the selection of the model, built-in with the mathematical relationships that best represent the variables of the system. There is no a priori a procedure to make this choice, so heuristic or experience-based criteria are generally used [2]. The second question refers to get the best parameters that complete the definition of the model for a given data series. For this latter case, there are several numerical procedures [3, 4] that allow to address the problem to optimize the parameter set considering the starting hypothesis.

The thermodynamic properties having the greatest influence on the simulation of separation processes are those related to phase equilibria (vapor-liquid equilibrium, VLE, in the scope of this work), as well as other thermodynamic quantities that arise in the mixing process. These properties are associated with the excess Gibbs energy g^E, which is written as g^E = g^E(x_i,p,T). As such, the goal of the modeling is to achieve a functional type of f = f(θ,x_i,p,T) that minimizes the norm |g^E − f|. The vector θ represents a set of parameters in the model, which must be optimized. However, due to the existing relation between phase and mixing properties, former approach may not be enough. In the first place, g^E values can only be obtained from VLE, which satisfies a dependency of the type, F(x_i,p,T) = 0, preventing us from obtaining individual relations of g^E with each one of the variables. On the other hand, defects in the experimental data could give place to incoherencies between experimental activity coefficients γ_i = γ_i(x_i,y_i,p,T) and the excess Gibbs function g^E. Thus, the VLE fitting process becomes a bi-objective optimization problem; hence, two error functions are included in the correlation problem. The complexity grows as other properties, such as h^E = h^E(x,p,T), are included in the modeling, since new error functions should also be minimized. However, one of the benefits of this approach is to use a unique thermodynamic model that avoids the issues caused by possible discontinuities or inconsistencies between partial models of some properties. Even more relevant is that this allows the model to describe better the physics of the system under study. This supposes a mean to verify the coherence of the mathematical formalism imposed by thermodynamics, validating the different properties. A drawback of this practice is the increase of the complexity of the procedure, but this is just a numerical issue of relative complexity. Therefore, addressing the global modeling of the thermodynamic behavior of solutions as a problem with multiple objectives is a notable contribution in chemical engineering.

It is known that the resolution of multi-objective problems does not produce a single result; on the contrary, a set of non-dominated results that constitute the so-called Pareto front [5] is obtained. There is no precise mathematical criterion that allows the selection of a unique result. However, the process simulation requires a single result to define the intended design, having to resort to external criteria different from those used to obtain the front.

This study evaluates the effect of choosing the different results from the Pareto front in the simulation task. To achieve this, some partial goals are proposed such as (a) to establish a rigorous methodology to carry out the optimization procedure with the suggested modeling and (b) to check the real impact of the chosen model on the simulation, with the purpose of proposing a selection criterion. Thus, the designed methodology should include different stages, like the data selection used in the procedure, obtaining the result front and the election of the final result. Two systems, considered as standard in many studies on thermodynamic behavior of solutions, are selected since the necessary experimental information (VLE and h^E) is available in literature. After checking the data sets [6, 7], making up the testing ↔ modeling steps, and having obtained the corresponding result fronts, each candidate model is analyzed on its suitability for process simulation.

2. Modeling procedures

2.1 Thermodynamic model

The correlation of the iso-p {VLE + h^E} data for each one of the two systems selected is carried out using a parametric equation already used by us [8], which is applied to the excess Gibbs function of a binary system:

$$g^{\mathrm{E}}\left(p,T,\mathbf{x}\right)=z_1\left(1-z_1\right)\sum _{\mathrm{i}=0}^2{g}_{\mathrm{i}}\left(p,T\right)z_1^{\mathrm{i}}{z}_1=\frac{x_1}{x_1+k_{\mathrm{g}}^{2‐1}{x}_2}$$ E1

and as described in [8]:

$$g_{\mathrm{i}}\left(p,T\right)=g_{\mathrm{i}1}+g_{\mathrm{i}2}{p}^2+g_{\mathrm{i}3}\mathit{pT}+g_{\mathrm{i}4}/T+g_{\mathrm{i}5}{T}^2$$ E2

where k_g^2–1 is a fitting parameter. In this study, several forms of Eq. (2) are tested, by adapting it according to the availability of experimental data. Thus, different “sub-models” are defined by neglecting terms in Eq. (2) which give rise to the following four cases:

$$\text{M1}\to \left\{g_{\mathrm{i}2},g_{\mathrm{i}3},g_{\mathrm{i}4},g_{\mathrm{i}5}\right\}=0;\text{M2}\to \left\{g_{\mathrm{i}2},g_{\mathrm{i}3},g_{\mathrm{i}5}\right\}=0;\text{M3}\to \left\{g_{\mathrm{i}2},g_{\mathrm{i}5}\right\}=0;\text{M4}\to g_{\mathrm{i}2}=0$$ E3

From Eq. (1), the expression that characterizes the activity coefficients is

$$\mathit{RT}ln{\gamma }_{\mathrm{i}}=z_1\left(1-z_1\right)\left(g_0+g_1{z}_1+g_2{z}_1^2\right)+\left(1-\mathrm{i}-x_1\right)\left[g_0+2\left(g_1-g_0\right)z_1+3\left(g_2-g_1\right)z_1^2-4g_2{z}_1^3\right]k_{\mathrm{g}}^{2‐1}{\left(z_1/x_1\right)}^2$$ E4

Likewise, the excess enthalpies h^E are calculated considering.

$$h^{\mathrm{E}}=g^{\mathrm{E}}-T{\left(\partial g^{\mathrm{E}}/\partial T\right)}_{\mathrm{x},\mathrm{p}};h^{\mathrm{E}}=z_1\left(1-z_1\right)\left(h_0+h_1{z}_1+h_2{z}_1^2\right)$$ E5

Eqs. (1) and (5) assume that the different sub-models established by Eq. (3) impose certain behavior hypothesis in the mixture. Thus, sub-model M1 implies that g^E = h^E. For the remaining assumptions, g^E ≠ h^E, allowing different functions to express the h^Es.

2.2 Calculation of the efficient fronts

To obtain the optimal parameters of the different sub-models that represent the data set iso-p VLE and h^E, the multi-objective optimization (MOO) procedure is used, which is characterized by the vector

$$\mathbf{\text{OF}}=\left[s\left(g^{\mathrm{E}}/\mathit{RT}\right)s\left(h^{\mathrm{E}}\right)\right]$$ E6

where s(y^E) is the root of the mean square error (or RMSE) calculated by the model when representing the generic property y^E:

$$s\left(y^{\mathrm{E}}\right)={\left[\sum {\left(y_{\exp }^{\mathrm{E}}-y_{\mathrm{cal}}^{\mathrm{E}}\right)}^2/m\right]}^{0.5};m=\text{number of experimental points}$$ E7

The result of such a problem is a discrete set of values, taken from the efficient front, s(g^E/RT) = f[s(h^E)], where each one of these give rise to different estimates of the thermodynamic behavior. Here the ε-constraint algorithm [3] is used to solve the MOO optimization problem. This procedure converts the multi-objective nonlinear problem (MNLP) into multiple single-objective problem (constrained nonlinear problem (CNLP)), thence minimizing only the first objective in Eq. (6), that is,

$${\mathit{\text{OF}}}_{\epsilon }=s\left(g^{\mathrm{E}}/\mathit{RT}\right)$$ E8

using the other vector element to establish a constraint in the calculation, that is,

$$\mathit{\text{constraint}}\to s\left(h^{\mathrm{E}}\right)<\epsilon$$ E9

The value of ε is limited, from 0 to 500 J mol⁻¹, since greater errors in the h^E are not acceptable. The efficient front is then achieved by solving the CNLP (Eqs. (8) and (9)) for different values of ε in the indicated range. To do so, we used a hybrid evolutionary algorithm [3], coupled with Nelder-Mead algorithm [9], for local refinement.

3. Verification and selection of data series

3.1 Acetone + ethanol

There are 10 useful references in literature containing VLE data for the acetone(1) + ethanol(2) system [10, 11, 12, 13, 14, 15, 16, 17, 18, 19] showing a total of 15 experimental series. In Figure 1, seven sets are iso-p (Figure 1(a)), and eight sets are iso-T (Figure 1(b)).

Figure 1.

Plot of VLE data of the binary: acetone(1) + ethanol(2), iso-p ≈ 101.32 kPa. (a) T,x,y; (c) γ,x; (d) g^E/RT,x. (×) Ref. [10]; (∇) Ref. [11]; (○) Ref. [12]; (◊) Ref. [13]; (+) Ref. [14]; (△) Ref. [15]; (□) Ref. [16]; (b) p,x,y. iso-T: (▷) Ref. [14]; (○) Ref. [17]; Ref. [18] (×) T = 372.67 K; (▽) T = 422.56 K, (△) T = 397.67 K; Ref. [19], (□) T = 344.19 K, (+) T = 363.19 K, (◊) T = 358.18 K.

For use we refer to data series with information for the two phases. This is a requirement to check the thermodynamic consistency. Figure 1(b) shows an azeotrope at high pressures, which moves from the acetone-rich region toward pure ethanol as pressure increases; the separation between the compositions in the two phases is reduced. Figure 1(c) and (d) contains, respectively, the variation of γ_i and g^E/RT with the liquid composition for the system at 101 kPa. Ethanol’s activity coefficient shows some errors that surely affect the global consistency of the data series. High slopes are observed as x₁ → 0 along with negative values of the natural logarithm of γ₁ for data series other than those referenced [14, 15, 16]. This is a clear sign of inconsistency.

Table 1 shows the results obtained in the application of different consistency methods. The observation of the said table produces some comments which are interesting in the work development. The Wisniak test and the direct van Ness test accept most of the experimental series. The first of the methods rejects two (noted as n° 10,12) while the second fails to validate only one, n° 8. The Kojima test rejects five data series, with the error observed in the system n° 12 being critical. A new methodology recently proposed by us [6, 7] was also applied, which rejects nine systems: three by the integral form (n° 9, 11, 12) and six by the differential form (n° 1, 3, 4, 8, 9, 12). These systems present obvious signs of inconsistency that are also observed by at least one of the other methods. Therefore, the overall assessment of these series (in other words, the quality of their experimental data) is not positive, ruling out the use of these series for other subsequent tasks, such as modeling and, specially, simulation.

N°	Ref.	Type	Area		Fredenslund		Wisniak		Kojima		Direct-Van Ness		Proposed test
			d%	V	×100	V	Dw	V	M(Ii)	V	δln(γ1/γ2)	V	f-int.	V	f-dif.	V	VF
1	10	Iso-101 kPa	15	nv	1.9	nv	2.5	v	24	v	0.09	v	0.06	v	−0.06	nv	nv
2	11	Iso-101 kPa	13	nv	1.8	nv	2.6	v	18	v	0.09	v	0	v	0.004	v	v
3	12	Iso-101 kPa	15	nv	1.8	nv	1.7	v	10	v	0.09	v	0.14	v	−0.05	nv	nv
4	13	Iso-101 kPa	23	nv	2.2	nv	2	v	15	v	0.09	v	0.11	v	−0.08	nv	nv
5	14	Iso-101 kPa	5	nv	1.2	nv	0.5	v	27	v	0.04	v	0	v	0.003	v	v
6	15	Iso-101 kPa	7.1	nv	0.5	v	2.1	v	16	v	0.02	v	0.16	v	0.005	v	v
7	16	Iso-101 kPa	1.4	v	0.6	v	2.8	v	40	nv	0	v	0.02	v	0.006	v	v
8	14	Iso-328 K	1.7	v	0.9	v	2.1	v	74	nv	0.01	v	0.36	v	−0.050	nv	nv
9	17	Iso-298 K	41	nv	2.7	nv	1.5	v	82	nv	0.18	nv	−0.09	nv	−0.370	nv	nv
10	18	Iso-372 K	3.3	nv	0.4	v	3.3	nv	14	v	0.02	v	0.11	v	0.009	v	v
11	18	Iso-398 K	0.4	v	0.6	v	3	v	69	nv	0.02	v	−2	nv	0.001	v	nv
12	18	Iso-423 K	12	nv	0.6	v	4.4	nv	159	nv	0.01	v	−24.5	nv	−0.01	nv	nv
13	19	Iso-344 K	2.8	nv	0.3	v	2.1	v	14	v	0.01	v	0.42	v	0.016	v	v
14	19	Iso-358 K	3.3	nv	0.3	v	3	v	9	v	0.01	v	0.64	v	0.014	v	v
15	19	Iso-363 K	0.2	v	0.2	v	2.5	v	19	v	0.01	v	0.83	v	0.029	v	v

Table 1.

Values obtained in the application of the thermodynamic consistency-test to VLE data of acetone + ethanol at different conditions.

Limiting values: areas-test: d < 2; Fredenslund-test: δy×100 < 1; Wisniak-test: D_w < 3; Kojima-test: M(I_i) < 30; Direct Van Ness-test: δ(ln γ₁/ln γ₂) < 0.16; proposed-test: f-int > 0, f-dif > 0. Header notation: V → verified; FV → fully verified. Test outcome: v → verified; nv → not verified; nd/− → not available.

Regarding h^E, six references were found [20, 21, 22, 23, 24, 25]. All extracted values are shown in Figure 2(a and b), where the different series show an acceptable coherence, having (dh^E/dT)_p,x > 0.

Figure 2.

h^E values for the binary acetone(1) + ethanol(2). (a) p = 101.32 kPa; (b) p = 400 kPa. (○) Ref. [20]; (×) Ref. [21]; (□) Ref. [22]; (+) Ref. [23]; (△) Ref. [24]; (◊) Ref. [25].

3.2 Benzene + hexane

In studies related to the thermodynamics of solutions, the experimental information generated by the binary benzene + hexane, for different properties, has been, and still is, used by researchers in that area as a reference of their investigations; therefore, the choice of this system is justified. Twenty-three useful VLE data series were found in the bibliography [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42], 16 iso-p and 7 iso-T. Figure 3 shows important errors in some of the series [29, 31, 37] at 101 kPa, with data missing the trend observed in the other series. These are propagated to the T vs x,y representations (Figure 4(a)). The random error for these series is serious as evidenced in Figure 4(b).

Figure 3.

Plot of iso-p VLE data of the binary benzene + hexane at 101 kPa. (a) T,x,y; (b) γ,x; g^E/RT,x. (○) Ref. [34] Ref. [35] (+) Ref. [26] (◊) Ref. [37]; (×) Ref. [29]; (∇) Ref. [32]; (▷) Ref. [31]; (◁) Ref. [37].

Figure 4.

Plot of iso-p VLE of the binary benzene(1) + hexane(2). (a) T,x,y; (b) g^E/RT,x. Ref. [33]: (○) p = 26 kPa; (□) p = 40 kPa; (+) p = 53 kPa; (◊) Ref. [26]; Ref. [28]: (×) p = 405.4 kPa; (△) p = 610.8 kPa; (▽) p = 810.8 kPa; (▷) p = 1013.5 kPa.

Figure 4 depicts the iso-p data set at pressures other than the atmospheric. Figure 5(b) indicates that three of the published series [28] (the ones that correspond to p = 610.8 kPa, 810.8 kPa, and 1013.5 kPa) show systematic errors, giving rise to g^E/RT < 0; thus γ_i < 1 as x₁ → 1. However, the representation of T vs x,y does not reflect this anomalous behavior, so that probably they will be caused by a systematic deficiency in the monitored temperature.

Figure 5.

Plot of iso-T VLE data of the binary benzene(1) + hexane(2). (a) p,x,y; (b) g^E/RT,x. (○) Ref. [40]; (□) Ref. [42] Ref. [41]: (+) T = 314 K; (◊) T = 323 K; (×) T = 333 K; (△) T = 303 K; (▽) Ref. [39].

Available iso-T data series are represented in Figure 5. Data are found for the whole composition interval only at 333 and 345 K. Several series were measured at 333 K [40, 41, 42], allowing their comparison. From the representation of p,x,y (Figure 5(a)) and related mixing quantities (Figure 5(b)), an acceptable coincidence is observed. Two of the data sets (Ref. [41]) show negative values of g^E/RT in the extreme points, at infinite dilution, which is a clear sign of inconsistency.

Conclusions drawn from the visual inspection are confirmed by the consistency analysis presented in Table 2. Data series n° 11 is discarded since consistency cannot be evaluated due to the lack of data in the zone corresponding to x₁ > 0.2. All consistency tests, except the direct van Ness, did not accept the series measured at the highest pressures (n° 13–16). This is due to the wrong relationship between the pure component saturation temperature (Appendix Table A1) and equilibrium temperature. However, it is interesting to recognize that when using other vapor pressure parameters this defect can be solved. The data series (n° 3, 4, 6–9, 19, 21–23) accepted by the proposed test [7] offer the highest consensus among all the battery of test about the quality of data, so these series could be used for subsequent tasks.

N°	Ref.	Type	Area		Fredenslund		Wisniak		Kojima		Direct-Van Ness		Proposed test
			d%	V	×100	V	Dw	V	M(Ii)	V	δln(γ1/γ2)	V	f-int.	V	f-dif.	V	FV
1	33	Iso-26 kPa	23.2	nv	1.2	nv	1.45	v	74	nv	0.123	v	−0.06	nv	−0.073	nv	nv
2	33	Iso-40 kPa	16.3	nv	1.2	nv	1.02	v	19	v	0.068	v	−0.02	nv	−0.013	nv	nv
3	33	Iso-53 kPa	17.5	nv	0.8	v	1.93	v	31	nv	0.063	v	0.04	v	0.001	v	v
4	26	Iso-97 kPa	19.2	nv	0.6	v	1.7	v	54	nv	0.02	v	0.05	v	0.005	v	v
5	31	Iso-101 kPa	36.5	nv	0.6	v	1.65	v	10	v	0.012	v	0.02	v	−0.047	nv	nv
6	26	Iso-101 kPa	3.6	nv	0.4	v	1.6	v	24	v	0.020	v	0.11	v	0.006	v	v
7	34	Iso-101 kPa	2.5	nv	0.3	v	1.63	v	24	v	0.004	v	0.04	v	0.007	v	v
8	35	Iso-101 kPa	1.5	v	0.1	v	1.7	v	27	v	0.003	v	0.01	v	0.014	v	v
9	37	iso-101 kPa	9	nv	0.4	v	1.75	v	44	nv	0.159	v	0.05	v	0.004	v	v
10	29	iso-101 kPa	12.3	nv	0.6	v	1.66	v	22	v	0.048	v	−0.01	nv	−0.013	nv	nv
11	32	Iso-101 kPa	nd	—	nd	—	nd	—	nd	—	nd	—	nd	—	nd	—	nv
12	37	Iso-101 kPa	9	nv	0.7	v	2.77	v	52	nv	0.043	v	0.02	v	−0.018	nv	nv
13	28	Iso-405 kPa	59	nv	1.2	nv	5.9	nv	1319	nv	0.2	nv	−0.21	nv	−0.145	nv	nv
14	28	Iso-610 kPa	58.9	nv	1.1	nv	7.51	nv	343	nv	0.043	v	−0.14	nv	−0.079	nv	nv
15	28	Iso-810 kPa	96.3	nv	2	nv	11.5	nv	291	nv	0.124	v	−0.16	nv	−0.08	nv	nv
16	28	Iso-1010 kPa	95.4	nv	3	nv	14.97	nv	6776	nv	0.364	nv	−0.07	nv	−0.192	nv	nv
17	41	Iso-303 K	16.8	nv	1.6	nv	0.61	v	670	nv	0.167	nv	0.03	v	−0.096	nv	nv
18	41	Iso-314 K	17.5	nv	2.7	nv	0.68	v	1149	nv	0.312	nv	0	nv	−0.122	nv	nv
19	41	Iso-323 K	7.9	nv	0.6	v	0.97	v	270	nv	0.048	v	0.13	v	0.006	v	v
20	40	Iso-333 K	31.1	nv	1.8	nv	2.42	v	41	nv	0.082	v	−0.33	nv	−0.048	nv	nv
21	42	Iso-333 K	1.7	v	0.1	v	1.31	v	12	v	0.002	v	0.5	v	0.017	v	v
22	41	Iso-333 K	1.2	v	0.4	v	1.22	v	88	nv	0.017	v	0.41	v	0.009	v	v
23	39	Iso-343 K	7.4	nv	0.5	v	1.53	v	19	v	0.027	v	0.29	v	0.001	v	v

Table 2.

Values obtained in the application of the thermodynamic consistency-test to VLE data of benzene + hexane at different conditions.

Limiting values: areas-test: d < 2; Fredenslund-test: δy×100 < 1; Wisniak-test: D_w < 3; Kojima-test: M(I_i) < 30; direct Van Ness-test: δ(ln γ₁/ln γ₂) < 0.16; proposed-test: f-int > 0, f-dif > 0. Header notation: V → verified; FV → fully verified. Test outcome: v → verified; nv → not verified; nd/− → not available.

A total of 17 references were found in literature [36, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58] of h^E measured in the range of [293–323] K and atmospheric pressure, except those data of Yi et al. [58], at 80 kPa. Figure 6(a) presents the data series considered in this study. It refers to a solution with endothermic effects (h^E ≈ 900 J mol⁻¹ at x = 0.5 and T = 298 K), decreasing as temperature increases; however, the high observed dispersion of experimental values is not sufficient evidence to ensure this effect. Inspection of the h^E does not recommend excluding any of these data series; thus they will be considered in the data treatment.

Figure 6.

Plot of h^E values for the binary benzene + hexane. (a) h^E vs x, (b) h^E (at x = 0.5) vs T. (○) Ref. [43], (□) Ref. [44], (+) Ref. [45], (◊) Ref. [46], (×) Ref. [47], (△) Ref. [48], (▽) Ref. [36], (▷) Ref. [49], (◁) Ref. [50], (□) Ref. [51], Ref. [52], Ref. [53], Ref. [54], Ref. [55], Ref. [56], Ref. [57], Ref. [58].

4. Results and discussion

4.1 Modeling: efficient front and models

The VLE data series previously verified has been modeled according to the procedure described in Section 2. The efficient fronts obtained for each system are shown in Figure 7.

Figure 7.

Efficient fronts for s(g^E/RT) = f(s(h^E)), obtained for acetone + ethanol. Arrows and labels indicate the chosen results on each of the fronts, using the models: (▴) M2; (×) M3; and (▾) M4.

Acetone(1) + ethanol(2): the efficient fronts of results for the binary acetone + ethanol were obtained from sub-models M2, M3, and M4. Sub-model M1 did not produce an acceptable representation of the thermodynamic behavior, so it was ignored. The obtained fronts shown in Figure 7 reveal that the three sub-models produced similar results when acceptable error in h^E is greater than 140 J mol⁻¹. This implies that, from this limit on, the problem becomes mono-objective, regardless of the sub-model used. For smaller errors of the h^E, s(h^E) < 140 J mol⁻¹, differences between sub-models M2, M3, and M4 to reproduce the VLE data are significant, reducing the maximum error by approximately δs(g^E/RT) ≈ 0.15 between them. The efficient front achieved with model M4 shows an almost constant s(g^E/RT). The other two sub-models reveal an exponential behavior as s(h^E) decreases. From the set of results obtained, three of them are chosen (see Figure 7) to carry out a more detailed analysis on their ability to describe simultaneously the VLE and h^E. Some particular comments on those three results are:

• The result labeled as “1” (by M4) describes well the two properties under study, as seen in Figure 8, but at the expense of using more parameters. Specifically, this model links the state variables at equilibrium (Figure 8(a)), describing the observed folding at the lowest temperatures. The estimate of activity coefficients (Figure 8(b)) is acceptable. Estimations of g^E/RT for this system are close to the set that displays the highest values. The estimate for the h^E is acceptable, although the model does not reproduce the data series extracted from Ref. [20].

• Result 2 (M3), which displays the highest error in the VLE estimate, wrongly estimates an azeotrope and a whole range immiscible region (Figure 8(a)). This poor description of the liquid phase is due to the high values of γ_i and the g^E/RT (Figure 8(b and c)). At first, the h^E is correctly represented by this parametrization, although the presence of the immiscibility truncates the validity region of the model (Figure 8(d)).

• Result 3 (M3) produces an unsatisfactory representation in the h^Es, which also presents an inversion of the thermal gradient of this property.

Figure 8.

Plot of VLE at 101 kPa and h^E estimates for binary acetone(1) + ethanol(2) at 101.32 kPa. Models drawn from Figure 7(a). Results (——) 1-(M4), 2-(M3), 3-(M3). (a) T vs x₁,y₁; (b)γ_i vs x₁; (c) g^E/RT vs x₁; (d) h^E vs x₁.

Final model selection should ensure that selected parametrization describes, at least qualitatively, the thermodynamic behavior closest to the real one. In this case, result 2 reproduces almost exactly the h^E, although it produces an incorrect g^E behavior, hence limiting its subsequent applicability. The choice between results 1 and 3 will depend on the influence of the h^E on the calculations in which the model is involved.

This analysis will serve as a basis to inspect the results of the remaining binary (benzene + hexane) on the results emitted by each of the models.

Benzene(1) + hexane(2): sub-model M4 was not used for the binary benzene + hexane because of the similarity with the efficient front produced by M3 model. Thus, sub-models M1, M2, and M3 were only applied, as shown in Figure 9. Sub-models M2 and M3 produce very similar results, with almost constant error, s(g^E/RT), lower than 5.5 × 10⁻³ J mol⁻¹ for all s(h^E). Consequently, those result fronts of sub-models M1 and M2 are evaluated. The efficient front of sub-model M1 produces a quasi-linear behavior with the variation of s(h^E). The difference between this front and that of sub-model M2 is δs(g^E/RT) ≈ 0.14, when ε = 0. The fronts for models M1 and M2 do not intersect, unlike the earlier case, showing a maximum value of s(h^E) at 170 J mol⁻¹ and 450 J mol⁻¹, respectively. The VLE diagrams produced by the four selected results are shown in Figure 10, along with one of the validated data series for this result. The best description of this system is achieved with result 4 (M1), which reproduces the behavior of T-x-y experimental data (Figure 10(a)) and the other quantities calculated (Figure 10(b and c)). Nevertheless, the description of h^E with this model is not good. Result 1 (M2) produces an azeotrope at x₁ < 0.2, which does not occur experimentally. This poor estimation occurs even though the greater number of parameters, increasing the model’s capacity to reproduce the h^E, as proven in Figure 11. Results 2 and 3 overestimate γ_i, hence g^E/RT (see Figure 10(b and c)).

Figure 9.

Efficient fronts for s(g^E/RT) = f(s(h^E)), obtained for benzene + hexane. Arrows and labels indicate the chosen results on each of the fronts, using the models: (●) M1; (▴) M2; and (×) M3.

Figure 10.

Plot of VLE at 101 kPa estimates for binary benzene(1) + hexane(2). Models drawn from Figure 9. Results (——) 1 (M2), 2 (M1), 3 (M1) and 4 (M1). (a) T vs x,y; (b)γ vs x; (c) g^E/RT vs x.

Figure 11.

Estimation of h^E values at 101.32 kPa using the different results indicated in the fronts of Figure 9 for benzene(1) + hexane(2). (a) h^E vs x₁ (T = 298.15 K), (b) h^E vs T (x₁ = 0.5). Results (——) 1 (M2); 2 (M1); 3 (M1); 4 (M1).

This discrepancy gives rise to the formation of minimum boiling point azeotropes, which are not in accordance with experimental data. Of all the results chosen, only result 1 (belonging to sub-model M2) shows a h^E that varies significantly with temperature, since sub-model M1 is independent of this variable. The use of either result 2 or 3 is discouraged since their estimations of h^E, especially at temperatures other than 298 K, are not correct, in addition to the described issues in representing the phase equilibria. This being said, result 3 represents the average behavior of this property in the working range.

From previous observations, result 4 (M1) is recommended for those cases where the h^E plays a secondary role in comparison to the reliable reproduction of the VLE diagram. Otherwise, result 1 (M2) is preferred, despite the qualitative misfit to experimental data. Table 3 presents an overview of the present section with a summary of the models to be used in the simulation task.

System	M1	M2	M3	M4
Acetone + ethanol	✖	✔	✔✔	✔✔
Benzene + hexane	✔✔	✔✔	✔	✖

Table 3.

List of sub-models applied to each system.

✔✔ → used for modeling and simulation;✔ → used for modeling; ✖ → not used.

4.2 Simulation results of a rectification process for each of the studied dissolutions

The models obtained previously are used in the design of a simulation operation comparing their capacity in terms of some operation variables such as composition and temperature profiles, as well as energy consumption. General conditions for the simulations are summarized in Table 4. In all cases, columns are fed with a 1 kmol/h at equimolar composition of the corresponding solution, at 298.15 K and 101.32 kPa. Simulations are performed using the RadFrac block of AspenPlus^© V8.8 (AspenOne^©, [59]).

Binary system	Reflux ratio	Distillate rate (kmol/h)	n° stages	Feed stage
Acetone + ethanol	6	0.5	22	16
Benzene + hexane	10	0.6	30	20

Table 4.

Operation data for the rectification columns to separate the binaries.

Acetone(1) + ethanol(2): the simulation of a rectification column to purify the acetone+ethanol binary system, using the result labeled as “2” in Figure 7, does not provide a coherent resolution because it estimates the presence of two immiscible liquid phases. The final values obtained with results 1 and 3 are detailed in Table 5, while the composition profiles are shown in Figure 12. In both cases the composition of the distillate is higher than 99% in acetone and at the same temperature in the stage.

Stage	Result 1 (M4)			Result 3 (M3)
	x1	y1	T/K	x1	y1	T/K
1	0.990	0.993	329.0	0.989	0.991	329.0
5	0.974	0.981	329.1	0.970	0.978	329.1
9	0.941	0.957	329.3	0.935	0.953	329.4
13	0.853	0.900	330.1	0.846	0.894	330.2
17	0.518	0.709	334.4	0.526	0.710	334.2
22	0.010	0.031	350.6	0.011	0.035	350.5
Qc/kJ h^−1	1.035E5			1.034E5
Qr/kJ h^−1	1.038E5			1.044E5

Table 5.

Quantities obtained in the simulation of a separation process for the binary acetone(1) + ethanol(2), using different values from the efficient front shown in Figure 7.

Figure 12.

Plot of composition and temperature profiles obtained in the simulation of a separation operation for acetone(1) + ethanol(2) system, using the different parametrizations proposed: (—) 1 (M4). L, liquid stream profile; V, vapor stream profile.

The exact purity is slightly higher with result 1 than with result 3, the difference being 0.2%. The calculation in the bottoms of the tower reveals differences between the two parametrizations, giving place to an effluent somewhat purer in the case of result 1. The difference between both models is 0.001 in molar fraction. These observations directly affect the calculation of the energy balance and, therefore, to the consumed energy. Thus, the consumption in the condenser is estimated similarly with both parametrizations, while that of the reboiler is significantly higher with result 3, due to the greater quantity of ethanol and the incorrect estimate of other quantities, such as the mixing enthalpies.

The composition profiles and temperature gradient of the two tested solutions present similar qualitative behaviors. Most of the column is used as an enriching region which requires a high number of stages in both cases.

Benzene(1) + hexane(2): values obtained in the simulation of the distillation operation of this binary dissolution are in Table 6, while composition and temperature profiles are plotted in Figure 13. For the first three results (1, 2, and 3), the presence of an azeotrope limits the distillate composition. Thus, result 1 produces an effluent with a benzene composition of 16.7% (v/v), while for results 2 and 3, the compositions are, respectively, 37.2 and 30.9%. The simulation carried out with the parametrization of result 4 produces a purer distillate of 15.6%, due to the absence of the azeotrope. However, the folding effect observed between the equilibrium curves in experimental data as well as the diagram estimated by result 4 complicates the separation beyond this point. This justifies that the composition profiles are similar to results 1 and 4 (see Figure 13).

Stage	Result 1 (M2)			Result 2 (M1)			Result 3 (M1)			Result 4 (M1)
	x1	y1	T/K	x1	y1	T/K	x1	y1	T/K	x1	y1	T/K
1	0.167	0.161	341.3	0.372	0.372	337.3	0.310	0.309	340.1	0.167	0.156	341.7
5	0.201	0.189	341.3	0.372	0.372	337.3	0.312	0.311	340.1	0.220	0.202	341.9
10	0.280	0.250	341.6	0.372	0.372	337.3	0.319	0.316	340.1	0.310	0.278	342.2
15	0.505	0.414	343.1	0.374	0.373	337.3	0.337	0.330	340.1	0.473	0.407	343.2
20	0.842	0.740	348.5	0.405	0.388	337.4	0.392	0.369	340.2	0.772	0.663	346.8
25	0.991	0.982	352.6	0.409	0.390	337.4	0.422	0.389	340.3	0.988	0.975	352.4
30	1.000	0.999	352.8	0.692	0.510	338.8	0.785	0.631	344.0	1.000	0.999	352.8
Qc/kJ h^−1	1.829E5			1.988E5			1.686E5			1.804E5
Qr/kJ h^−1	1.886E5			2.050E5			1.802E5			1.898E5

Table 6.

Quantities obtained in the simulation of a distillation process for the binary benzene(1) + hexane(2), using different values from the efficient front shown in Figure 9.

Figure 13.

Plot of composition and temperature profiles obtained in the simulation of a separation operation for the binary benzene(1) + hexane(2), using the different parametrizations proposed: (——) 1 (M2); (- - - -) 2 (M1); (······) 4 (M1). L, liquid stream profile; V, vapor stream profile.

The residual streams obtained by results 1 and 4 contain benzene with a purity higher than 99.9%, while the other two results do not produce the separation of the binary dissolution with high purities (no more than 79%), due to the presence of the azeotrope at intermediate composition. The differences noted in compositions in the head and bottom of the column cause that the corresponding temperatures are also different. Between results 1 and 4, there is a difference in temperature of almost 0.5 K, while results 2 and 3 estimate lower temperatures by more than 1 K with respect to the others. In this case, the modeling errors in another of properties, such as excess enthalpies, as well as the differences in the output composition, also affects to produce noteworthy differences in the energetic consumptions of both condenser and boiler. In this sense, similar results are obtained for results 1 and 4, although, a priori, it is not possible to indicate which of these is closer to the true behavior of the column.

5. Conclusions

The aim of this work is to present the practice carried out on a set of operations constituting a procedure involving the following sequences: experimentation ↔ verification ↔ modeling ↔ simulation. A description of the last three operations is proposed since the necessary experimental information (iso-p and iso-T VLE data and h^E) was extracted from the available publications. The proposed methodology is applied to two significant binary systems in the dissolution thermodynamic field, such as acetone-ethanol and benzene-hexane. The data checking is carried out with different classical methods, which is recommended in the recent literature [6]. In addition, a rigorous method recently designed by the authors [7] was also used in order to guarantee the quality of the experimental data used. This method has the advantage of offering some information about the origin of deficiencies in the experimental data.

A polynomial expression was used in the modeling step (see Eqs. (1) and (3)), used in the correlation of thermodynamic properties, from which four sub-models (M1–M4) are established in Section 2.1, depending on the availability of data for each system to avoid overfitting. Modeling in all cases was performed on the Gibbs dimensionless function, g^E/RT (VLE), and on h^E data, two-objective optimization approach addressed as a sequence of mono-objective subproblems according to the ε-constraint algorithm. A set of models are selected from the attained efficient fronts to provide a rationale on the trade-off decision task that is supposed to yield the final model for later uses, such as design/simulation. For the two studied cases, these fronts reveal a quasi-linear trend (with a negative slope) when the number of parameters used in the model is small. Efficient fronts’ slope decreases as the number of parameters increases, until the error of VLE delinks to that of the h^E for highly flexible models.

The rectification of the selected systems was simulated using the RadFrac block of AspenPlus^© with the selected thermodynamic parametrizations. For the binary acetone + ethanol, some models showed an immiscible region, not reflected by the experimental information, giving rise to incoherent simulations. Two models gave rise to inexistent azeotropes in the benzene + hexane dissolution, with incorrect results in the simulation. On the other hand, errors in the excess enthalpy estimates did not influence on the procedure simulation, but it should be noted that it has an important role in the modeling of binary systems. In summary, the influence of a certain set of parameters on the simulation varies depending on each particular system. Besides, it is observed that, as the number of parameters grows in a model, the optimization problem mutes toward a mono-objective one since the considered criteria are invariant in one another. In these cases, taking the set of parameters that present the lowest error on h^E is the best option. However, increasing the number of parameters might lead to overfitting if not enough attention is paid to the model extrapolation capabilities.

Acknowledgments

This work was supported by the Ministerio de Economía y Competitividad (MINECO) of the Spanish Government, Grant CTQ2015-68428-P. One of us (AS) is also grateful to the ACIISI (from Canaries Government, No. 2015010110) for the support received. This work is a result of the Project “AIProcMat@N2020—Advanced Industrial Processes and Materials for a Sustainable Northern Region of Portugal 2020,” with the reference NORTE-01-0145-FEDER-000006, supported by Norte Portugal Regional Operational Programme (NORTE 2020), under the Portugal 2020 Partnership Agreement, through the European Regional Development Fund (ERDF); Associate Laboratory LSRE-LCM—UID/EQU/50020/2019—funded by national funds through FCT/MCTES (PIDDAC).

A.

Vapor pressures of pure compounds are calculated using Antoine equation:

$$log{p}_{\mathrm{i}}^{\mathrm{o}}/\mathit{kPa}=A-B/\left[T/K-C\right]$$ E10

Parameters A, B, and C, from literature [60, 61, 62], are shown in Table A1.

	A	B	C	Ref.
Acetone	6.42448	1312.25	32.45	[60]
Ethanol	7.24677	1598.67	46.42	[61]
Benzene	6.03053	1211.03	52.36	[62]
Hexane	5.87891	1089.49	60.55	[62]

Table A1.

Parameters of Antoine equation for pure compounds in this work.