Interfacial tensions of industrial ﬂuids from a molecular-based square gradient theory

This work reports a procedure for predicting the interfacial tension of pure ﬂuids. It is based on scaling arguments applied to the inﬂuence parameter of the van der Waals theory of inhomogeneous ﬂuids. The molecular model stems from the application of the Square Gradient Theory to the SAFT-VR Mie equation of state. The theory is validated against computer simulation results for homonuclear pearl-necklace linear chains made up to six Mie ( λ − 6) beads with repulsive exponents spanning from λ = 8 to 44 by combining the theory with a corresponding states correlation to determine the intermolecular potential parameters. We provide a predictive tool to determine interfacial tensions for a wide range of molecules including hydrocarbons, ﬂuorocarbons, polar molecules, among others. The proposed methodology is tested against comparable existing correlations in the literature, proving to be vastly superior, exhibiting an average absolute deviation of 2.2 %.


INTRODUCTION
Interfacial tension (IFT) is arguably the key thermophysical property that governs the behaviour of inhomogeneous fluids. Its relevance is rooted in the fact that the magnitude of the IFT and its relationship to other state variables (i.e., temperature, pressure, composition) controls several interfacial phenomena such as wetting transitions, interfaces at the vicinity of critical states, nucleation of new phases, etc. The physical understanding and modelling of IFT also provides a route to link tensions with the inhomogeneous behaviour of fluids at molecular level, such as concentration of species along the interfacial zone, interfacial width, etc. 1,2 An additional distinctive characteristic of the IFT is that its value can be obtained from experimental measurements, 3 molecular simulations, 4 and theoretical approaches. 1,3,5,6 Specifically, experimental determinations can be carried out by using tensiometers. In these devices, the IFT is indirectly measured from the force needed to detach an object from a free surfaces (e.g., Wilhelmy plate and du Noüy ring tensiometers) or by combining Laplace's equation with some characteristic dimensions of the system, such as the liquid height in a capillary tube (e.g., capillary rise tensiometer) or the silhouette of a pendant or ellipsoid drop (e.g., pendant drop and spinning drop tensiometers, respectively). For further discussions related to tensiometers and the experimental techniques the reader is redirected to Refs. 3and7 . From a molecular simulation perspective, inhomogeneous fluids can be simulated in the canonical ensemble by using both Molecular Dynamics and Monte Carlo schemes. 4 In both schemes, the IFT can be computed from the mechanical and/or the thermodynamic route. In the mechanical route or IK method 8 , the IFT is computed from the integration of the difference between the normal and tangential pressure (Hulshof's integral 9 ), which are described by the diagonal components of the Irving and Kirkwood tensor. 10 In the thermodynamic route or Test Area method 11 , the IFT is computed from the change in the Helmholtz energy in the limit of an infinitesimal perturbation in the interfacial area. Empirically, IFT can be related to the difference between the liquid and vapor densities through a phenomenological relationship know as the Parachor. 12,13 Such an approach is useful from a practical standpoint, buy its lack of rigour precludes any meaningful extrapolation. From a more fundamental viewpoint, the calculation of the IFT can be based on corresponding states principles 14,15 and statistical mechanics perturbation methods 6 , where the Square Gradient Theory (SGT) 16,17 stands out as one of the most widely used. From a formal perspective, classical density functional theory (DFT) also provides a route to determine density profiles and IFT in simple scenarios, but it us yet to be fully developed for non-spherical fluids. [18][19][20][21] In SGT, the Helmholtz energy density of the interfacial fluid is described by the sum of two contributions. The first part takes into account the Helmholtz energy density for the homogeneous fluid at a local-density, while the second part represents the inhomogeneous contribution of Helmholtz energy by a product of square local-density gradients and some characteristic parameters. These latter parameters have been historically called influence parameters since their values govern the stability and characteristic length scales of the interfaces. The popularity of SGT can be attributed to its relative simplicity and to the unique proposition of using the same equation of state to model simultaneously the homogeneous Despite the success of SGT for describing interfacial properties of pure fluids and fluid mixtures, this theory depends crucially on the independent determination of the influence parameter. Theoretically, the influence parameter can be computed from its molecular definition (i.e., integration of the direct correlation function of the homogeneous fluid), but the available theories for the two-body direct correlation function between two species in homogeneous fluids are not completely developed, as the results still exhibit poor performance when compared to experimental or molecular dynamics results. [72][73][74][75] To circumvent this problem, Carey 22,76 proposed to invest the problem and back-calculate the influence parameter using experimental data of IFT and SGT to later correlate the results to the EoS parameters. This semi empirical approach has been broadly used for pure fluids and nowadays quite refined correlations are available. For instance, Zuo and Stenby 27 , Miqueu et al. 29 and Lin et al. 34 , have used the Peng -Robinson EoS 77 and its volume translated version in SGT to correlate the influence parameter as a function of temperature and the acentric factor. These correlations have shown a remarkably good agreement between SGT estimations and experimental data. The same procedure has been also used starting from molecular based EoS, such as SAFT EoS and its variants, where both experimental data and Molecular Dynamics simulations have been used to correlated the influence parameter (see for instance Refs. 48,53,71,and78 ). Mixtures add another dimension of complexity, whereas the corresponding binary (cross) influence parameter must then be determined, usually in an empirical fashion through simple geometric mixing rules or by fitting to binary experimental data. This approach seems to work for most simple cases and is trivially extendable to multicomponent mixtures.
In summary, while SGT is a powerful theory for describing the interfacial tension of pure fluids and fluid mixtures its main limitation for it to be used as a predictive theory is the lack of generality and limited transferability of the influence parameters. The main goal of this work is to develop a flexible, transferable and universal set of relations for the influence parameter for pure fluids. A rather long-standing affort has been made to produce a molecular-based equation of state that can faithfully represent in a quantitative fashion the macroscopic thermodynamic properties of fluids with these potentials. The latest version of these theories, the SAFT-Mie equations 79,80 , which will be discussed later, has been sucessfully employed both as a tool for fitting, correlating and subsequent prediction of fluid phase equilibria in a wide range of scenarios (e.g. vapor-liquid equilibria, water-octanol partition coefficients, liquid-liquid equilibria, etc.) for a wide range of industrially relevant fluids including, but not limited to polar fluids, refrigerants, crude oils, polymers, etc. However, possibly, the most interesting feature of these equations is the direct and quantitative link to the underlying potential, such that information gathered experiments can be incorporated into intermolecular potentials of interest here, is to garner experience from the molecular simulation of vapor-liquid interfaces to directly feed into this framework, in order to build a robust and transferable model capable of predicting the properties of an interfacial system from a minimal amount of commonly available experimental information, such as critical constants (see Refs. [65][66][67][81][82][83][84][85][86][87] for a complete discussion).
This paper is organized as follows: We summarize the main working expressions of the SGT and the SAFT-VR Mie EoS in Section II. In Section III we briefly consider the Molecular simulation methodology used in this work. The main results obtained from the corresponding states correlations for interfacial tension and applications for selected fluids are discussed in Section IV. Finally, the main conclusions are summarized in Section V. In fact, according to our records there are more than 150 scientific papers related to SGT and this theory has been used as base for several PhD thesis (see for example Refs. 76,94-98 ).
As SGT has been broadly discussed in the literature, this section only condenses the main working expressions for modeling the interfacial behavior for the case of pure fluids in vapor -liquid equilibrium. The reader is directed Refs. 99and100 and the corresponding PhD thesis for a complete deduction of SGT.
In the SGT, the interfacial density of a pure fluid, ρ (z), varies continuously from the bulk density of a vapor (ρ (z → −∞) = ρ V ) to the bulk density of a liquid (ρ (z → +∞) = ρ L ). In order to describe this continuous evolution, van der Waals proposed to express the Helmholtz energy (A) of an interfacial or inhomogeneous fluid as a second order Taylor expansion about the homogenous Helmholtz energy density, a 0 , at the local density ρ. For the case of pure fluids characterized by flat interfaces between adjacent phases, the Taylor expansion may be performed along the interface width by considering a normal z-coordinate (perpendicular to the plane of the interface) as follows: where S corresponds to the interfacial area and c denotes the influence parameter. In this expression, the first term within the integral refers to the homogeneous fluid contribution and the second term corresponds to the inhomogeneous part expressed as gradient term multiplied by the influence parameter (c). The minimization of Equation (1) for a closed system leads to the following second order differential equation of ρ (z): in Equation (2), Ω represents the grand thermodynamic potential, which is defined as where the superscript denotes that the term is evaluated at phase equilibrium conditions, and a 0 is the molar Helmholtz energy.
Within the SGT (cf. Equation (1)), the interfacial tension, γ, between vapor -liquid phases can be computed from the following expression: Alternatively, the interfacial tension of pure fluids can be also calculated by using the following integral expression, which is obtained by combining Equation (3) and Equation (4): In this work, we select the most update version of SAFT, the SAFT-VR-Mie EoS, 79 which represents a significative advance in SAFT models. 102 The main improvement of this version versus older encarnations is the expression up to third order in the residual Helmholtz energy of the monomer term and the flexibility brought about by being based on the Mie potential 103 , u Mie , which can be represented by In Equation (6), λ r and λ a are the repulsion and attraction parameters of the intermolecular potential, respectively, r ij is the center-to-center distance of the interacting segments, ε is the energy scale corresponding to the potential well depth, σ is the length scale, corresponding loosely with an effective segment diameter, and C is a constant defined as: The Mie potential reverts to the well-known Lennard-Jones model 104 if the repulsive and attractive exponents are taken as 12 and 6, respectively. The expression of the Helmholtz energy density of SAFT-VR Mie EoS for a non associating chain fluid is given by 79 where a = A/(Nk B T ) and A is the total Helmholtz energy, N is the total number of molecules, N av is the Avogadro constant, T is the temperature, k B is the Boltzmann con- , and ρ is the molar density of the fluid. a MONO represents monomer (unbounded) contribution for a chain composed of m s tangential segments, a CHAIN accounts for the formation of chain molecules and a IDEAL is the ideal gas contribution. For a complete overview of this model the reader is referred to Ref. 79 .
Coarse-graining of fluid potential using SAFT In this work, the pure fluids are modelled as freely jointed tangential non associating spheres (pearl-necklace model) characterized by five parameters: m s , λ r , λ a , ε, and σ, which can be found from several routes. One could be tempted to fit these parameters to the properties of a lower resolution model, e.g. a fully atomistic classic molecular model of the acentric factor, and liquid density) and the force field parameters. This latter procedure is followed herein. The number of beads, m s , that describe a molecular model is determined beforehand by observation of the molecule geometry. The underlying model requires the bead to be tangent to each other and in a linear configuration (pearl-necklace model). Ramrattan et al. 107 have shown that there is a conformality relationship between the exponents of the Mie potential; an infinite number of exponent pair (λ a , λ r ) will provide essentially the same field phase behaviour. Following this, we chose to fix the attractive potential λ a = 6 108 leaving the repulsive exponent λ r = λ as the lone parameter that defines the range of the intermolecular potential.
Once the EoS parameters have been fixed, Equation (8) Equation (9) corresponds to the mechanical equilibrium condition (P 0 = P L = P V ), Equation (10) expresses the chemical potential constraint ((∂a 0 /∂ρ) T,V ≡ µ, µ 0 = µ L = µ V ), and Equation (11) is a differential stability condition for interfaces, comparable to the Gibbs energy stability constraint of a single phase. 36 The influence parameter In the original van der Waals theory, the influence parameter, c, is defined as a constant, but modern versions of this theory reflect that this parameter should be a function of the direct correlation function of the homogeneous fluid. According to Bongiorno et al., 73 and Yang et al., 93 the rigorous definition of c is given by the following integral expression: where c 0 (r; ρ) is the direct correlation function of homogeneous fluid and r is a spatial coordinate. Since c 0 (r; ρ) is intractable from an analytic viewpoint, some models have been developed to estimate the influence parameters from other measurable or computable quantities. According to Rowlinson and Widom, 1 one of the most successful approximations for c 0 is to consider c 0 (r; ρ) ≈ c 0 (r; T ) where c 0 (r; T ) can be described by the Percus -Yevick approximation: 109 where g (r) is the radial distribution function of a fluid in the homogeneous state and u (r) is the intermolecular potential, respectively. In a mean field approximation, a locally uniform fluid distribution can be assumed, hence, g (r) ≈ 1. Linearizing Equation (13), a meanspherical approximation for the direct correlation function of homogeneous fluid can be obtained: which if replaced in Equation (12), and considering an isotropic fluid becomes Equation (15) represents the simplest approximate model for c and it acquires a final form once the intermolecular potential is defined. For the case of the Mie potential (cf. Equation (6) and Equation (15)) simplifies to: From Equation (16) it follows that c can be treated as a constant once the intermolecular potential exponents (λ r , λ a ) and fluid parameters (ε, σ) have been defined. Some particular cases of Equation (16) (16) predicts c ≃ 7.181 N 2 av εσ 5 for a Lennard-Jones pure fluid which is 1.6 times the value reported by Duque et al., 48 .
In summary, the mean-spherical approximation for the direct correlation function of homogeneous fluid provides a route to obtain an analytical expression for c. This expression can be used for qualitative description of interfacial properties from SGT, but requires refinement if it is to be used as a predictive tool. In this work, the refinement of the theory is carried out by exploiting the direct link between EoS and the underlying intermolecular enough molecules to ensure a representative bulk phase. To guarantee this, L z was fixed to be much larger than L x and L y . Typically L x = L y = 20σ and L z is set to be 5 to 7 times larger. In all cases, the two interfaces spontaneously appear in the x − y plane.
In order to reduce the potential truncation and system size effects involved in the phase equilibrium and interfacial properties calculations, the cut-off radius (r c ) has been taken equal to a relatively large value of 10σ. It has been shown 53,114,115 that a cut-off above six segment diameters provides a reliable description for the pressure and interfacial properties of the Lennard-Jones fluid and we expect that to translate to Mie fluids. In order to characterize the bulk phase and interfacial behavior, density profiles are calculated by dividing the system in 400 slabs along the z direction. The molecular density profiles, ρ i (z), are obtained by assigning the position of each bead, z i , to the corresponding slab, and constructing the molecular density from mass balance considerations. Additionally, these profiles are displaced so that the center of mass of the liquid slab lies at the center of the simulation cell, this displacement helps to avoid smearing of the profiles due to fluctuations of the center of mass location. In order to estimate errors on the variables computed, the sub-blocks average method has been applied. 119 In that approach, the production period is divided into n independent blocks. The statistical error is then deduced from the standard deviation of the average divided by n 1/2 . The equilibrium pressure and interfacial tension are obtained using the Irving-Kirkwood method, where the profiles of the pressure tensor diagonal elements are calculated employing the virial expression: 10 where P kk is the pressure tensor elements, the subscript kk represents the spacial coordinate, either x, y, or z, k B is Boltzmann's constant, T is the absolute temperature, S is the interfacial area, N is the number of molecules, f ij is the force on molecule i due to molecule j, and r ij represents the distance between molecules i and j. f ij , r ij contributions have been equally distributed among the slabs corresponding to each molecule and all the slabs between them. The term k B T ρ (z) in Equation (17) takes into account the kinetic contribution, which is proportional to the ideal gas pressure. The term into the ⟨· · · ⟩ brackets corresponds to the configurational part which is evaluated as ensemble averages and not at instantaneous values. From the pressure elements of Equation (17), the vapor pressure can be determined, corresponding to the P zz element, while the interfacial tension, γ, can be calculated as 9 In this expression, the additional factor 1/2 comes from having two interfaces in the system.  The advantage of this combination is that the fluid is described by the same five parameters: m s , λ r , λ a , ε, and σ in both approaches. As explained previously, the attractive exponent can be fixed (λ a = 6) leaving the repulsive exponent λ r = λ as the lone parameter that defines the range of the intermolecular potential, without any loss of generality. 106,107 In order to combine the MD results of IFT (Equation (18)) to SGT (Equation (5)), we postulate that c is a constant for each particular fluid (i.e., c = ℑ (m s , ε, σ, λ, 6)). If so, the influence parameter can be regressed from the integral (Equation (5)). It proves valuable to express the IFT from MD (Equation (18)) and SGT (Equation (5)) in reduced variables: Of course, both equations, Equation (19) and Equation (20), must provide the same numerical results as we compare the exact results from MD to those predicted from the theory.
The integral on the right hand side of Equation (20) can be evaluated without recourse to specifying a particular value of the influence parameter, hence a plot of the tension calculated through MD (γ * MD ) as a function of the´ρ * ,L ρ * ,V 2 (Ω * + P * ,0 )dρ * provides a means of evaluating the prefactor of Equation (20), namely 1 ms c εσ 5 . Since m s , ε and σ are specified a priori, the procedure provides an explicit value of the influence parameter.
Considering a constant value for the attractive exponent (λ a = 6), the value of the repulsion exponent is λ r = λ and the latter expression reduces to: It is interesting to point out the differences between the estimation of influence parameters presented here to the results published using previous versions of SAFT-VR Mie 126,127 coupled with SGT. In this work, we obtain a temperature-independent influence parameter for very soft potentials such as for instance the Mie (8-6) fluid, even when chain length increases (c.f. Figure 1(a)). For the same soft potential, Galliero et al., 53 have obtained a different behaviour using a previous SAFT-VR Mie model. In fact, they reported a strongly temperature dependent influence parameter, especially when the fluid is approaching the critical region. Presumably this thermal dependence of the influence parameter is an artifact product of the inability of previous SAFT models to represent accurately the critical and near-critical region 128 . The success of the current version of the SAFT theory relies on the extraordinary ability of the third-order expansion term proposed by Lafitte et al. 79 to produce a satisfactory description of VLE near the critical region.

