- Email: [email protected]
Non-Lorentz–Berthelot Lennard-Jones mixtures: A systematic study
Fluid Phase Equilibria 277 (2009) 42–48
Contents lists available at ScienceDirect
Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid
Non-Lorentz–Berthelot Lennard-Jones mixtures: A systematic study Michael Rouha a , Ivo Nezbeda a,b,∗ a b
Faculty of Science, J. E. Purkinje University, 400 96 Ústí nad Labem, Czech Republic E. Hála Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, Academy of Sciences, 165 02 Prague 6 – Suchdol, Czech Republic
a r t i c l e
i n f o
Article history: Received 21 October 2008 Received in revised form 5 November 2008 Accepted 6 November 2008 Available online 18 November 2008 Keywords: Binary mixtures Lennard-Jones fluids Non-Lorentz–Berthelot rules Excess mixing properties Partial molar properties Tikhonov regularization
a b s t r a c t Binary mixtures of two identical Lennard-Jones fluids with non-Lorentz–Berthelot combining rules have been simulated over the entire concentration range at three state points in order to examine the effect of cross interactions on the mixing properties, excess volumes and enthalpies, and partial molar volumes, and on the structure. Various combinations of deviations, in both the energy and size cross parameters, from the Lorentz–Berthelot rule, have been considered and the results analyzed using a recently proposed method based on the Tikhonov regularization. It is found that the most important is the size cross parameter and that by varying the combining rules a variety of qualitatively different behavior can be produced including, e.g., a minimum in the partial molar volumes observed otherwise only in complex mixtures of associating fluids; occurrence of this minimum seems to be associated with changes in the second coordination shell as witnessed by the pair correlation function. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Knowledge of the thermodynamic behavior of fluid mixtures provides not only the fundamental prerequisite for design and operation of industrial processes but also fundamental information on the intermolecular interactions. According to the contemporary level of understanding and modeling of molecules and their interactions, molecules are viewed as bodies with two types of interaction sites: neutral ones, typically interacting as Lennard-Jones (LJ) or exp-6 sites, and the Coulombic sites [1–5]. Parameters of these sites are then determined by their fitting to a set of selected experimental data of pure fluids. During this parametrization, interactions between the unlike sites are estimated by means of certain combining rules. Whereas there is no problem with the interaction between the Coulombic sites, the interaction between, and hence the combining rules for the unlike neutral sites, is rather problematic. In general, the geometric mean rule (the Bethelot rule) is used for the energy parameter (the depth of the potential well), and the arithmetic rule (the Lorentz rule) for the length-scaling (excluded volume) parameter. Although there is only a very vague theoretical justification of the Lorentz–Berthelot (LB) rule [6–8], due to the adjustment of the results to experimental data its inadequacy does not cause principal problems in developing pure fluid interaction models. However,
∗ Corresponding author at: Faculty of Science, J. E. Purkinje University, 400 96 Ústí nad Labem, Czech Republic. Tel.: +420 220390 296; fax: +420 220920 661. E-mail address: [email protected] (I. Nezbeda). 0378-3812/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2008.11.007
the situation is quite different when dealing with mixtures. In this case the developed pure fluid interaction models are used to predict the properties of mixtures and it is evident that the entire concept of combining rules is, in general, ill-footed: even the used pure fluid potentials cannot adequately describe the interaction between the like species in the mixture because they do not account for the presence of the other species. It has been thus well known for decades that the mixture properties are very sensitive to the cross interactions but this is only a consequence of the entire concept. Since a re-parametrization of the pure fluid potentials with respect to more realistic cross interactions seems quite an unrealistic goal, an improvement of the description of mixtures rests on an improvement of the combining rules. One such a straightforward, and physically fully justified possibility provides inclusion of polarizability which, however, would cause nearly unsurmountable problems for theory and makes also molecular simulations very complex and hence very time consuming. Simple physical considerations show that the main problem lies in the cross LJ interactions and no wonder therefore that a number of papers on the LJ mixtures and combining rules have thus been published (for references see, e.g., the recent papers [9–11] and references therein). Nonetheless, the fundamental question whether deviations from the LB rule may produce qualitatively different mixture properties (and what type) does not seem to have been addressed at all. Typically, either certain modifications of the LB rule have been considered and examined or even different rules have been tested for a certain class of mixtures focusing only on mixing or equilibrium properties and seeking always only a better numerical agreement [9–12]. To our best knowledge we are thus not aware of any fundamental
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
and systematic study on the effect of the combining rules on the qualitative changes in the behavior of mixtures. To study effects of deviations of the cross interactions from those given by the LB rule, Duh et al. [13] considered some time ago two identical LJ fluids at a supercritical temperature. However, they considered in simulations only two compositions and only one set of the cross parameters. Adopting the same approach, i.e., to consider the mixture of two identical LJ fluids, we carried out recently a feasibility study considering a large number of non-Lorentz–Berthelot (nLB) cross parameters but only at one state point [14]. Furthermore, we focused on and reported results for excess properties only. Nonetheless, the shape of the excess volume for certain combinations of the nLB cross parameters indicated that the mixture may also exhibit qualitative changes in its behavior, e.g., exhibiting the minimum in the partial molar volume instead of a monotonous dependence on the concentration. We have therefore continued in investigations of mixtures along the same line focusing on the partial molar quantities and structure. Since the evaluation of the partial molar quantities requires a large number of accurate data on excess functions in dependence on concentrations (and hence lengthy time consuming simulations), in the present paper we consider only selected combinations of the cross parameters looking particularly for qualitative changes in the thermodynamic behavior. These changes must result from and reflect changes in the structure (in the case of low concentrations changes in the solvation shell) and we therefore compute also the structure factors and correlation functions which may reveal them. We also consider three different state points to find the impact of thermodynamic conditions (if there is any) on the behavior of the mixture. Finally, from the numerical point of view, there is a novelty in elaboration of the results: to analyze the obtained simulation data and evaluate derivatives (i.e., partial molar quantities) from noisy experimental data we use (and test) a recently proposed method based on the Tikhonov regularization [15] instead of the common Redlich–Kister parametrization [16].
Standard NPT Monte Carlo (MC) simulations [17,18] on binary mixtures of LJ fluids, uLJ (r) = 4εij
ij
r
−
6 ij
r
,
(1)
were performed. Both components of the mixture were identical, 11 = 22 ≡ ;
ε11 = ε22 ≡ ε,
(2)
but the cross interactions were varied in order to find out how they affect the thermodynamic properties of the mixture; in particular, the excess mixing properties and their derivatives. The cross interaction parameters, 12 and ε12 , are given by 11 + 22 ≡ (1 + ı ), 2
(3)
ε12 = (1 + ıε )(ε11 ε22 )1/2 ≡ (1 + ıε )ε.
(4)
12 = (1 + ı ) and
If both ı and ıε are set to zero, the LB rule is recovered. The size and energy parameters, and ε, are used henceforth to scale distances and energies, respectively, and we set therefore = 1 and ε/kB = 1, where kB is the Boltzmann constant. Excess property Y of mixing is given, in general, by Y = Ymix (P, T ) −
Fig. 1. The considered state points shown in the phase diagram of the pure LennardJones fluid. The triple and critical points are marked as tp and cp, respectively. Table 1 Thermodynamic parameters of the three simulated state points. State point
A
B
C
T∗ P∗
0.7000 0.0017
1.10 0.08
1.45 1.00
the number of moles of component i. In the case of the considered LJ mixture we get Y = Ymix (P, T ) − YLJ (P, T ).
(6)
Whereas the total volume is obtained directly during the simulation, enthalpy, H, is evaluated from the relation H U P = + − 1, NkB T NkB T kB T
2. Basic definitions and computational details
12
43
(0)
ni yi (P, T ),
(5) (0)
where Ymix is the total measured property of the mixture, yi is the molar property of pure compound i [at the same (P, T )], and ni is
(7)
where U is the internal energy and is the number density, = N/V . Partial molar quantities yi are given, in general, by
yi =
∂Ymix ∂ni
.
(8)
/ ni T,P,nj =
In the case of binary mixtures, the partial molar function of either component, let us say 1, can be expressed by means of the excess function as
(0)
y1 = y1 +
∂Y ∂n1
(0)
= y1 + y − x2 T,P,n2
∂y ∂x2
,
(9)
T,P
where xi denotes the mole fraction of component i, and y = Y/(n1 + n2 ). A common way to handle the excess properties is by smoothing the data by means of the Redlich–Kister parametrization [16], y = x1 x2
Ai (x2 − x1 )i ,
(10)
where the parameters Ai are obtained by the least-square minimization. However, one cannot be sure that the derivative obtained from this fitted polynomial will provide an accurate description of the derivative of the measured data. An alternative is to use a more general method proposed recently by Lubansky et al. [15]. The method overcomes the problem of numerical differentiation of noisy (or unreliable) experimental results by converting the problem into the solution of an integral equation of the first kind, using Tikhonov regularization, that provides the second derivative at a large number of uniformly spaced points; the required first
44
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
Fig. 2. The excess volumes and enthalpies for state points A (upper row), B (middle row), and C (bottom row), and ı = −0.2 (full line), ı = +0.2 (dashed line), and ı = −0.2 (circles), ı = 0.0 (triangles), and ı = +0.2 (squares).
derivative is then obtained by their integration. An advantage of Tikhonov regularization is that it does not require any assumption about the functional form of the data. Furthermore, by using the built-in regularization parameter the noise amplification is kept under control. (For further details we refer the reader to the original papers [15,19].) We used three different values of to verify that the minimum in the partial molar volume, if there is any, is not caused only by noise in the data but that it reflects the real behavior of the data. The values used were opt , 10opt , and 100opt , where opt was the optimal value found by using the general cross vali-
dation technique, as suggested in [15]. A higher value of means a smoother derivative and less effect of the noise of the data. To confirm correctness of the results, we have also checked the obtained results at selected points by comparing them with those obtained using the common Redlich–Kister parametrization. The studied systems consisted of N = 1372 particles in a cubic simulation box. The LJ potential cutoff was set to rc = 5 and this non-shifted cutoff potential with the usual long-range corrections was used. Every MC sweep consisted of 3N attempted displacements followed by one attempt to change the volume of the system.
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
45
Fig. 3. The partial molar volumes at state point B.
The simulation parameters were set so as to maintain, approximately, the acceptance ratio of the displacements about 0.3 and that of the volume changes about 0.5. Both negative and positive values of ıε and ı were considered for concentrations covering the entire range (0, 1). Significant deviations of the combining rules from those of the LB ones may push the mixture into a metastable or even solid phase. To be sure that the reported results do correspond to a liquid mixture, various control quantities were monitored during the simulation to keep the development of the systems under control [20].
3. Results and discussion The mixtures were considered at three state points (see Table 1): one at an ambient-like condition near the triple point (A), one at a subcritical temperature near the critical point (B), and one at a supercritical temperature far away from the critical point (C), in order to get an impression of the overall effects of the combining rules and thermodynamic conditions on the mixing properties. These state points are indicated in the LJ phase diagram shown in Fig. 1 for better comprehension. Concentrations covered the range
46
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
Fig. 4. The pair correlation functions at state point B for ı = −0.2 (left column) and ı = +0.2 (right column) and different ı and concentrations (pure fluid, full line; x = 0.25, short dashed line; x = 0.5, long dashed line; x = 0.75, dotted line).
from x ∈ 0.0 − 0.5 (because of the symmetry of the mixture, results for the other half of the range are identical) with a step of 0.05. As shown in our previous paper [14], the effects of ıε and ı on the excess volume are quite predictable, intuitively, so that only the limiting values of the largest deviations from the LB rule, ı = ±0.2, have been preferably considered and are reported in this paper. The excess volumes and enthalpies of mixing are shown in Fig. 2, and the partial molar volumes are shown in Fig. 3. As regards the excess volume, consequences of the deviations are rather straightforward, and the trends found previously for
state point A are valid also at other conditions. The sign of V is determined by the sign of ı , where a negative deviation leads to a negative excess volume and vice versa. As for the effect of ıε , positive values of ıε give, in general, negative excess properties and vice versa, but the overall changes are only numerical and not significant. A bit different and more complex situation exists for the excess enthalpy. In this case we face a competition between energetic and excluded volume effects and can observe certain qualitative changes due to changes in ı . Examining H in dependence on ı
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
47
Fig. 5. Temperature dependence of the partial molar volume for ı = −0.2 and ı = 0.0 at the considered state points A, B, and C.
for a given ı we find that while for ı > 0 the excess enthalpy is either a concave or convex function of concentration, negative value of ı leads to the change of the curvature to an S-shaped function. As a consequence of this change we may get the excess enthalpy with three extrema (over the entire composition range), i.e., the well-known feature observed experimentally, for instance, for the benzonitrile–toluene mixture [21]. Partial molar volumes at state B are shown in Fig. 3 in dependence on ı for ı = ±0.2. Regardless of the changes in the energetic interaction, for ı < 0 we observe a minimum in v whereas for ı > 0 the partial molar volume is a monotonously decreasing function. These changes in the thermodynamic behavior must result from certain changes in the structure. To detect them, we have therefore computed also the structure factor and pair correlation functions, gij . Because of the symmetry of the components, the A–B pairs are always exposed to the same environment and the most interesting functions are thus the correlation functions between the like species. Their behavior in dependence on concentration is shown in Fig. 4 for state point B; results obtained at state points A and C are qualitatively similar and will be shortly discussed below. In this figure we compare gAA functions for two different combinations of ı parameters: one, for which the minimum in the partial molar volume occurs, and the other with the ‘common’ behavior of the partial molar volume. The former one results from setting ı = −0.2 and the latter to +0.2. Different values of ı modify only numerically the curves but do not change their character. The difference between the two sets is discernible. As we see, in the former case a considerable change in the course of gAA is found in the second coordination shell. The curve gets deformed around the first minimum where a sort of a secondary peak start to evolve. This is somewhat understandable because ı < 0 and particles may get closer one to another. This feature may resemble formation of an amorphous structure but this is not the case. The formation of the amorphous structure is characterized by development of a bump on the rising part of gAA [22] and may be also detectable on the structure factor by appearance of additional peaks. No such behavior has been observed. What we observe is first a certain restructuring in the second coordination shell followed then by a shift of the second peak to closer separations. For the other set of ı parameters the shape of gAA remains LJ-like. The above finding may be quite surprising. The binary mixture of non-additive hard spheres exhibits demixing for ı > 0 and a similar behavior one could expect also for the LJ mixture with this ı and not for ı < 0. An explanation why this is not the case is rather simple: it is a competition between the energy and excluded volume effects. To be more specific, in Fig. 5 we show the partial
molar volume for states A, B, and C, and ı = −0.2 and ı = 0.0. At the lowest temperature, state A, a rather shallow minimum in v occurs but gets more pronounced in state B. When temperature is further increased the minimum gets milder again and at very high temperatures the hard sphere interaction dominates and the minimum disappears. In the case of ı > 0 and temperatures considered, the hard sphere interaction is not dominating and it becomes important only upon further increase of temperature. 4. Conclusions It has been known for decades that the mixing properties cannot be estimated purely on the basis of the pure fluid properties and LB rules. However, when seeking better numerical agreement with experimental data of excess functions, it has been the cross energy interaction that has been primarily manipulated with the excluded volume parameters being kept, in most cases, intact [11,23]. Studying binary mixtures of identical LJ fluids with different cross interaction parameters we have found that the manipulation with the cross energy parameter is important for a fine tuning of results but this does not give rise to a qualitative change in the behavior of the mixture. It is the excluded volume, i.e., deviations from the Lorentz rule, which affects the results most and which may also cause qualitative changes. After submitting this paper, there appeared on line a paper by Boda and Henderson [24] reporting the same conclusions based on a similar study of the LJ mixture. As regards an experimental support for deviations from the Lorentz rule for the geometrical parameter, we are able to refer only to a case study conducted by Schnabel et al. [11] who investigated the effect of the combining rules on the two phase equilibrium properties of CO+C2H6 and N2+C3H6 binary mixtures and found that they strongly depend on both the size and the energy parameters. The above finding seems to apply to all thermodynamic functions. We have focused on the partial molar volume whose behavior for aqueous solutions of alcohols was only recently explained on the molecular level by diffraction experiments [25,26]. Using molecular simulations, this feature is obtained only after incorporation of polarizability [27]; pair-wise potential models have not been able to reproduce this phenomenon [28]. A similar behavior we have observed in this paper for a binary mixture of identical LennardJones fluids for a certain combination of the cross parameters. There may be thus an alternative to polarizability to change qualitatively mixture properties. These findings may help to better understand the behavior of mixtures on the molecular level and to develop combining rules for different classes of mixtures for practical applications.
48
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
Acknowledgements This research was supported by the Czech National Research Program “Information Society”, project No. 1ET400720409, and by the Grant Agency of the Academy of Sciences of the Czech Republic, project No. IAA400720802. References [1] W.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. (1984) 6638–6646. [2] W.L. Jorgensen, J. Phys. Chem. 90 (1986) 1276–1284. [3] W. Damm, A. Frontera, J. TiradoRives, W.L. Jorgensen, J. Compt. Chem. 18 (1997) 1955–1970. [4] W.L. Jorgensen, J. Chem. Phys. (1982) 4156–4163. [5] J.R. Errington, A.Z. Panagiotopoulos, J. Phys. Chem. B 102 (1998) 7470–7475. [6] E.A. Mason, T.H. Spurling, The Virial Equation of State, Pergamon, Oxford, 1969. [7] D.W. Calvin, T.M. Reed III, J. Chem. Phys. 54 (1971) 3733–3738. [8] K.R. Hall, G.A. Iglesias-Silva, G.A. Mansoori, Fluid Phase Equilib. 91 (1993) 67–76. [9] J. Delhommelle, P. Millie, Mol. Phys. 99 (2001) 619–625. [10] F.J. Blas, I. Fujihara, Mol. Phys. 100 (2002) 2823–2838. [11] T. Schnabel, J. Vrabec, H. Hasse, J. Mol. Liq. 135 (2007) 170–178.
[12] J.N.C. Lopes, Mol. Phys. 96 (1999) 1649–1658. [13] D.M. Duh, D. Henderson, R.L. Rowley, Mol. Phys. 91 (1997) 1143–1147. [14] M. Rouha, F. Mouˇcka, I. Nezbeda, Collect. Czech. Chem. Commun. 73 (2008) 533–540. [15] A.S. Lubansky, Y.L. Yeow, Y.-K. Leong, S.R. Wickramasinghe, B. Han, AIChE J. 52 (2006) 323–332. [16] O. Redlich, A.T. Kister, Ind. Eng. Chem. 40 (1948) 345–348. [17] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [18] D. Frenkel, B. Smit, Understanding Molecular Simulation. From Algorithms to Applications, Academic Press, London, 1996. [19] Y.L. Yeow, Y.-K. Leong, J. Solution Chem. 36 (2007) 1047–1061. [20] I. Nezbeda, J. Kolafa, Mol. Simul. 14 (1995) 153–163. [21] J.S. Rowlinson, F.L. Swinton, Liquids and Liquid Mixtures, Butterworths, London, 1982. [22] I. Nezbeda, W.R. Smith, Mol. Phys. 75 (1992) 789–803. [23] T. Schnabel, J. Vrabec, H. Hasse, Fluid Phase Equil. 263 (2008) 144–159. [24] D. Boda, D. Henderson, Mol. Phys. (2008). [25] S. Dixit, A.K. Soper, J.L. Finney, J. Crain, Europhys. Lett. 59 (2002) 377–383. [26] A. Soper, Invited talk at the 15th International Conference on the Properties of Water and Steam, Berlin, Germany, September 7–11, 2008. [27] F. Mouˇcka, I. Nezbeda, Coll. Czech. Chem. Commun., in press. [28] D. Gonzales-Salgado, I. Nezbeda, Fluid Phase Equilib. 240 (2006) 161–166.
Contents lists available at ScienceDirect
Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid
Non-Lorentz–Berthelot Lennard-Jones mixtures: A systematic study Michael Rouha a , Ivo Nezbeda a,b,∗ a b
Faculty of Science, J. E. Purkinje University, 400 96 Ústí nad Labem, Czech Republic E. Hála Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, Academy of Sciences, 165 02 Prague 6 – Suchdol, Czech Republic
a r t i c l e
i n f o
Article history: Received 21 October 2008 Received in revised form 5 November 2008 Accepted 6 November 2008 Available online 18 November 2008 Keywords: Binary mixtures Lennard-Jones fluids Non-Lorentz–Berthelot rules Excess mixing properties Partial molar properties Tikhonov regularization
a b s t r a c t Binary mixtures of two identical Lennard-Jones fluids with non-Lorentz–Berthelot combining rules have been simulated over the entire concentration range at three state points in order to examine the effect of cross interactions on the mixing properties, excess volumes and enthalpies, and partial molar volumes, and on the structure. Various combinations of deviations, in both the energy and size cross parameters, from the Lorentz–Berthelot rule, have been considered and the results analyzed using a recently proposed method based on the Tikhonov regularization. It is found that the most important is the size cross parameter and that by varying the combining rules a variety of qualitatively different behavior can be produced including, e.g., a minimum in the partial molar volumes observed otherwise only in complex mixtures of associating fluids; occurrence of this minimum seems to be associated with changes in the second coordination shell as witnessed by the pair correlation function. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Knowledge of the thermodynamic behavior of fluid mixtures provides not only the fundamental prerequisite for design and operation of industrial processes but also fundamental information on the intermolecular interactions. According to the contemporary level of understanding and modeling of molecules and their interactions, molecules are viewed as bodies with two types of interaction sites: neutral ones, typically interacting as Lennard-Jones (LJ) or exp-6 sites, and the Coulombic sites [1–5]. Parameters of these sites are then determined by their fitting to a set of selected experimental data of pure fluids. During this parametrization, interactions between the unlike sites are estimated by means of certain combining rules. Whereas there is no problem with the interaction between the Coulombic sites, the interaction between, and hence the combining rules for the unlike neutral sites, is rather problematic. In general, the geometric mean rule (the Bethelot rule) is used for the energy parameter (the depth of the potential well), and the arithmetic rule (the Lorentz rule) for the length-scaling (excluded volume) parameter. Although there is only a very vague theoretical justification of the Lorentz–Berthelot (LB) rule [6–8], due to the adjustment of the results to experimental data its inadequacy does not cause principal problems in developing pure fluid interaction models. However,
∗ Corresponding author at: Faculty of Science, J. E. Purkinje University, 400 96 Ústí nad Labem, Czech Republic. Tel.: +420 220390 296; fax: +420 220920 661. E-mail address: [email protected] (I. Nezbeda). 0378-3812/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2008.11.007
the situation is quite different when dealing with mixtures. In this case the developed pure fluid interaction models are used to predict the properties of mixtures and it is evident that the entire concept of combining rules is, in general, ill-footed: even the used pure fluid potentials cannot adequately describe the interaction between the like species in the mixture because they do not account for the presence of the other species. It has been thus well known for decades that the mixture properties are very sensitive to the cross interactions but this is only a consequence of the entire concept. Since a re-parametrization of the pure fluid potentials with respect to more realistic cross interactions seems quite an unrealistic goal, an improvement of the description of mixtures rests on an improvement of the combining rules. One such a straightforward, and physically fully justified possibility provides inclusion of polarizability which, however, would cause nearly unsurmountable problems for theory and makes also molecular simulations very complex and hence very time consuming. Simple physical considerations show that the main problem lies in the cross LJ interactions and no wonder therefore that a number of papers on the LJ mixtures and combining rules have thus been published (for references see, e.g., the recent papers [9–11] and references therein). Nonetheless, the fundamental question whether deviations from the LB rule may produce qualitatively different mixture properties (and what type) does not seem to have been addressed at all. Typically, either certain modifications of the LB rule have been considered and examined or even different rules have been tested for a certain class of mixtures focusing only on mixing or equilibrium properties and seeking always only a better numerical agreement [9–12]. To our best knowledge we are thus not aware of any fundamental
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
and systematic study on the effect of the combining rules on the qualitative changes in the behavior of mixtures. To study effects of deviations of the cross interactions from those given by the LB rule, Duh et al. [13] considered some time ago two identical LJ fluids at a supercritical temperature. However, they considered in simulations only two compositions and only one set of the cross parameters. Adopting the same approach, i.e., to consider the mixture of two identical LJ fluids, we carried out recently a feasibility study considering a large number of non-Lorentz–Berthelot (nLB) cross parameters but only at one state point [14]. Furthermore, we focused on and reported results for excess properties only. Nonetheless, the shape of the excess volume for certain combinations of the nLB cross parameters indicated that the mixture may also exhibit qualitative changes in its behavior, e.g., exhibiting the minimum in the partial molar volume instead of a monotonous dependence on the concentration. We have therefore continued in investigations of mixtures along the same line focusing on the partial molar quantities and structure. Since the evaluation of the partial molar quantities requires a large number of accurate data on excess functions in dependence on concentrations (and hence lengthy time consuming simulations), in the present paper we consider only selected combinations of the cross parameters looking particularly for qualitative changes in the thermodynamic behavior. These changes must result from and reflect changes in the structure (in the case of low concentrations changes in the solvation shell) and we therefore compute also the structure factors and correlation functions which may reveal them. We also consider three different state points to find the impact of thermodynamic conditions (if there is any) on the behavior of the mixture. Finally, from the numerical point of view, there is a novelty in elaboration of the results: to analyze the obtained simulation data and evaluate derivatives (i.e., partial molar quantities) from noisy experimental data we use (and test) a recently proposed method based on the Tikhonov regularization [15] instead of the common Redlich–Kister parametrization [16].
Standard NPT Monte Carlo (MC) simulations [17,18] on binary mixtures of LJ fluids, uLJ (r) = 4εij
ij
r
−
6 ij
r
,
(1)
were performed. Both components of the mixture were identical, 11 = 22 ≡ ;
ε11 = ε22 ≡ ε,
(2)
but the cross interactions were varied in order to find out how they affect the thermodynamic properties of the mixture; in particular, the excess mixing properties and their derivatives. The cross interaction parameters, 12 and ε12 , are given by 11 + 22 ≡ (1 + ı ), 2
(3)
ε12 = (1 + ıε )(ε11 ε22 )1/2 ≡ (1 + ıε )ε.
(4)
12 = (1 + ı ) and
If both ı and ıε are set to zero, the LB rule is recovered. The size and energy parameters, and ε, are used henceforth to scale distances and energies, respectively, and we set therefore = 1 and ε/kB = 1, where kB is the Boltzmann constant. Excess property Y of mixing is given, in general, by Y = Ymix (P, T ) −
Fig. 1. The considered state points shown in the phase diagram of the pure LennardJones fluid. The triple and critical points are marked as tp and cp, respectively. Table 1 Thermodynamic parameters of the three simulated state points. State point
A
B
C
T∗ P∗
0.7000 0.0017
1.10 0.08
1.45 1.00
the number of moles of component i. In the case of the considered LJ mixture we get Y = Ymix (P, T ) − YLJ (P, T ).
(6)
Whereas the total volume is obtained directly during the simulation, enthalpy, H, is evaluated from the relation H U P = + − 1, NkB T NkB T kB T
2. Basic definitions and computational details
12
43
(0)
ni yi (P, T ),
(5) (0)
where Ymix is the total measured property of the mixture, yi is the molar property of pure compound i [at the same (P, T )], and ni is
(7)
where U is the internal energy and is the number density, = N/V . Partial molar quantities yi are given, in general, by
yi =
∂Ymix ∂ni
.
(8)
/ ni T,P,nj =
In the case of binary mixtures, the partial molar function of either component, let us say 1, can be expressed by means of the excess function as
(0)
y1 = y1 +
∂Y ∂n1
(0)
= y1 + y − x2 T,P,n2
∂y ∂x2
,
(9)
T,P
where xi denotes the mole fraction of component i, and y = Y/(n1 + n2 ). A common way to handle the excess properties is by smoothing the data by means of the Redlich–Kister parametrization [16], y = x1 x2
Ai (x2 − x1 )i ,
(10)
where the parameters Ai are obtained by the least-square minimization. However, one cannot be sure that the derivative obtained from this fitted polynomial will provide an accurate description of the derivative of the measured data. An alternative is to use a more general method proposed recently by Lubansky et al. [15]. The method overcomes the problem of numerical differentiation of noisy (or unreliable) experimental results by converting the problem into the solution of an integral equation of the first kind, using Tikhonov regularization, that provides the second derivative at a large number of uniformly spaced points; the required first
44
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
Fig. 2. The excess volumes and enthalpies for state points A (upper row), B (middle row), and C (bottom row), and ı = −0.2 (full line), ı = +0.2 (dashed line), and ı = −0.2 (circles), ı = 0.0 (triangles), and ı = +0.2 (squares).
derivative is then obtained by their integration. An advantage of Tikhonov regularization is that it does not require any assumption about the functional form of the data. Furthermore, by using the built-in regularization parameter the noise amplification is kept under control. (For further details we refer the reader to the original papers [15,19].) We used three different values of to verify that the minimum in the partial molar volume, if there is any, is not caused only by noise in the data but that it reflects the real behavior of the data. The values used were opt , 10opt , and 100opt , where opt was the optimal value found by using the general cross vali-
dation technique, as suggested in [15]. A higher value of means a smoother derivative and less effect of the noise of the data. To confirm correctness of the results, we have also checked the obtained results at selected points by comparing them with those obtained using the common Redlich–Kister parametrization. The studied systems consisted of N = 1372 particles in a cubic simulation box. The LJ potential cutoff was set to rc = 5 and this non-shifted cutoff potential with the usual long-range corrections was used. Every MC sweep consisted of 3N attempted displacements followed by one attempt to change the volume of the system.
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
45
Fig. 3. The partial molar volumes at state point B.
The simulation parameters were set so as to maintain, approximately, the acceptance ratio of the displacements about 0.3 and that of the volume changes about 0.5. Both negative and positive values of ıε and ı were considered for concentrations covering the entire range (0, 1). Significant deviations of the combining rules from those of the LB ones may push the mixture into a metastable or even solid phase. To be sure that the reported results do correspond to a liquid mixture, various control quantities were monitored during the simulation to keep the development of the systems under control [20].
3. Results and discussion The mixtures were considered at three state points (see Table 1): one at an ambient-like condition near the triple point (A), one at a subcritical temperature near the critical point (B), and one at a supercritical temperature far away from the critical point (C), in order to get an impression of the overall effects of the combining rules and thermodynamic conditions on the mixing properties. These state points are indicated in the LJ phase diagram shown in Fig. 1 for better comprehension. Concentrations covered the range
46
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
Fig. 4. The pair correlation functions at state point B for ı = −0.2 (left column) and ı = +0.2 (right column) and different ı and concentrations (pure fluid, full line; x = 0.25, short dashed line; x = 0.5, long dashed line; x = 0.75, dotted line).
from x ∈ 0.0 − 0.5 (because of the symmetry of the mixture, results for the other half of the range are identical) with a step of 0.05. As shown in our previous paper [14], the effects of ıε and ı on the excess volume are quite predictable, intuitively, so that only the limiting values of the largest deviations from the LB rule, ı = ±0.2, have been preferably considered and are reported in this paper. The excess volumes and enthalpies of mixing are shown in Fig. 2, and the partial molar volumes are shown in Fig. 3. As regards the excess volume, consequences of the deviations are rather straightforward, and the trends found previously for
state point A are valid also at other conditions. The sign of V is determined by the sign of ı , where a negative deviation leads to a negative excess volume and vice versa. As for the effect of ıε , positive values of ıε give, in general, negative excess properties and vice versa, but the overall changes are only numerical and not significant. A bit different and more complex situation exists for the excess enthalpy. In this case we face a competition between energetic and excluded volume effects and can observe certain qualitative changes due to changes in ı . Examining H in dependence on ı
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
47
Fig. 5. Temperature dependence of the partial molar volume for ı = −0.2 and ı = 0.0 at the considered state points A, B, and C.
for a given ı we find that while for ı > 0 the excess enthalpy is either a concave or convex function of concentration, negative value of ı leads to the change of the curvature to an S-shaped function. As a consequence of this change we may get the excess enthalpy with three extrema (over the entire composition range), i.e., the well-known feature observed experimentally, for instance, for the benzonitrile–toluene mixture [21]. Partial molar volumes at state B are shown in Fig. 3 in dependence on ı for ı = ±0.2. Regardless of the changes in the energetic interaction, for ı < 0 we observe a minimum in v whereas for ı > 0 the partial molar volume is a monotonously decreasing function. These changes in the thermodynamic behavior must result from certain changes in the structure. To detect them, we have therefore computed also the structure factor and pair correlation functions, gij . Because of the symmetry of the components, the A–B pairs are always exposed to the same environment and the most interesting functions are thus the correlation functions between the like species. Their behavior in dependence on concentration is shown in Fig. 4 for state point B; results obtained at state points A and C are qualitatively similar and will be shortly discussed below. In this figure we compare gAA functions for two different combinations of ı parameters: one, for which the minimum in the partial molar volume occurs, and the other with the ‘common’ behavior of the partial molar volume. The former one results from setting ı = −0.2 and the latter to +0.2. Different values of ı modify only numerically the curves but do not change their character. The difference between the two sets is discernible. As we see, in the former case a considerable change in the course of gAA is found in the second coordination shell. The curve gets deformed around the first minimum where a sort of a secondary peak start to evolve. This is somewhat understandable because ı < 0 and particles may get closer one to another. This feature may resemble formation of an amorphous structure but this is not the case. The formation of the amorphous structure is characterized by development of a bump on the rising part of gAA [22] and may be also detectable on the structure factor by appearance of additional peaks. No such behavior has been observed. What we observe is first a certain restructuring in the second coordination shell followed then by a shift of the second peak to closer separations. For the other set of ı parameters the shape of gAA remains LJ-like. The above finding may be quite surprising. The binary mixture of non-additive hard spheres exhibits demixing for ı > 0 and a similar behavior one could expect also for the LJ mixture with this ı and not for ı < 0. An explanation why this is not the case is rather simple: it is a competition between the energy and excluded volume effects. To be more specific, in Fig. 5 we show the partial
molar volume for states A, B, and C, and ı = −0.2 and ı = 0.0. At the lowest temperature, state A, a rather shallow minimum in v occurs but gets more pronounced in state B. When temperature is further increased the minimum gets milder again and at very high temperatures the hard sphere interaction dominates and the minimum disappears. In the case of ı > 0 and temperatures considered, the hard sphere interaction is not dominating and it becomes important only upon further increase of temperature. 4. Conclusions It has been known for decades that the mixing properties cannot be estimated purely on the basis of the pure fluid properties and LB rules. However, when seeking better numerical agreement with experimental data of excess functions, it has been the cross energy interaction that has been primarily manipulated with the excluded volume parameters being kept, in most cases, intact [11,23]. Studying binary mixtures of identical LJ fluids with different cross interaction parameters we have found that the manipulation with the cross energy parameter is important for a fine tuning of results but this does not give rise to a qualitative change in the behavior of the mixture. It is the excluded volume, i.e., deviations from the Lorentz rule, which affects the results most and which may also cause qualitative changes. After submitting this paper, there appeared on line a paper by Boda and Henderson [24] reporting the same conclusions based on a similar study of the LJ mixture. As regards an experimental support for deviations from the Lorentz rule for the geometrical parameter, we are able to refer only to a case study conducted by Schnabel et al. [11] who investigated the effect of the combining rules on the two phase equilibrium properties of CO+C2H6 and N2+C3H6 binary mixtures and found that they strongly depend on both the size and the energy parameters. The above finding seems to apply to all thermodynamic functions. We have focused on the partial molar volume whose behavior for aqueous solutions of alcohols was only recently explained on the molecular level by diffraction experiments [25,26]. Using molecular simulations, this feature is obtained only after incorporation of polarizability [27]; pair-wise potential models have not been able to reproduce this phenomenon [28]. A similar behavior we have observed in this paper for a binary mixture of identical LennardJones fluids for a certain combination of the cross parameters. There may be thus an alternative to polarizability to change qualitatively mixture properties. These findings may help to better understand the behavior of mixtures on the molecular level and to develop combining rules for different classes of mixtures for practical applications.
48
M. Rouha, I. Nezbeda / Fluid Phase Equilibria 277 (2009) 42–48
Acknowledgements This research was supported by the Czech National Research Program “Information Society”, project No. 1ET400720409, and by the Grant Agency of the Academy of Sciences of the Czech Republic, project No. IAA400720802. References [1] W.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. (1984) 6638–6646. [2] W.L. Jorgensen, J. Phys. Chem. 90 (1986) 1276–1284. [3] W. Damm, A. Frontera, J. TiradoRives, W.L. Jorgensen, J. Compt. Chem. 18 (1997) 1955–1970. [4] W.L. Jorgensen, J. Chem. Phys. (1982) 4156–4163. [5] J.R. Errington, A.Z. Panagiotopoulos, J. Phys. Chem. B 102 (1998) 7470–7475. [6] E.A. Mason, T.H. Spurling, The Virial Equation of State, Pergamon, Oxford, 1969. [7] D.W. Calvin, T.M. Reed III, J. Chem. Phys. 54 (1971) 3733–3738. [8] K.R. Hall, G.A. Iglesias-Silva, G.A. Mansoori, Fluid Phase Equilib. 91 (1993) 67–76. [9] J. Delhommelle, P. Millie, Mol. Phys. 99 (2001) 619–625. [10] F.J. Blas, I. Fujihara, Mol. Phys. 100 (2002) 2823–2838. [11] T. Schnabel, J. Vrabec, H. Hasse, J. Mol. Liq. 135 (2007) 170–178.
[12] J.N.C. Lopes, Mol. Phys. 96 (1999) 1649–1658. [13] D.M. Duh, D. Henderson, R.L. Rowley, Mol. Phys. 91 (1997) 1143–1147. [14] M. Rouha, F. Mouˇcka, I. Nezbeda, Collect. Czech. Chem. Commun. 73 (2008) 533–540. [15] A.S. Lubansky, Y.L. Yeow, Y.-K. Leong, S.R. Wickramasinghe, B. Han, AIChE J. 52 (2006) 323–332. [16] O. Redlich, A.T. Kister, Ind. Eng. Chem. 40 (1948) 345–348. [17] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [18] D. Frenkel, B. Smit, Understanding Molecular Simulation. From Algorithms to Applications, Academic Press, London, 1996. [19] Y.L. Yeow, Y.-K. Leong, J. Solution Chem. 36 (2007) 1047–1061. [20] I. Nezbeda, J. Kolafa, Mol. Simul. 14 (1995) 153–163. [21] J.S. Rowlinson, F.L. Swinton, Liquids and Liquid Mixtures, Butterworths, London, 1982. [22] I. Nezbeda, W.R. Smith, Mol. Phys. 75 (1992) 789–803. [23] T. Schnabel, J. Vrabec, H. Hasse, Fluid Phase Equil. 263 (2008) 144–159. [24] D. Boda, D. Henderson, Mol. Phys. (2008). [25] S. Dixit, A.K. Soper, J.L. Finney, J. Crain, Europhys. Lett. 59 (2002) 377–383. [26] A. Soper, Invited talk at the 15th International Conference on the Properties of Water and Steam, Berlin, Germany, September 7–11, 2008. [27] F. Mouˇcka, I. Nezbeda, Coll. Czech. Chem. Commun., in press. [28] D. Gonzales-Salgado, I. Nezbeda, Fluid Phase Equilib. 240 (2006) 161–166.