A Robust Physics-Based Calculation of Evolving Gas–Liquid Interfaces

Lukáš Šatura; Mária Minichová; Michal Pavelka; Juraj Kosek; Alexandr Zubov

doi:10.1515/jnet-2021-0080

Article Open Access

A Robust Physics-Based Calculation of Evolving Gas–Liquid Interfaces

Lukáš Šatura , Mária Minichová , Michal Pavelka , Juraj Kosek and Alexandr Zubov

Published/Copyright: February 4, 2022

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Journal of Non-Equilibrium Thermodynamics Volume 47 Issue 2

Abstract

Density gradient theory describes the evolution of diffuse interfaces in both mixtures and pure substances by minimization of the total free energy, which consists of a non-convex bulk part and an interfacial part. Minimization of the bulk free energy causes phase separation while building up the interfacial free energy (proportional to the square of gradients of the species’ densities) and it results in the equilibrium shape of the interface. However, direct minimization of the free energy is numerically unstable and the coefficients in the interfacial part of the free energy are often estimated from experimental data (not determined from the underlying physics). In this paper we develop a robust physics-based numerical approach that leads to the interface density profiles for both pure substances and mixtures. The model is free of fitting parameters and validated by available experimental data.

Keywords: density gradient theory; Peng–Robinson equation of state; gas–liquid interface; surface tension

1 Introduction

Understanding the mechanism of particle organization at phase boundaries and surfaces is vital for many branches of industry, medicine, engineering, and science, ranging from separation processes (e. g., liquid–liquid extraction, distillation), bubble flow reactors, shale gas extraction, and production of hetero-phase materials with tailored morphology [1], [2], [3], [4], [5] to the study of age-related human diseases [6].

Probably the first rigorous formulation of diffuse interface dynamics was created by Van der Waals [7], whose theory was later revisited and refined by Landau and Lifshitz [8], Mitsui and Furuichi [9], and many others [10]. Cahn and Hilliard [11] applied the theory to investigate evolution of phase boundary in vapor–liquid and liquid–liquid systems. The objective is to find the equilibrium profiles of densities of particular species ρ i ( r ) by minimizing the total free energy, which consists of non-convex bulk terms and non-local interfacial terms. But what does the total free energy of a system look like?

A classical way is to take an equation of state (EoS), for instance the Van der Waals EoS, and add the squares of the gradients of the densities of particular species multiplied by some constants (influence parameters) which are estimated from experimental data [12], [13], [14], recently even for ternary [15] and quaternary systems [16]. In order to make the theory more predictive, the influence parameters can also be obtained from the first principles. Indeed, the EoS itself comes from a bulk free energy obtained by a local approximation of an interaction potential [17]. When a more precise approximation involving expansion of the interaction potential in space is carried out, the influence parameters can be calculated [18]. In other words, both the bulk and interfacial parts of the free energy can be obtained from the underlying physics. Here we extend this approach to the Peng–Robinson–Stryjek–Vera (PRSV) EoS, improving the predictions for non-polar substances. It should be noted that although the interfacial influence parameters can be rigorously derived from the direct correlation function of the homogeneous fluids [19], [20], [21], mathematically simpler and more practical expressions are usually used.

Minimization of the free energy results in profiles of densities of the species across the phase boundary. In the spatially one-dimensional case both density functional theory (DFT) and density gradient theory (DGT) typically arrive at formulation of a boundary value problem (BVP) with second-order ordinary differential equations accompanied by fixed boundary values (pre-calculated values of density in the bulk of the two distinct phases) [12], [22], [23], [24]. The equations are then discretized, which leads to a set of non-linear algebraic equations [25]. However, the high non-linearity of the problem, unknown length of the phase boundary (and thus of the proper computational domain), stability issues, and non-uniqueness of the desired solution make the problem computationally challenging.

Therefore, iterative schemes [23] adding a time derivative (relaxation) term to the equations have been proposed [26], [27], and the vapor–liquid interface profile is then calculated by solving a one-dimensional Cahn–Hilliard equation. Here, we apply this method to non-polar and polar substances and compare the results with experimental data. Although the relaxation terms come from the mutual friction between the components in the case of mixtures, such reasoning does not apply in the case of pure substances (no friction) because then there is no other component to have friction with. In this paper we provide physical and mathematical reasoning for the relaxation terms in pure substances based on reduction of evolution equations for Korteweg fluids.

Section 2 contains the PRSV EoS and the corresponding free energy, its physics-based non-local extension due to the gradients of densities, the general form of the Cahn–Hilliard equation, and a description of the numerical technique. Section 3 then contains particular results of the simulations and their interpretation. The developed model is validated using available experimental data, most importantly by means of the values of surface/interfacial tension predicted for various systems. The numerical method gives the equilibrium density profiles for non-polar and polar pure substances and mixtures in a robust and physically motivated way.

2 Methodology

2.1 Peng–Robinson–Stryjek–Vera equation of state

The expression for free energy of the vapor–liquid interface is in our work built upon the Peng–Robinson (PR) EoS [28]. This equation predicts liquid densities, vapor pressures, and vapor–liquid equilibrium ratios for pure non-polar chemicals. It is actually a model with three parameters that depend on the critical temperature and pressure, T c and p c , and an acentric factor, ω. The EoS reads

(1) p = R T ρ 1 − b ρ − a α ρ 2 1 − b 2 ρ 2 + 2 b ρ ,

with

