Stefan
Eckelsbach
and
Jadran
Vrabec
*
Thermodynamics and Energy Technology, University of Paderborn, 33098 Paderborn, Germany. E-mail: jadran.vrabec@upb.de; Fax: +49 5251 60 3522; Tel: +49 5251 60 2421
First published on 21st September 2015
Vapor–liquid equilibria (VLE) of the pure substances acetone, oxygen and nitrogen as well as their binary mixtures are studied by molecular dynamics (MD) simulation with a direct approach. Thereby, particular attention is paid to the interface behavior on the molecular level, yielding total and partial density profiles as well as surface tension data. The classical approach by van der Waals is used to analyze the total density profiles. It is found that an extended function is needed to describe those profiles for the mixtures containing acetone, due to the strong adsorption of the volatile component at the vapor side of the interface. Based on these representations the interface thickness is studied. The surface tension results are compared to experimental data, correlations thereof and results from other molecular approaches. Due to the scarcity of experiments, the parachor method is employed to obtain predictive surface tension data for the mixtures. Following the same approach, the present surface tension results are correlated for the mixtures containing acetone.
In contrast to phenomenological approaches, molecular modeling and simulation is based on a sound physical foundation and is therefore well suited for the prediction of such properties and processes under extreme ambient conditions. In this work a model fuel (consisting of acetone, oxygen and nitrogen) was investigated, using mainly molecular modeling and simulation. MD simulations of VLE interfaces were carried out, focusing on high temperatures in the region that is not too far from the critical line of the binary mixtures containing acetone. Thereby, interface behavior was studied on the molecular level with respect to important thermodynamic properties such as the surface tension. The results are provided as reference values for larger scale models, i.e. computational fluid dynamics, within the collaborative research center SFB-TRR 75 (“Droplet dynamics under extreme ambient conditions”), which investigates this model fuel, particularly the injection of acetone droplets into a gaseous oxygen + nitrogen environment under high pressure.2
The present simulation volumes were created as rectangular cuboids, elongated in z direction, with a liquid phase surrounded by the coexisting vapor phase under periodic boundary conditions. The saturated densities of the bulk phases that are needed for this initial setup were taken from earlier simulations, likewise the equilibrium composition of the phases in case of mixtures.6 Due to the inhomogeneity of the systems, a long range correction has to be used, which takes the varying density along the z axis into account. The approach of Janeček8 solves this problem and different implementations based on his method have recently been suggested.9–12 Here, the implementation by Werth et al.13 was used for molecular interactions exceeding the cutoff radius of 2.2 nm. This method was extended for multi-site Lennard-Jones models,13 as they were used in this work, by the conjunction with a centre-of-mass cutoff scheme based on angle averaging.14 It was implemented in the simulation code to correct the potential energy, virial and force. This approach is based on the continuous calculation of the density profiles throughout the simulation, which were also evaluated in this work. However, another formulation of the long range correction, based on Janeček's method, was recently suggested by MacDowell and Blas9 as well as Martínez-Ruiz et al.15 avoiding the continuous evaluation of the density profile. Based on previous work of our group on the regarded systems,6,16 including the binary mixtures, the knowledge of the VLE data was used to straightforwardly conduct simulations in the canonical ensemble. After sufficient equilibration, the molecular systems were sampled with a time step length of 1 fs over 2 ns containing 32000 and for state points closer to the critical region 48000 molecules. Further details of the simulation setup are given as ESI.†
(1) |
The density profile across the VLE interface of pure acetone is exemplarily plotted in Fig. 1 for a temperature of 450 K. Note that the critical temperature is about Tc = 508.1 K.23 The simulation data can well be described with function (1) and show only small sinusoidal deviations on the vapor side as seen in prior work for a simple Lennard-Jones model fluid.21 The relative deviations on this side are larger than on the liquid side of the interface mainly because of statistical reasons. First, there are significantly less molecules in the vapor phase, such that the statistical uncertainties are higher, second, the relative deviations are larger due to the lower vapor density. Nonetheless, the relative deviations are small in the entire considered spatial range, i.e. below 1%, and the data in the interfacial region are well represented showing deviations on roughly the same level.
Fig. 1 Density profile and deviations from the fitting function in case of acetone at 450 K through the VLE interface. ○ simulation results; fit of eqn (1) to the simulation results. The deviations are given as δρ = (ρ − ρfit)/ρfit. |
For the other simulations of acetone as well as pure oxygen and nitrogen, the qualitative behavior of the density profiles was found to be very similar over the whole temperature range and can also be accurately represented by eqn (1). The quantitative results appear as expected. The interfacial thickness D widens and the difference between the saturated bulk densities ρl and ρv declines for higher temperatures closer to the critical point due to the monotonic decrease of the liquid and the monotonic increase of the vapor density. The numerical simulation results are provided as ESI.†
For the binary mixture oxygen + nitrogen at 120 K and an equimolar liquid bulk composition, the total and partial density profiles are depicted in Fig. 2. In the interface region, the partial densities of oxygen and nitrogen show a similar slope in their transition from liquid to vapor density, but with a spatial shift of about 0.28 nm. Thus nitrogen, being the lower boiling component, shows an adsorption layer at the vapor side of the interface, which was also found for different binary mixtures in other work.24–27 This behavior can be explained by the properties of the molecular force field models.24,26,28 The attraction to the liquid phase is larger for both components, due to its higher density. For Lennard-Jones models with similar size parameter, which is the case here (σO = 3.1062 Å and σN = 3.3211 Å), this leads to an easier transition into the liquid phase for the substance with stronger dispersive attraction and an accumulation of the substance with weaker dispersive attraction at the interface. The according parameters of the molecular models are εO/kB = 43.183 K and εN/kB = 34.897 K. The total density profile data were fitted with eqn (1) and show sinusoidal relative deviations as seen for the pure substances. Here, the deviations in the interfacial range are more pronounced, but nonetheless they are only slightly larger than the scatter in the vapor phase. Thus, it was concluded that the total density profile of this mixture can be represented with a sufficient accuracy by function (1).
Fig. 2 Total and partial density profiles as well as deviations from the fitting function through the VLE interface of the mixture oxygen + nitrogen at 120 K and an equimolar liquid bulk composition. ○ total density and deviations from the fitting function; partial density of oxygen; – – partial density of nitrogen; — fit of eqn (1) to the simulation results. The deviations are given as δρ = (ρ − ρfit)/ρfit. |
For a given temperature, the interfacial thickness D of oxygen + nitrogen was found to be between the interface widths of the pure constituent substances. However, regarding the temperature reduced by its critical value of the considered fluid, this binary mixture forms a wider interface than its constituent pure components, which is shown in Fig. 3. This is due to the mixture's more pronounced increase of the interfacial width in relation to its critical temperature, compared to the pure oxygen data, and is caused by the adsorption of nitrogen at the interface (cf.Fig. 2). It can also be seen that the slopes of the simulation results are very similar for oxygen, nitrogen and their mixture. This likewise applies to the acetone data, which are plotted in the ESI,† although their progression towards the critical point is slightly steeper.
Fig. 3 Width of the interface over the reduced temperature TTc−1 of oxygen and nitrogen as well as their binary mixture with equimolar liquid bulk composition. ◑ oxygen; ◐ nitrogen; ○ oxygen + nitrogen with equimolar liquid bulk composition. The lines represent the empirical correlation21D = a(1 − T/Tc)b + c(1 − T/Tc)−0.63 fitted to the data. For the determination of the critical temperatures of the pure substances, the parameters of the corresponding equations of state were applied.29,30 Based on these, the GERG-2004 equation of state31 was used to calculate the critical temperature of the mixture. The statistical uncertainties of the simulation data were considered as the standard error on a 95% confidence interval and are only visible if they exceed the symbol size. |
The two mixtures containing acetone show a qualitatively different density and partial density behavior than the mixture oxygen + nitrogen. The attraction between acetone molecules is much more pronounced, not only due to its high Lennard-Jones energy parameters, but also due to its strong dipole moment μ = 3.4 D. Therefore, the adsorption of the volatile component at the interface is more distinct in this case. For simulations not too close to the critical line, this leads to a maximum of the partial density in the interface, which is a known phenomenon.24,26–28,32 This is exemplarily shown in Fig. 4 for oxygen + acetone at 400 K with an oxygen mole fraction in the liquid bulk phase xO2 = 0.05 mol mol−1. In the upper subplot, the density and partial density profiles across the VLE interface are depicted over the whole density range of the liquid and vapor phase. The liquid phase mainly consists of acetone and to a much lesser extent of oxygen, which is represented by the difference in the partial densities. On the contrary, the vapor phase is composed mostly of oxygen. The lower subplot shows the same simulation data scaled up for the interesting region on the level of the vapor density, which clearly exhibits the enrichment of oxygen in the interface, exceeding its partial liquid density. Of course, this accumulation also influences the total density profile.
Due to this behavior, the density profile across the interface is not represented very well by function (1), especially on the vapor side. Hence, the extended equation21
(2) |
A comparison of the two functions (1) and (2) is exemplarily shown in Fig. 5 for this mixture (oxygen + acetone at 400 K and xO2 = 0.05 mol mol−1). In the upper panel, the simulation values for the total density and both fitting functions are plotted. It can be seen that eqn (2) provides a significantly better description of the MD data. This is confirmed by the lower subplot, where the deviations of the simulation results from the functions are shown. The largest deviations of 2.6% occur for eqn (1) in the transition region from the interface towards the vapor phase, which are much higher than those in the remaining spatial range, indicating its inability to adequately describe the profile data. On the other hand, the deviations from eqn (2) are on a consistent level over the whole spatial range and do not exceed 0.6%. For the residual sum of squares in the considered region, this means an improvement from 1.9 × 10−2 to 7.1 × 10−4.
Fig. 5 Total density profile and deviations from the fitting functions for the mixture oxygen + acetone at 400 K and xO2 = 0.05 mol mol−1 on the vapor side of the VLE interface. ○ simulation results; , ◑ fit of eqn (1) to the simulation results and their deviation; , ● fit of eqn (2) to the simulation results and their deviation. The deviations are given as δρ = (ρ − ρfit)/ρfit. |
This behavior can be seen over a wide range of composition and for both simulated isotherms of the mixture oxygen + acetone. It can be identified by the difference between the interfacial widths Dl and Dv, which are plotted in Fig. 6. For pure acetone it was discussed above that eqn (1) suffices in describing the density profile (cf.Fig. 1) so that Dl and Dv as well as zl and zv are equal for xO2 = 0 mol mol−1 if eqn (2) is used to represent the simulation data. In the range xO2 = 0.04–0.20 mol mol−1 with pressures reaching from 4.6 to about 20 MPa,6 the simulations show a strong accumulation of oxygen on the vapor side of the interface, thus eqn (2) is needed to describe the total density profiles accurately. Accordingly, the widths show significant differences here. For compositions closer to the critical line, the accumulation of oxygen in the interface vanishes and Dl and Dv converge. This can be understood by considering the convergence of the saturated bulk densities approaching the critical region and finally the disappearance of the phase differences at the critical line.
The widths of the interface of the mixture nitrogen + acetone are shown in Fig. 7. Compared to the results for oxygen + acetone, the slope of the widths is steeper and the lenticular shape is more distinct. Qualitatively, the appearance of the data is similar, also showing increasing differences starting from the pure substance and the agreement of Dl and Dv within their error bars close to criticality.
The convergence of Dl and Dv indicates that eqn (1) can be used to accurately describe the total density profile for high liquid mole fraction values of the more volatile component. In Fig. 8 a state point in the extended critical region is shown for the mixture nitrogen + acetone with a liquid bulk mole fraction of 0.32 mol mol−1 and a corresponding pressure of 43 MPa.6 The proximity to the critical line is obvious considering the small difference between the saturated bulk densities, which is about 3 mol l−1, i.e. only 22% of the liquid density. The plotted total density data show that the accumulation of nitrogen at the vapor side of the interface occurs only slightly and the total density profile is well represented by function (1).
Fig. 8 Total and partial density profiles as well as deviations from the fitting function through the VLE interface of the mixture nitrogen + acetone at 400 K and xN2 = 0.32 mol mol−1. ○ total density and deviations from the fitting function; partial density of nitrogen; – – partial density of acetone; fit of eqn (1) to the simulation results. The deviations are given as δρ = (ρ − ρfit)/ρfit. |
The results from the evaluation of the total density profiles for the mixtures oxygen + nitrogen, oxygen + acetone and nitrogen + acetone were obtained by fitting the discussed functions to the block averaged profiles from present MD simulations with the Levenberg–Marquardt algorithm,33 followed by statistically analyzing the data. For fitting the data with eqn (1), all parameters were adjusted at once. However, due to the increase in the degrees of freedom, the fitting procedure with eqn (2) was found to be less stable. Here, the saturated bulk densities ρl and ρv were determined prior to fitting the other parameters. The numerical results are given in Tables 1–3, respectively.
T (K) | ρ l (mol l−1) | ρ v (mol l−1) | D (nm) | γ (mN m−1) |
---|---|---|---|---|
70 | 33.801 (4) | 0.020 (2) | 0.592 (6) | 14.7 (3) |
80 | 32.239 (5) | 0.123 (2) | 0.697 (6) | 12.2 (3) |
90 | 30.601 (5) | 0.346 (5) | 0.823 (6) | 9.7 (2) |
100 | 28.854 (6) | 0.721 (9) | 0.988 (8) | 7.3 (3) |
110 | 26.891 (9) | 1.35 (1) | 1.22 (1) | 5.1 (3) |
120 | 24.63 (2) | 2.39 (2) | 1.61 (2) | 3.2 (2) |
130 | 21.73 (3) | 4.19 (3) | 2.37 (4) | 1.3 (2) |
T (K) | x O2 (mol mol−1) | ρ l (mol l−1) | ρ v (mol l−1) | D l (nm) | D v (nm) | γ (mN m−1) |
---|---|---|---|---|---|---|
400 | 0.00 | 11.211 (5) | 0.284 (4) | 1.34 (1) | 1.35 (5) | 11.0 (4) |
400 | 0.05 | 11.682 (5) | 1.659 (5) | 1.45 (2) | 1.22 (3) | 9.8 (3) |
400 | 0.10 | 12.139 (4) | 2.989 (7) | 1.62 (2) | 1.1 (2) | 8.8 (3) |
400 | 0.15 | 12.605 (5) | 4.443 (5) | 1.70 (2) | 1.41 (6) | 7.2 (3) |
400 | 0.20 | 13.026 (4) | 5.717 (6) | 1.80 (3) | 1.53 (7) | 6.0 (3) |
400 | 0.24 | 13.378 (7) | 6.821 (6) | 1.88 (5) | 1.78 (6) | 5.1 (4) |
400 | 0.26 | 13.519 (7) | 7.218 (9) | 1.98 (5) | 1.84 (6) | 4.4 (4) |
450 | 0.00 | 9.656 (9) | 0.839 (9) | 1.90 (9) | 1.94 (7) | 5.0 (3) |
450 | 0.04 | 9.88 (1) | 1.613 (6) | 2.19 (4) | 1.78 (9) | 4.4 (2) |
450 | 0.05 | 9.952 (8) | 1.785 (9) | 2.24 (3) | 1.6 (2) | 4.3 (3) |
450 | 0.10 | 10.21 (1) | 2.716 (8) | 2.4 (1) | 1.9 (1) | 3.7 (2) |
450 | 0.15 | 10.49 (1) | 3.602 (9) | 2.66 (5) | 2.1 (3) | 2.9 (3) |
450 | 0.20 | 10.69 (3) | 4.300 (6) | 2.79 (5) | 2.4 (2) | 2.4 (3) |
450 | 0.24 | 10.87 (5) | 4.96 (2) | 3.5 (2) | 3.5 (3) | 2.3 (2) |
450 | 0.26 | 10.74 (8) | 4.91 (3) | 3.6 (2) | 3.4 (5) | 2.1 (2) |
T (K) | x N2 (mol mol−1) | ρ l (mol l−1) | ρ v (mol l−1) | D l (nm) | D v (nm) | γ (mN m−1) |
---|---|---|---|---|---|---|
400 | 0.00 | 11.211 (5) | 0.284 (4) | 1.34 (1) | 1.35 (5) | 11.0 (4) |
400 | 0.05 | 11.518 (2) | 1.983 (4) | 1.58 (1) | 1.0 (2) | 9.1 (3) |
400 | 0.10 | 11.857 (4) | 3.801 (4) | 1.70 (2) | 0.8 (3) | 7.3 (4) |
400 | 0.16 | 12.231 (6) | 5.72 (1) | 1.95 (3) | 1.6 (1) | 5.2 (4) |
400 | 0.24 | 12.824 (7) | 8.309 (6) | 2.36 (8) | 2.3 (1) | 3.1 (4) |
400 | 0.32 | 13.342 (8) | 10.37 (1) | 3.2 (1) | 2.9 (2) | 2.1 (4) |
450 | 0.00 | 9.656 (9) | 0.839 (9) | 1.90 (9) | 1.94 (7) | 5.0 (3) |
450 | 0.04 | 9.783 (9) | 1.743 (9) | 2.32 (4) | 1.7 (2) | 4.1 (2) |
450 | 0.05 | 9.80 (1) | 1.92 (1) | 2.43 (4) | 1.7 (2) | 3.9 (3) |
450 | 0.10 | 10.01 (1) | 3.151 (6) | 2.79 (4) | 1.6 (3) | 3.0 (3) |
450 | 0.15 | 10.18 (1) | 4.296 (9) | 2.96 (6) | 2.3 (1) | 2.0 (3) |
450 | 0.22 | 10.37 (1) | 5.759 (6) | 3.7 (1) | 3.0 (2) | 1.1 (3) |
450 | 0.26 | 10.42 (4) | 6.60 (3) | 4.8 (2) | 4.6 (3) | 0.7 (3) |
(3) |
(4) |
The surface tension of acetone is given in Fig. 9. Even for this common pure substance, the lack of experimental data is evident. Although there is quite a number of data points in the lower temperature region,37–42 the highest measurements reach up to only 351 K due to the flammability of acetone. The results from density functional theory27,43 of Klink and Gross44 slightly overestimate these laboratory data and a correlation of experimental data over the whole temperature range. For low temperatures, the results of the present work show a similar behavior as those of Klink and Gross, yielding deviations to experimental data of approximately 10%. This coincides with the results by Werth et al.,45 who have shown in their work that an overestimation of the surface tension on this magnitude is typical for molecular simulation data of many pure substances. At higher temperatures, the results are in better agreement with the correlation by Mulero et al.,46 which lies within the statistical uncertainties of the present MD simulation data.
Fig. 9 Surface tension over temperature for acetone. ○ molecular dynamics, this work; ● density functional theory, Klink and Gross;44 correlation of experimental data, Mulero et al.;46 + experimental data.37–42 The statistical uncertainties of the simulation data were considered as the standard error on a 95% confidence interval and are only visible if they exceed the symbol size. |
Experimental data for the considered mixtures are scarce, especially there are no experimental surface tension data at all. For the mixture oxygen + nitrogen, the surface tension can be determined by correlations46 and binary parameters for the underlying mixing rule.31 However, for the considered mixtures with acetone, there are no correlations available. A simple ansatz to estimate the surface tension is given by the parachor method,47 which can also be used to predict γ for multicomponent mixtures48
(5) |
In Fig. 10 the simulated surface tension of oxygen and nitrogen as well as their mixture with an equimolar liquid bulk composition from this work is compared to experimental reference data and to results by Neyt et al.,51 who used direct Monte Carlo simulations. For the pure substances, the data from the simulations show deviations to experimental data and correlations of about 12%, which is similar to the outcome for pure acetone. The present results for oxygen match very well with those by Neyt et al. In fact, the molecular model for oxygen used in their work was adapted from the one used here5 by substituting the point quadrupole of the original model with three electrostatic charges.52 For nitrogen, different force field models were used, yielding some slight differences between the results, however, they still agree within their statistical uncertainties. For the equimolar mixture, the parachor method was used to obtain additional data for comparison. Despite the very good input values for this mixture due to accurate EOS29–31 and γ correlations,46 this approach underestimates the surface tension over the whole temperature range compared to the correlation by Mulero et al.46 by an almost constant value of 0.6 mN m−1. The surface tension results from present MD simulations match with the correlation46 within their statistical uncertainties.
Fig. 10 Surface tension over temperature for oxygen and nitrogen and their mixture with equimolar liquid bulk phase. ○ molecular dynamics, this work; ● molecular dynamics, Neyt et al.;51 parachor method; correlation of experimental data, Mulero et al.;46 + experimental data.53–55 The statistical uncertainties of the simulation data were considered as the standard error on a 95% confidence interval and are only visible if they exceed the symbol size. |
In Fig. 11 the surface tension of the mixture nitrogen + acetone is plotted over the composition of the saturated liquid bulk phase for the isotherms 400 and 450 K, respectively. The values from present MD simulations as well as their slope nearly coincide with the results from density functional theory by Klink and Gross.44 For small xN2, the surface tension from MD simulation is slightly lower, yielding a good agreement with the correlation of experimental data for pure acetone.46 The prediction with the parachor method starts at exactly this correlation value, because its parameters were adjusted to it. However, with increasing xN2 this method underestimates the surface tension compared to the other results.
Fig. 11 Surface tension over the nitrogen mole fraction of the liquid bulk phase xN2 at 400 and 450 K for the mixture nitrogen + acetone. ○ molecular dynamics, this work; ● density functional theory, Klink and Gross;44 × correlation of experimental data, Mulero et al.;46 parachor method; fit of eqn (6) to the present molecular simulation data. The statistical uncertainties of the simulation data are given as the standard error on a 95% confidence interval. |
For the mixture oxygen + acetone, the surface tension over the liquid bulk mole fraction xO2 is shown in Fig. 12 for the isotherms 400 and 450 K. The results of the present MD simulations are higher than those of the mixture nitrogen + acetone, but resemble their shape. The surface tension values predicted by the parachor method are again significantly smaller than those from MD simulation.
Fig. 12 Surface tension over the oxygen mole fraction of the liquid bulk phase xO2 at 400 and 450 K for the mixture oxygen + acetone. ○ molecular dynamics, this work; × correlation of experimental data, Mulero et al.;46 parachor method; fit of eqn (6) to the present molecular simulation data. The statistical uncertainties of the simulation data are given as the standard error on a 95% confidence interval. |
Even though the parachor method is not very reliable for the prediction of the surface tension for the two mixtures containing acetone, it still provides a model for the correlation of the simulation results. For the mixtures oxygen + acetone and nitrogen + acetone, the present MD surface tension data were used to parametrize the parachor-like function
(6) |
Component | T (K) | K 1 ((mN m−1)1/n l mol−1) | K 2 ((mN m−1)1/n l mol−1) | n | |
---|---|---|---|---|---|
1 | 2 | ||||
O2 | C3H6O | 400 | 0.0015103 | 0.19868 | 3.6662 |
O2 | C3H6O | 450 | 0.18508 | 0.33430 | 1.7121 |
N2 | C3H6O | 400 | 0.15825 | 0.32052 | 2.1214 |
N2 | C3H6O | 450 | 0.084631 | 0.23123 | 2.8564 |
The numerical surface tension results from this work are listed in Tables 1–3 for the three different binary mixtures (oxygen + nitrogen, oxygen + acetone and nitrogen + acetone). For the calculation of the statistical uncertainties in this work, the data of the productive runs with a simulation time of 2 ns were divided into 32 blocks and the errors were estimated from the standard deviation of the mean.
For the pure substances it was found that the present MD simulations overestimate the surface tension compared to experimental data and correlations by approximately 10 to 12%. However, the results are in line with the data from other simulations and deviations on this order were found to be typical for molecular simulations of many pure substances in earlier work. Due to the scarcity of experiments for the three considered binary mixtures (oxygen + nitrogen, oxygen + acetone and nitrogen + acetone), predictions from the parachor method were considered for comparison. Where available, the correlations of experimental data agree well with the present MD simulation results, however the predictions from the parachor method underestimate the other data consistently. For the mixture nitrogen + acetone, results from density functional theory were available and show a very good agreement with the surface tension data from present MD simulations. For the two mixtures containing acetone, parachor-like correlations were parametrized to describe the surface tension results from this work.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5cp03415a |
This journal is © the Owner Societies 2015 |