Interfacial Properties for Molecular Fluids
The correlation proposed above is validated by comparing the results obtained from the SGT + SAFT-VR Mie EoS to MD simulation results for the same molecular models (i.e., Mie chains of variable repulsive exponent value and number of segments). Specifically, we test the accuracy of Equation (23) for use within the theory for predicting interfacial properties, such as interfacial density, ρ * s = ρ s σ 3 , profile along the interfacial region, z * = z/σ and interfacial tension, γ * = γσ 2 /ε. Figure 2 shows the ρ * s − z * profiles for the case of molecular chain fluids (m s = 2, 3, 4, and 5) interacting by a Mie (10-6) potential. As expected, as the temperature increases, the interfacial region becomes wider, however the take-home message from these figures is the very good agreement between MD and SGT with the influence parameter calculated from Equation (23). In addition to the interfacial profiles in Figure 3 we display the interfacial tension as a function of temperature (γ * − T * ) for molecular chain fluids (m s = 1 to 6) interacting through a Mie (λ, 6) potential with λ = 8, 10, 12, 20. From In the last expression, T br denotes the reduced normal boiling temperature (T br = T b /T c ).
In Equations (24) and (25), the temperature is in Kelvin and the pressure is in bars.
A second popular correlation has been proposed by Curl  (1 − T r ) 11/9 As an extension of the scaling proposed initially by Guggenheim 14 , in the Pitzer expression, ω is the acentric factor, which is related to the deviation between the vapor pressure of a given fluid and that of a noble gas. A further correlation for γ has been proposed by Zuo and Stenby. 131 In this work, the authors interpolate between two well-defined reference fluids.
The final expression for γ is given by the following expressions: In Equation ( In addition to the previous correlations, Miqueu et al., 135 have proposed the following expression: In Equation (30), V c is the critical volume of the fluid and t = 1 − T/T c . This correlation has been successfully applied for petroleum fluids (i.e. light gases, saturated hydrocarbons, aromatics) and polar compounds (i.e. refrigerants).

Interfacial Properties for Industrial Fluids
In this section, we test the proposed correlation for the case of industrial fluids. It is important to note that the approach used here is able to only to calculate the variation of the interfacial tension to the temperature but also to calculate the interfacial density profiles.
As an example, we retake the model for hexane discussed in detail in Ref. 106 (23)) and MD simulations carried out by using the same Mie parameters than the theory. From Figure 4(c), it is seen that these profiles display the expected behavior, (i.e., they decrease monotonically across the interface following s hyperbolic tangent shape, that spans from the liquid to the vapor bulk phase).
Noticeably there is a very good agreement between theory and simulations over a broad temperature range. Figure 4(d) represents the interfacial tension behavior as a function of temperature. From the latter figures, it is evident that here is a remarkable quantitative agreement to the MD results as well as to experimental tensiometry results all the way from low temperatures to the critical temperature. In the Supporting Information, we include a workbook written in Mathematica code that performs all calculations described above of n-hexane.
In order to evaluate the performance of SAFT-VR Mie + SGT and the new expression for the influence parameter (see Equation (23)), we selected some test fluids (e.g. hydrocarbons, N 2 , refrigerants, etc.), and applied a two-step predictive approach. First, the fluids are idealized as chains of CG tangential spheres interacting with each other through a Mie (λ − 6) potential, whose parameters (m s , ε, σ, λ) are obtained by using the corresponding state principia described by Mejía et al. 106 . Basically, in this step once the number of beads in the chain is defined (m s ) by examining the overall molecular geometry (i.e. its length to beadth ratio), the value of λ is calculated from the acentric factor of the fluid (ω), expressing the relationship between the range of the potential and the vapor pressure of the fluid. In an analogous fashion, the energy parameter (ε) is obtained from the critical temperature of the fluid (T c ), and the value of σ is calculated from the liquid density evaluated at 0.7 of T c .
Once the Mie (λ − 6) parameters have been identified, the van der Waals constant, α, can be computed from Equation (22) and the influence parameter is obtained from Equation (23).
A second final step is to calculate the phase equilibrium from Equations (9) to (11) and the corresponding interfacial tension from Equation (5). Table I  Absolute Deviation (%AAD γ) of the calculated interfacial tensions is 2.8 % in the case of hydrocarbons (see Figure 5(a)), and 3.7 % for the other fluids selected (see Figure 5(b)).
In order to evaluate the performance of the proposed methodology to other correlations,  135 . It is seen that the proposed method is not only more accurate than other available correlations, but it is also broader in terms of applicability range. Uniquely, the procedure could also be inverted: given the interfacial tension of a fluid, a set of Mie parameters can be specifically found by a simple analytical approach without the need of performing simulations. This is particularly useful in developing top-down coarse-grain potential for simulations of surfactants and interfacial fluids. This extension is not pursued here but will be the subject of future work.
The proposed molecular model is both robust and transferable, producing both as set of molecular parameters amenable to be used in molecular simulation and a fully consistent theory for predicting the interfacial properties of Mie fluids. It is applicable in as much as a homonuclear non-associating chain remains a good model of the pure fluid; i.e. if is not expected to work for strongly associating fluids (e.g. water and small alcohols). However, where applicable, the predictions show a remarkable accuracy when compared to traditional corresponding state principa correlations. In addition, the proposed approach not only gives an excellent model to predict the interfacial tension over a wide temperature range, even very close to the critical region, but also provides a route to obtain microscopic information of the interfacial region, such as interfacial density profiles.    (stars) density profiles obtained from MD simulations, (squares) results from MD simulations; (filled black circle) recomended experimental data from DECHEMA 136 .

FIGURE 5
Comparison between calculated (lines) with the SAFT-VR Mie+SGT and Equation (23) for the influence parameters and experimental 136 (symbols) interfacial tensions as a function of temperature for various components. Information about the SAFT-VR Mie parameters that were used in the calculations and the experimental data can be found in