(2) a ( T c ) = 0.457235 R 2 T c 2 p c , b ( T c ) = b ( T ) = 0.077796 R T c p c , a ( T ) = a ( T c ) α ( T r , ω ) , α = [ 1 + κ ( 1 − T r ) ] 2 , κ = 0.37464 + 1.54226 ω − 0.26992 ω 2 , and T r = T T c .

Polar compounds need a modification to the PR EoS, the PRSV EoS [29], [30], where κ is redefined as

(3) κ = κ 0 + [ κ 1 + κ 2 ( κ 3 − T r ) ( 1 − T r ) ] ( 1 + T r ) ( 0.7 − T r ) , κ 0 = 0.378893 + 1.4897153 ω − 0.17131848 ω 2 + 0.0196544 ω 3 ,

with κ 1 , κ 2 , and κ 3 being tabulated adjustable parameters for pure compounds satisfying κ = κ 0 for T r ≥ 0.7 while κ = κ 0 for T r ≥ 1.0 for water and alcohols.

The Helmholtz free energy density corresponding to the PR and PRSV EoSs is

(4) f = R T ( ∑ i C ρ i ln ( ρ i − 1 ) − ρ ln ( 1 − b ρ ) ) + a ρ 2 2 b ln ( 1 + ( 1 − 2 ) b ρ 1 + ( 1 + 2 ) b ρ ) ,

where ρ = ∑ i C ρ i is the total molar density and a = ∑ i ∑ j x i x j a i j and b = ∑ i x i b i represent the mixing rules [31], [32]. Note that x i = ρ i / ρ stands for molar fraction of the ith component and the cross term a i j ( i ≠ j ) follows the conventional one-binary-parameter mixing rule from pure-component parameters, defined by equations (2),

(5) a i j = a i a j ( 1 − k i j ) .

The binary parameters, k i j and k j i , which do not depend on the mixture composition, can be obtained from the equality of fugacities on both sides of the interface for each mixture component [30]. Although the EoS describes the bulk thermodynamic properties of the system, it can be used also to estimate the energy stored in the interface, as shown in the following section.

2.2 Weakly non-local free energy and surface tension

Surface tension expresses how free energy is stored in the interface between the phases. The interface is characterized by steep changes in densities while far from the interface the densities approach their bulk values. Assuming a spatially one-dimensional interface of length L along the spatial coordinate z, the total free energy then consists of the bulk contribution and weakly non-local terms (dependent on the first gradients) [17]

(6) F ( ρ ( z , t ) ) = ∫ 0 L ( f ( z , t ) + ∑ i j c i j 2 ∇ ρ i ( z , t ) · ∇ ρ j ( z , t ) ) d z .

The chemical potential of component i is the variational derivative of the free energy with respect to ρ i and has thus also a bulk part and a gradient part,

(7) μ i = μ i b + μ i j ∇ with μ i b = ( ∂ f ∂ ρ i ) T , ρ i ≠ j and μ i j ∇ = − c i j d 2 ρ j d z 2 for i , j = 1 , 2 , … , C .

The influence parameter c i j is a key parameter in the calculation of equilibrium profiles. Physically, it measures the interfacial tension at a given density gradient and affects the steepness of the interfacial density profile. It stems from the DGT developed by Van der Waals and its theoretical expression for pure components ( i = j) reads

(8) c i i a i b i 2 / 3 = d i ,

where d i = 0.305 according to the original DGT formulation [33], or

(9) d i = A i t i + B i ,

(10) A i = − 10 − 16 1.2326 + 1.3757 ω i , B i = 10 − 16 0.9051 + 1.5410 ω i ,

where t i = 1 − T r i , see Miqueu and co-workers [34] for non-polar pure fluids when described by the PR EoS. The cross terms c i j ( i ≠ j ) follow the conventional mixing rule,

(11) c i j = ( 1 − β i j ) c i i c j j for i ≠ j .

The binary parameter β i j and the values of the pure-component influence parameters c i i may be correlated with experimental values of surface/interfacial tension, γ, of a pertinent system of spatial length L via

(12) γ = ∫ 0 L ∑ i C ∑ j C c i j ∇ ρ i · ∇ ρ j d z or γ = ∫ 0 L c ( d ρ d z ) 2 d z

for a single component. This can be seen by transformation of free energy (6) to the grand potential and by subsequent differentiation with respect to the surface area [17]. Solutions to these equations provide boundary conditions for the evolution equations minimizing the free energy, which are discussed in the following section. Nevertheless, the parameter β i j might not be needed in the model at all, provided that the chosen EoS models the mixture properties well [35], [36].

2.3 Cahn–Hilliard equation

Densities of the mixture components are governed by their conservation laws, which lead to the BVPs

(13a) ∂ ρ i ∂ t ( z , t ) = − ∇ · ( ρ i ( z , t ) u i ( z , t ) ) , z ∈ ( 0 , L ) , t > 0 ,

(13b) ∂ ρ i ∂ z ( z , t ) = ∂ μ i ∂ z ( z , t ) = 0 , z = 0 , z = L , t > 0 ,

(13c) ρ i ( z , 0 ) = ρ i init , z ∈ ⟨ 0 , L ⟩ , t = 0 .

In order to close this system of equations, velocities u i need to be defined as products of mobilities M i ( z , t ) and gradients of generalized chemical potentials, which in vector notation reads

(14) u ( z , t ) = M ( z , t ) ∇ μ ( z , t ) R T ,

where M ( z , t ) denotes the mobility tensor [37]. In a single-component system, mobility is equivalent to self-diffusivity, M ( z , t ) ≡ D [11], [38].

For a single-component system with the PR EoS, we obtain the Cahn–Hilliard equation

(15) ∂ ρ ∂ t = D R T ∇ · [ ρ ∇ ( ∂ f ∂ ρ − c ∇ 2 ρ ) ] or ∂ ρ ˜ ∂ τ = 1 R T ∂ ∂ z ˜ [ ρ ˜ ∂ ∂ z ˜ ( ∂ f ∂ ρ − c ρ c L 2 ∂ 2 ρ ˜ ∂ z ˜ 2 ) ]

in the non-dimensional form with ρ ˜ = ρ ρ c , τ = t D L 2 , z ˜ = z L , and with ρ c denoting the molar density of the species in the critical state. Note that the characteristic dimensionless time τ actually encodes the self-diffusivity D. However, this scaling is not feasible in multi-component mixtures because it may lead to disproportion among the magnitudes of densities of the components, reducing the stability of the numerical scheme.

The evolution equations for mixtures (13a) can be derived by reduction of the hydrodynamic equations with momenta of the species among the state variables and the mobilities related to the friction between the species [39]. But what is the physical motivation for a diffusion equation in the case of pure species? The free energy contains a non-local contribution even in the case of pure substances because such contribution is caused by the long-range part of interaction potential among molecules of the substance [17]. Evolution of density and momentum density of the fluid is then described by the Korteweg equations,

(16a) ∂ t ρ = − ∇ · ( ρ u ) ,

(16b) ∂ t ( ρ u ) = − ∇ · ( ρ u ⊗ u ) − ∇ p ( ρ ) + ∇ ( 1 2 c ( ∇ ρ ) 2 + ρ ∇ · ( c ∇ ρ ) ) − ∇ · ( c ∇ ρ ⊗ ∇ ρ ) ︸ = − ρ ∇ F ρ + ∇ · ( ν ρ ( ∇ u + ( ∇ u ) T ) ) ,

where u denotes the fluid velocity, c is the influence parameter, F ρ is the functional derivative of the total free energy with respect to the density, and ν is the kinematic viscosity [39]. In the high Mach number regime, where variations of density are typically observed [40], density in the Korteweg equations satisfies a fast diffusion equation,

(17) ∂ t ρ = 2 ν Δ ρ ,

which can be shown by analysis of convergence of weak solutions [41], [42]. This diffusion equation is called fast because it represents the initial stage of the evolution (flattening the density profile). After the fast diffusion equation has relaxed to its equilibrium, the fluid behaves according to the incompressible Navier–Stokes equations. In other words, the fast diffusion equation can be seen as a consequence of the Korteweg equations in the high Mach number regime.

Instead of the rigorous mathematical analysis, we can also employ a heuristical approach of dynamic maximum entropy reduction or asymptotic analysis [39]. First, however, we need to replace the viscous term in the Korteweg equations with a friction term − ζ u. To motivate such simplification, imagine a tube with no-slip boundary conditions with a viscous fluid passing through. In the stationary state the Hagen–Poiseuille flow is present, where the viscous term is proportional to the average velocity of the fluid. Asymptotic analysis of the momentum Korteweg equation with ζ being large then leads to

(18) − ρ ∇ F ρ = − ζ u .

When this equation is plugged into the equation for density, we obtain a Cahn–Hilliard-type equation for pure substances.

2.4 Numerical solution procedure

Integration of the Cahn–Hilliard equation as an initial value problem was employed numerically via the explicit Euler method. The profile of each component’s molar density field was discretized into N spatial steps within one dimension of length L. In each jth iteration step, the right-hand side term of equation (15) was calculated gradually from the inner terms to the outside, i. e., applying only first and second derivatives numerically, in the form of finite differences. Formulae of these derivatives – the gradient and the Laplace operator – were chosen of the fourth-order accuracy (using four- and five-point formulae, resp.), to account for a possibly coarse grid.

With respect to the boundary conditions, the model was implemented following two alternative methods. In the non-periodic formulation, Dirichlet boundary conditions with fixed bulk densities of individual species pre-calculated from conditions of thermodynamic equilibrium were applied to the boundary grid points, accompanied by fixing the second derivative in the boundary grid points to zero, as this setting enhanced the stability of the calculation dramatically when compared to the use of the standard Neumann conditions (i. e., first derivative equal to zero), cf. equation (13b). In the periodic formulation of the model, periodic boundary conditions were applied. An advantage of using periodic boundary conditions is the necessity of using only the initial conditions to start the simulation from. Remarkably, the thermodynamically consistent system should then naturally converge into domains with densities as calculated in a non-periodic model formulation.

The overall simulation procedure can be summarized in the following points:

In the case of a non-periodic model, the bulk molar densities of individual species need to be pre-calculated for the subsequent use as Dirichlet boundary conditions. Conditions of thermodynamic equilibrium, i. e., equality of pressure and chemical potential of each component in both phases, need to be solved in this context, while in two-component systems, pressure needs to be specified as an additional degree of freedom. The solution was performed in MATLAB R2021a, using the Levenberg–Marquardt algorithm of the non-linear least-squares constrained lsqnonlin solver from the Optimization Toolbox.
The length of the computational domain was in all simulations set to L = 10 nm. In the non-periodic model case, the pre-calculated Dirichlet boundary conditions were combined with fixing the second derivative of the state variables equal to zero at both domain boundaries. As an initial condition in the non-periodic formulation, a linear density profile between the equilibrium boundary points was used, while in the periodic formulation, a noise-like initial profile was imposed, using pseudo-random numbers to generate perturbations in the homogeneous system:
(19) ρ i init = ρ i ( v ) + ρ i ( ℓ ) 2 [ 1 + rand ( 0 , 1 ) × sin ( N z L ) ] , i = 1 , … N .
Numerical integration of the dimensionless equation (15) for single-component systems, or eqs. (13a) in the case of two-component systems, was carried out using an explicit Euler scheme.
From the steady-state density profiles, interfacial tension was finally estimated using equation (12) by employing the trapezoidal rule. The binary parameters β i j in equation (11) are set to zero in our model; therefore no parameter estimation from experimental data is done, rendering the predictions free of any fitting parameters. Apart from the first step, the whole solution procedure is implemented in FORTRAN language to achieve the highest possible computational performance.

3 Results and discussion

The PRSV EoS leads to the equilibrium densities, vapor pressures, and equilibrium compositions of vapors and liquids. Figure 1 shows how calculated densities of three aliphatic alkanes, benzene, water, and ethanol depend on temperature. Table 1 summarizes critical properties of the compounds together with other input parameters of the PRSV EoS. While for non-polar compounds the predicted and measured densities are in excellent agreement, the predictive power of the PRSV EoS is lower for less symmetric (n-decane) or polar (water, ethanol) compounds, as expected [30]. On the other hand, prediction of vapor pressure is accurate for all calculated compounds across the whole temperature range, including the polar molecules.

Figure 1

Experimental data (discrete points) and prediction by the PRSV EoS (curves) for: (a) and (b) bulk liquid densities of three aliphatic alkanes, benzene, water, and ethanol and (c) equilibrium vapor pressure, all as a function of temperature. Measured data for densities are sourced as follows: n-hexane [44], n-octane [45], n-decane [46], benzene [47], water [48], ethanol [49]. Vapor pressures were taken from [46]. The displayed quantities are scaled by their critical values, listed in Table 1 together with other parameters of the PRSV EoS.

The equilibrium densities serve as boundary conditions for the determination of the equilibrium density profiles. Dynamic simulations of the interface density profile evolution according to the scaled Cahn–Hilliard equation (15) were run at several temperatures for seven different compounds. Figure 3 shows an example of the interface evolution for n-hexane. In the lower part of the figure we can see equalization of the chemical potential.

Because the linear initial profile does not reflect a realistic event occurring at the interface, but rather a starting point for the numerical method, next we impose periodic boundary conditions. Figure 2 shows the results with the initial profile given by the average of the two bulk densities biased by a random perturbation (in the order of one tenth of the reduced density throughout the domain). We observe Ostwald ripening, where smaller domains disappear and larger domains (droplets) are formed.

Stability and the convergence rate of the numerical scheme with periodic boundary conditions are significantly lower when compared to the non-periodic model implementation, as noted also by Nauman and He [43]. However, both non-periodic and periodic versions of the model arrive at quantitatively equal shapes of the interface.

Once the equilibrium density profiles are obtained, the value of surface tension can be extracted from them using equation (12). Table 2 summarizes the numerically obtained values of surface tension together with the pertinent temperature ranges T ran (intersections of the temperature ranges for interfacial tension from the literature and ranges recommended for use of the PRSV EoS [29]). Figure 4 then illustrates the accuracy of the prediction of surface tension. For the selected non-polar organic substances, the prediction is in agreement with the experimental data. On the other hand, both PR and PRSV EoSs fail to predict surface tension in polar compounds (except for the nearly critical region). Therefore, another theoretical framework should be used in the case of polar compounds; for instance, the statistical association fluid theory (SAFT) [50] and the extended perturbed-chain SAFT (PC-SAFT) give accurate predictions in polar fluids [51]. However, it is unclear how to derive the influence parameters in the weakly non-local part of free energy from the EoS in the case of the SAFT model.

Table 1

Parameters of the Peng–Robinson–Stryjek–Vera EoS with temperature range T ran of applicability of the values [29], [30].

Compound	T c ( K )	p c ( Pa )	ω	κ 1	κ 2	κ 3	T ran ( K )
nitrogen (N₂)	126.20	3,400,000	0.03726	0.01996	0.3162	0.535	64–126
water	647.29	22,089,750	0.34380	−0.06635	0.0199	0.443	274–623
n-hexane	507.30	3,012,360	0.30075	0.05104	0.8634	0.460	232–503
n-octane	568.76	2,4864,90	0.39822	0.04460	0.6214	0.509	258–563
n-decane	617.50	2,103,490	0.49052	0.04510	0.8549	0.527	310–563
cyclohexane	553.64	4,075,000	0.20877	0.07023	0.6146	0.530	280–553
benzene	562.16	4,898,000	0.20929	0.07019	0.7939	0.523	279–543
ethanol	513.92	6,148,000	0.64439	−0.03374	−2.6846	0.592	293–485

$Figure 2 Evolution of the scaled density ρ / ρ c \rho /{\rho _{c}} and scaled chemical potential μ / ( R T )\mu /(R\hspace{0.1667em}T) profile of n-hexane in scaled time, τ, at T = 286.15 KT=286.15\hspace{0.1667em}\text{K} with periodic boundary conditions. The domain of one dimension, z, scaled by the length of L = 10 − 8 mL={10^{-8}}\hspace{0.1667em}\text{m}, was discretized into N = 300N=300 points. The scaled time is computed as τ = Y k\tau =Y\hspace{0.1667em}k, where Y is the number of iterations. Here, dimensionless time step k = 10 − 9 k={10^{-9}}.$

Figure 2

Evolution of the scaled density ρ / ρ c and scaled chemical potential μ / ( R T ) profile of n-hexane in scaled time, τ, at T = 286.15 K with periodic boundary conditions. The domain of one dimension, z, scaled by the length of L = 10 − 8 m, was discretized into N = 300 points. The scaled time is computed as τ = Y k, where Y is the number of iterations. Here, dimensionless time step k = 10 − 9 .

$Figure 3 Evolution of the scaled density ρ / ρ c \rho /{\rho _{c}} and scaled chemical potential μ / ( R T )\mu /(R\hspace{0.1667em}T) profile of n-hexane in scaled time, τ, at T = 286.15 KT=286.15\hspace{0.1667em}\text{K} with fixed boundary conditions – equilibrium bulk densities. The domain of one dimension, z, scaled by the length of L = 10 − 8 mL={10^{-8}}\hspace{0.1667em}\text{m}, was discretized into N = 500N=500 points. The scaled time is computed as τ = Y k\tau =Y\hspace{0.1667em}k, where Y is the number of iterations. Here, dimensionless time step k = 10 − 10 k={10^{-10}}.$

Figure 3

Evolution of the scaled density ρ / ρ c and scaled chemical potential μ / ( R T ) profile of n-hexane in scaled time, τ, at T = 286.15 K with fixed boundary conditions – equilibrium bulk densities. The domain of one dimension, z, scaled by the length of L = 10 − 8 m, was discretized into N = 500 points. The scaled time is computed as τ = Y k, where Y is the number of iterations. Here, dimensionless time step k = 10 − 10 .

Table 2

Interfacial tension prediction – evaluation. Values predicted at the selected temperature T are compared with those experimentally measured [46], [52]. The temperature range of experimental data T ran is provided, as well as the statistical measure R 2 of agreement between predictions and experimental data over a wide range of temperatures.

Compound	T ( K )	γ T calc ( mN m − 1 )	γ T exp ( mN m − 1 )	T ran ( K )	R 2
n-hexane	313	15.23	16.28	233–502	0.98
n-octane	339	16.70	17.24	258–528	0.98
n-decane	364	17.20	17.31	310–580	0.99
cyclohexane	388	12.07	14.08	282–553	0.96
benzene	279	29.87	30.77	279–549	0.96
water	355	81.60	62.28	276–544	0.67
ethanol	293	40.04	22.36	293–503	0.64

Figure 4

Interfacial tension of the liquid–vapor systems versus temperature – model predictions (curves) versus experimental data (discrete points). Each solid line connects 11 values obtained by integration of the equilibrium density profile predicted via the PR EoS for each compound. Compared data points from literature: n-hexane [52], n-octane, n-decane, water, and benzene [46].

$Figure 5 Binary vapor–liquid system of nitrogen and n-hexane at T = 298.15 KT=298.15\hspace{0.1667em}\mathrm{K} and p = 20 barp=20\hspace{0.1667em}\mathrm{bar}, according to the PRSV EoS. Evolution of density profiles at time step k = 10 − 20 sk={10^{-20}}\hspace{0.1667em}\mathrm{s}, i. e., with total Y = 6 × 10 8 Y=6\times {10^{8}} iteration steps. The degree of discretization was N = 100N=100.$

Figure 5

Binary vapor–liquid system of nitrogen and n-hexane at T = 298.15 K and p = 20 bar, according to the PRSV EoS. Evolution of density profiles at time step k = 10 − 20 s, i. e., with total Y = 6 × 10 8 iteration steps. The degree of discretization was N = 100.

Similarly as in the case of pure substances, we will now analyze binary mixtures, predicting also gas solubility in liquids. We use the non-dimensional Cahn–Hilliard system of equations with mobilities (14) for unitary relative sizes of the molecules and estimated self-diffusivity constants [53]. Figure 5 shows the density profiles of nitrogen and n-hexane with fixed boundary conditions. Although the prediction of the liquid composition via the PR EoS description is not accurate (error up to 20 % in the composition), the dynamic simulation reveals interesting features predicted by the model (Figure 5). The density of n-hexane looks similar as in the case of pure n-hexane in Figure 3, but the density of nitrogen has a peak at the interface. This corresponds to the adsorption of gaseous molecules at the gas–liquid interface, which decreases the surface tension. Similar profiles have already been theoretically analyzed [54], [55], [56] and later obtained from numerical simulations, for instance by Kou and Sun [32] or Mu and co-workers [27], obtaining the adsorption of methane or a surfactant at the interface. For more detailed information about this interesting phenomenon, we refer the reader to the comprehensive review by Stephan and Hasse [57]. Table 3 shows the trend of decreasing surface tension with increasing pressure.

Table 3

Interfacial tension in the vapor–liquid equilibrium in a mixture of n-hexane and nitrogen at constant temperature T = 298.15 K and different pressures p. Values of bulk nitrogen densities calculated for both the liquid and the vapor phase are provided as well.

p ( kPa )	ρ N 2 L , calc ( mol m − 3 )	ρ N 2 V , calc ( mol m − 3 )	γ T calc ( mN m − 1 )
303.975	44.16	114.54	16.43
506.625	75.69	196.51	16.27
810.600	122.97	319.72	16.03
2026.500	312.41	815.11	15.13

4 Conclusions

The DGT model for evolving vapor–liquid interfaces is based on the Cahn–Hilliard equations for both mixtures and pure substances and it allows for non-stationary simulation of phase boundary in a computationally robust and effective way.

To obtain expressions for the free (Helmholtz) energy of the evolving vapor–liquid system, we used the Stryjek–Vera modification of the PR EoS (PRSV EoS). Values of the influence parameters, which express how energy is stored in the presence of gradients of densities, were not fitted, but determined from the physics behind the EoS. The PRSV EoS provides good estimates for the bulk density values of non-polar substances, while for the selected polar species its predictive power is significantly weaker. This can also be seen when estimating the surface tension, which was carried out for seven distinct substances (five non-polar) at 11 different temperatures. In a two-component simulation of an n-hexane–nitrogen system, the model predicts adsorption of gaseous molecules on the gas–liquid interface as well as the reduction of the interfacial tension with increasing pressure.

The evolution equations were solved by the finite difference method with fixed and periodic boundary conditions. The periodic implementation does not require pre-calculated values of the bulk densities, but it is less numerically stable and more computationally demanding. Therefore, we carried out most of the simulations with fixed boundary conditions, which led to fast and reliable predictions of the values of surface tension (approximately 30 seconds on a common laptop).

In the future we would like to employ more advanced EoSs (for instance PC-SAFT), focus on the description of polymer-solvent systems, and develop a parallelized version of the model to enhance the computational efficiency. Moreover, we would like to attract attention to a more rigorous derivation of the Cahn–Hilliard equation from the equations for Korteweg fluids in the case of pure substances. Our conjecture is that the results could be refined by having a fast Cahn–Hilliard equation instead of the fast diffusion equation.

Funding source: Grantová Agentura České Republiky

Award Identifier / Grant number: 19-22834S

Funding statement: The authors acknowledge financial support from the Czech Science Foundation (GAČR), project 19-22834S.

References

[1] X. Liu and D. Zhang, A review of phase behavior simulation of hydrocarbons in confined space: Implications for shale oil and shale gas, J. Nat. Gas Sci. Eng. 63 (2019), 102901.10.1016/j.jngse.2019.102901Search in Google Scholar

[2] A. Nistor, M. Vonka, A. Rygl, M. Voclová, M. Minichová, J. Kosek, et al., Polystyrene Microstructured Foams Formed by Thermally Induced Phase Separation from Cyclohexanol Solution, Macromol. React. Eng. 11 (2017), 1600007.10.1002/mren.201600007Search in Google Scholar

[3] D. Anders and K. Weinberg, A Thermodynamically Consistent Approach to Phase-Separating Viscous Fluids, J. Non-Equilib. Thermodyn. 43 (2018), 185–191.10.1515/jnet-2017-0052Search in Google Scholar

[4] M. Rezakazemi and S. Shirazian, Gas-Liquid Phase Recirculation in Bubble Column Reactors: Development of a Hybrid Model Based on Local CFD – Adaptive Neuro-Fuzzy Inference System (ANFIS), J. Non-Equilib. Thermodyn. 44 (2019), 29–42.10.1515/jnet-2018-0028Search in Google Scholar

[5] L. Zou and X. Zhang, Non-Equilibrium Thermodynamic Analysis of Flash Evapouration Process in Vacuum Ice Making, J. Non-Equilib. Thermodyn. 46 (2020), 139–147.10.1515/jnet-2020-0085Search in Google Scholar

[6] S. Alberti, A. Gladfelter and T. Mittag, Considerations and challenges in studying liquid-liquid phase separation and biomolecular condensates, Cell 176 (2019), 419–434.10.1016/j.cell.2018.12.035Search in Google Scholar PubMed PubMed Central

[7] J. D. van der Waals, Thermodynamische theorie der capillariteit in de onderstelling van continue dichtheidsverandering, Verhandel. Konink. Akad. Weten. Amsterdam 1 (1893), 56.Search in Google Scholar

[8] L. Landau and E. Lifshitz, On the Theory of the Dispersion of Magnetic Permeability in Ferromagnetic Bodies, Phys. Z. Sowjetunion 8 (1935), 153–169.10.1016/B978-0-08-036364-6.50008-9Search in Google Scholar

[9] T. Mitsui and J. Furuichi, Domain Structure of Rochelle Salt and KH 2 PO 4 , Phys. Rev. 90 (1953), 193–202.10.1103/PhysRev.90.193Search in Google Scholar

[10] R. Mauri, Multiphase Flows, Springer-SBM, The Netherlands, 2013.10.1007/978-94-007-5461-4_9Search in Google Scholar

[11] J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, J. Chem. Phys. 28 (1958), 258–267.10.1002/9781118788295.ch4Search in Google Scholar

[12] X. Liang, M. L. Michelsen and G. M. Kontogeorgis, A density gradient theory based method for surface tension calculations, Fluid Phase Equilib. 428 (2016), 153–163.10.1016/j.fluid.2016.06.017Search in Google Scholar

[13] Y. -H. Li, K. H. Dillard and R. L. Robinson, Vapor-Liquid Phase Equilibrium for Carbon Dioxide–n-Hexanr at 40, 80, and 120 °C, J. Chem. Eng. Data 26 (1981), 53–55.10.1021/je00023a018Search in Google Scholar

[14] O. G. Niño-Amézquita and S. Enders, Phase equilibrium and interfacial properties of water + methane mixtures, Fluid Phase Equilib. 407 (2016), 143–151.10.1016/j.fluid.2015.05.005Search in Google Scholar

[15] A. Danzer and S. Enders, Liquid-Liquid Equilibrium and Interfacial Properties of the System Water + Hexylacetate + 1-Hexanol, Chem. Ing. Tech. 91 (2019), 1597–1605.10.1002/cite.201900013Search in Google Scholar

[16] A. Danzer and S. Enders, Thermodynamic Modeling of Time-Dependent Interfacial Properties in Reactive Liquid-Liquid Systems Close to the Critical Point, J. Chem. Eng. Data 65 (2020), 312–318.10.1021/acs.jced.9b00626Search in Google Scholar

[17] L. Landau and E. Lifshitz, Statistical Physics – Part 1, vol. 5, Course of Theoretical Physics, Butterworth-Heinemann, 1980.10.1016/B978-0-08-023039-9.50007-XSearch in Google Scholar

[18] L. M. Pismen, Nonlocal diffuse interface theory of thin films and the moving contact line, Phys. Rev. E 64 (2001), 021603.10.1103/PhysRevE.64.021603Search in Google Scholar

[19] V. Bongiorno, L. E. Scriven and H. T. Davis, Molecular Theory of Fluid Interfaces, J. Colloid Interface Sci. 57 (1976), 462–475.10.1016/0021-9797(76)90225-3Search in Google Scholar

[20] A. Y .M. Yang, P. D. Fleming and J. H. Gibbs, Molecular theory of surface tension, J. Chem. Phys. 64 (1976), 3732–3747.10.1063/1.432687Search in Google Scholar

[21] R. Evans, The nature of the liquid-vapour interfaceand other topics in the statistical mechanics of non-uniform classical fluids, Adv. Phys. 28 (1979), 143–200.10.1080/00018737900101365Search in Google Scholar

[22] J. Wu, Density Functional Theory for Chemical Engineering: From Capillarity to Soft Materials, AIChE J. 52 (2006), 1169–1193.10.1002/aic.10713Search in Google Scholar

[23] J. Gross, A density functional theory for vapour-liquid interfaces using the PCP-SAFT equation of state, J. Chem. Phys. 131 (2009), 204705.10.1063/1.3263124Search in Google Scholar PubMed

[24] J. Mairhofer, B. Xiao and J. Gross, A classical density functional theory for vapour-liquid interfaces consistent with the heterosegmented group-contribution perturbed-chain polar statistical associating fluid theory, Fluid Phase Equilib. 472 (2018), 117–127.10.1016/j.fluid.2018.05.016Search in Google Scholar

[25] X. Liang and M. L. Michelsen, General approach for solving the density gradient theory in the interfacial tension calculations, Fluid Phase Equilib. 451 (2017), 79–90.10.1016/j.fluid.2017.07.021Search in Google Scholar

[26] Z. Qiao and S. Sun, Two-Phase Fluid Simulation Using a Diffuse Interface Model with Peng-Robinson Equation of State, SIAM J. Sci. Comput. 36 (2014), B708–B728.10.1137/130933745Search in Google Scholar

[27] X. Mu, F. Frank, F. O. Alpak and W. G. Chapman, Stabilized density gradient theory algorithm for modeling interfacial properties of pure and mixed systems, Fluid Phase Equilib. 435 (2017), 118–130.10.1016/j.fluid.2016.11.024Search in Google Scholar

[28] D. -Y. Peng and D. B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1976), 59–64.10.1021/i160057a011Search in Google Scholar

[29] R. Stryjek and J. H. Vera, PRSV: An improved Peng-Robinson equation of state for pure compounds and mixtures, Can. J. Chem. Eng. 64 (1986), 323–333.10.1002/cjce.5450640224Search in Google Scholar

[30] R. Stryjek and J. H. Vera, PRSV2: A cubic equation of state for accurate vapour-liquid equilibria calculations, Can. J. Chem. Eng. 64 (1986), 820–826.10.1002/cjce.5450640516Search in Google Scholar

[31] I. H. Bell and A. Jäger, Helmholtz energy transformations of common cubic equations of state for use with pure fluids and mixtures, J. Res. NIST 121 (2016), 1712–1738.10.6028/jres.121.011Search in Google Scholar PubMed PubMed Central

[32] J. Kou and S. Sun, Thermodynamically consistent modeling and simulation of multi-component two-phase flow with partial miscibility, Comput. Methods Appl. Math. 331 (2018), 623–649.10.1016/j.cma.2017.11.023Search in Google Scholar

[33] P. M. W. Cornelisse, The Square Gradient Theory Applied – Simultaneous Modelling of Interfacial Tension and Phase Behaviour, Ph. D. Thesis, TU Delft, Netherlands, 1997.Search in Google Scholar

[34] C. Miqueu, B. Mendiboure, A. Graciaa and J. Lachaise, Modeling of the surface tension of multicomponent mixtures with the gradient theory of fluid interfaces, Ind. Eng. Chem. Res. 44 (2005), 3321–3329.10.1021/ie049086lSearch in Google Scholar

[35] S. Stephan, K. Langenbach and H. Hasse, Interfacial properties of binary Lennard-Jones mixtures by molecular simulation and density gradient theory, J. Chem. Phys. 150 (2019), 174704.10.1063/1.5093603Search in Google Scholar

[36] S. Stephan and H. Hasse, Interfacial properties of binary mixtures of simple fluidsand their relation to the phase diagram, Phys. Chem. Chem. Phys. 22 (2020), 12544–12564.10.1039/D0CP01411GSearch in Google Scholar

[37] R. Krishna and J. A. Wesselingh, The Maxwell-Stefan approach to mass transfer, Chem. Eng. Sci. 52 (1997), 861–911.10.1016/S0009-2509(96)00458-7Search in Google Scholar

[38] J. Kim, S. Lee, Y. Choi, S. -M. Lee and D. Jeong, Basic principles and practical applications of the Cahn–Hilliard equation, Math. Probl. Eng. 2016 (2016), 1–11.10.1155/2016/9532608Search in Google Scholar

[39] M. Pavelka, V. Klika and M. Grmela, Multiscale Thermo-Dynamics, de Gruyter, 2018.10.1515/9783110350951Search in Google Scholar

[40] L. Landau and E. Lifshitz, Fluid Mechanics, vol. 6, Course of Theoretical Physics, Butterworth-Heinemann, 1987.Search in Google Scholar

[41] B. Haspot and E. Zatorska, From the highly compressible Navier-Stokes equations to the Porous Medium equation – rate of convergence, preprint (2021), https://arxiv.org/abs/1504.04219v2.Search in Google Scholar

[42] M. Caggio and D. Donatelli, High Mach number limit for Korteweg fluids withdensity dependent viscosity, J. Differ. Equ. 277 (2021), 1–37.10.1016/j.jde.2020.12.017Search in Google Scholar

[43] E. B. Nauman and D. Q. He, Nonlinear diffusion and phase separation, Chem. Eng. Sci. 56 (2001), 1999–2018.10.1016/S0009-2509(01)00005-7Search in Google Scholar

[44] S. I. Mekhtiev, A. A. Mamedov, S. K. Khalilov and M. A. Aleskerov, Experimentelle Untersuchung des Einflusses von Octylmethacrylat auf Viskosität und Dicht der Kohlenwasserstoffe, Izv. Vyssh. Uchebn. Zaved. Neft Gaz (1975), 64–100.Search in Google Scholar

[45] C. C. Chappelow, P. S. Snyder and J. Winnick, Density of liquid n-octane, J. Chem. Eng. Data 16 (1971), 440–442.10.1021/je60051a036Search in Google Scholar

[46] E. W. Lemmon, Thermophysical Properties of Fluid Systems, NIST chemistry WebBook, NIST standard reference database number 69, http://webbook.nist.gov, National Institute of Standards and Technology, 2005.Search in Google Scholar

[47] I. A. Tugarev, Z. I. Avdus and V. F. Nozdrev, An Experimental Study of the P-V-T Parameters of the Liquid Phase in Acetone-Benzene Mixtures Along the Boundary Curve, Zh. Fiz. Khim. 49 (1975), 1256–1258.Search in Google Scholar

[48] V. A. Borzunov, V. N. Razumikhin and V. A. Stekolnikov, Bestimmung der Dichte von n-Hexan und Wasser bei Drücken bis 10000 kg/cm2, Teplofiz. Svoistva Vesh. Mater (1970), 146–152.Search in Google Scholar

[49] I. F. Golubev, T. N. Vasilkovskaya and V. S. Zolin, Experimentelle Untersuchung der Dichte von alyphatischen Alkoholen bei verschiedenen Temperaturen und Drücken, Inzh. Fiz. Zh. 38 (1980), 668–670.Search in Google Scholar

[50] W. G. Chapman, K. E. Gubbins, G. Jackson and M. Radosz, SAFT: Equation-of-State Solution Model for Associating Fluids, Fluid Phase Equilib. 52 (1989), 31–38.10.1016/0378-3812(89)80308-5Search in Google Scholar

[51] C. Klink and J. Gross, A density functional theory for vapour-liquid interfaces of mixtures using the perturbed-chain polar statistical associating fluid theory equation of state, Ind. Eng. Chem. Res. 53 (2014), 6169–6178.10.1021/ie4029895Search in Google Scholar

[52] B. A. Grigorev, B. V. Nemzer and G. D. Tatevosov, Experimentelle Untersuchung der Oberflächespannung von n-Pentan, n-Hexan und n-Heptan, Izv. Vyssh. Uchebn. Zaved. Neft Gaz (1985), 53–58.Search in Google Scholar

[53] A. B. Medvedev, Estimating the self-diffusion and mutual diffusion coefficients of binary mixtures on the basis of a modified van der Waals model, Combust. Explo. Schock+ 53 (2017), 420–432.10.1134/S0010508217040062Search in Google Scholar

[54] M. M. Telo da Gama and R. Evans, The structure and surface tension of the liquid-vapour interface near upper the ciritical end point of a binary mixture of Lennard-Jones fluids, Mol. Phys. 48 (1983), 229–250.10.1080/00268978300100181Search in Google Scholar

[55] A. H. Falls, L. E. Scriven and H. T. Davis, Adsorption, structure, and stress in binary interfaces, J. Chem. Phys. 78 (1983), 7300–7317.10.1063/1.444720Search in Google Scholar

[56] D. J. Lee, M. M. Telo da Gama and K. E. Gubbins, Adsorption and Surface Tension Reduction at the Vapor-Liquid Interface, J. Phys. Chem.-US 89 (1985), 1514–1519.10.1021/j100254a041Search in Google Scholar

[57] S. Stephan and H. Hasse, Enrichment at vapur-liquid interfaces of mixtures: establishing a link between nanoscopic and macroscopic properties, Int. Rev. Phys. Chem. 39 (2020), 319–349.10.1080/0144235X.2020.1777705Search in Google Scholar

Received: 2021-10-26

Revised: 2021-12-26

Accepted: 2022-01-21

Published Online: 2022-02-04

Published in Print: 2022-04-30

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/jnet-2021-0080

Keywords for this article

density gradient theory; Peng–Robinson equation of state; gas–liquid interface; surface tension

Creative Commons

BY 4.0