Mathematical modeling of planar and spherical vapor–liquid phase interfaces for multicomponent ﬂuids

. Development of methods for accurate modeling of phase interfaces is important for understanding various natural processes and for applications in technology such as power production and carbon dioxide separation and storage. In particular, prediction of the course of the non-equilibrium phase transition processes requires knowledge of the properties of the strongly curved phase interfaces of microscopic droplets. In our work, we focus on the spherical vapor–liquid phase interfaces for binary mixtures. We developed a robust computational method to determine the density and concentration proﬁles. The fundamentals of our approach lie in the Cahn-Hilliard gradient theory, allowing to transcribe the functional formulation into a system of ordinary Euler-Langrange equations. This system is then split and modiﬁed into a shape suitable for iterative computation. For this task, we combine the Newton-Raphson and the shooting methods providing a good convergence speed. For the thermodynamic properties, the PC–SAFT equation of state is used. We determine the density and concentration proﬁles for spherical phase interfaces at various saturation factors for the binary mixture of CO 2 and C 9 H 20 . The computed concentration proﬁles allow to the determine the work of formation and other characteristics of the microscopic droplets.


Introduction
The phase interface is the layer separating bulk thermodynamic phases, the liquid and the gas phases in the present case. The thickness of the phase interface is typically on the order of nanometres. Yet, understanding the physics of vapor-liquid phase interfaces is important for numerous natural processes and applications in technology. In particular, the present study aims at modeling the phase transition phenomena in systems relevant for carbon capture and storage technologies. The subject of interest of this study are microscopic droplets, whose diameter is comparable with the thickness of the phase interface. Such droplets occur in the process of nucleation of a liquid phase from the metastable supersaturated gas phase.
The phase equilibrium state of the system has to be investigated first. For this purpose, an equation of state (EoS) is needed. The usual approach is to choose one of the simple low parameter equation (Van der Waals, Redlich-Kwong, Peng-Robinson or their modifications).
Here we employed an advanced PC-SAFT EoS. This equation belongs to a family of EoSs derived from the statistical fluid association theory (SAFT) [1]. The perturbed chain modification PC-SAFT was developed by Gross and Sadowski [2]. Various approaches are available to characterize the phase interface. The Gibbsian approach, treating the phase interface as infinitesimally a e-mail: : celnydav@fjfi.cvut.cz thin, allows to properly describe the integral properties of the phase interface such as the surface tension, molecular adsorption and related thermodynamic properties. However, it does not provide a tool to predict the dependence of the surface tension on composition of the adjacent phases and on the curvature of the phase interface. The so-called capillarity approximation, which is one of the main assumptions of the classical nucleation theory (CNT) [3], assigns properties of the macroscopic planar phase interface to the strongly curved interface of the microscopic droplets. Modeling the effect of the surface curvature is possible with the gradient theory(GT) [4,5] or density functional theory (DFT) [6]. These theories enable to obtain the density profiles (dependence of the local density on the coordinate perpendicular to the surface) and other characteristics of the phase interface.
This study is a continuation of the research of phase interfaces using GT previously done by Hrubý et al. [7,8] and Planková [9]. In this work, we propose a new approach to the mathematical solution of the problem obtained from GT. This solution is iterative and offers a way of splitting the problem into two consecutive steps that can be solved separately.

The modeled system
We consider an n-component mixture (two components in this study) in the state of metastable supersaturated gas phase. We consider a microscopic droplet in the thermal and chemical equilibrium with the gas phase, meaning that the temperature and chemical potentials are homogeneous throughout the system. The modeled droplet corresponds to the so-called critical cluster of the molecular theory. We note that the opposite situation with a bubble in the equilibrium with a superheated or stretched liquid would be very much alike.
When the system is composed of two components, the phase equilibrium state with a planar phase interface can be easily computed for the given temperature T and pressure p based on the conditions of thermodynamic equilibrium of the liquid and gaseous phases: equality of temperatures, chemical potentials for both components, and equality of pressures in both phases. As a result, we obtain the densities and compositions of both phases.
To proceed from the system with a planar phase interface to the system with a droplet, we increase the concentration of the less volatile component (here n-nonane) in the mother phase (gas) by multiplying it by the saturation ratio α > 1, while keeping the concentration of the more volatile component (carbon dioxide) and temperature constant. The saturation ratio directly influences the size of the droplet: a bigger saturation ratio results in a smaller droplet. We note that choosing a too high saturation ratio might result into an unstable state (a state beyond the vapor-liquid spinodal). For the supersaturated gas phase obtained in this way, we determine the liquid phase in thermal and chemical equilibrium. From the condition of equality of chemical potentials of each components in both phases we determine the concentrations of the components in the liquid phase and the corresponding pressure of the liquid phase, which is bigger than for the mother gas phase. The magnitude of difference is also directly influenced by the saturation factor in a manner that higher α means higher pressure difference. The thermodynamic state of the liquid obtained in this way is the state of idealized liquid droplet in the Gibbsian picture. The difference of pressure is related to the Laplace pressure due to the surface tension of the curved interface. The "equilibration" of the droplet with the gas phase is iterative process. The iteration is continued even a few steps after the stop criterion is reached and the best result is chosen among the computed after-steps. The purpose of this measure is avoiding wrong results due to premature termination of the iteration.
The thermodynamic state obtained in this way is hypothetical because even the center of the droplet is influenced by the phase interface. However, the density and concentrations obtained in this way are used to provide the first guess of local conditions in the center of the droplet in the gradient theory modeling.

Gradient theory
The gradient theory, also called Cahn-Hilliard theory [4], was developed by Cahn and Hilliard in 1957, although its foundations were already given by van der Waals. The basis of the GT is an expansion of the Helmholtz energy into a Taylor series with respect to the gradient of density or gradients of concentrations. This expansion is then used to obtain the grand potential of the system, which is treated as a functional of the density or concentration profiles. These profiles are then determined by minimizing the functional, or more precisely, finding its stationary point by solving the corresponding Euler-Lagrange (E-L) differential equations.
In this study, we define the functional as the work of formation given by Eq. (1): Here, Ω is the grand potential. The subscripts in this equation indicate homogeneous system (only the gaseous mother phase) and inhomogeneous system (system including a droplet with some gaseous surroundings at the same temperature and chemical potentials). The grand potential Ω can be expressed in terms of the Helmholtz free energy of the inhomogeneous system. The form of the functional can be transformed into the spherical coordinates. The changed coordinate system allows using the symmetry of spherical droplet and simplifies the expression into the following form: (2) We can see from Eq. (2) that integral is one dimensional and all concentrations ρ i depend on the radial coordinate only. c i,k in Eq. (2) denotes the so-called influence parameter. And the term Δω is a composite of the Helmholtz energy density and term i ρ i μ i . Also this term depends on density profile function ρ(r). We denote density profile function in a vector shape because it is composed of concentration profiles for both components.
For given density profile the functional (2) returns the work needed to form an inhomogeneity represented by this density profile. Our task is now to obtain the most probable density profile for the given system conditions. Such demand according to the paper [4] requests a minimization of the formation work. The solution yields saddle point and minimal points of functional W. We are interested only in the saddle point because it corresponds to the required density profile for the droplet. The minimum points on the other hand represent a homogeneous system (gaseous or liquid) without any phase interface.
The saddle point is found using E-L system of two second order differential equations Eq. (3): Where ρ denotes differential of dρ/dr. Expressing the work of formation in Eq. (3) in terms of the GT, we obtain Eq. (4): Here, the right hand side of Eq. (4) represents differences between chemical potentials of a hypothetical homogeneous phase corresponding to local concentrations and the chemical potentials of the mother phase. In case that the interaction of the unlike molecules is not qualitatively different from the interactions of the like molecules, it is possible to approximate the influence parameter of mutual interaction of components c i,k as a geometric mean of the influence parameters for the pure components Eq. (5) is the basis for the present computations. Various methods can be developed to solve this equation. The present method provides benefits of substantial speed-up, modifiability of computation and result precision.

Solution of the Euler-Lagrange equations
Our method uses elementary results of the gradient theory and modifies it into a two step iterative method. We notice the similarity in right hand side of equation within Eq. system (4) and subtract the first equation from the others. This operation gives us system composed of one differential equation (6) and n − 1 non-linear algebraic equations (7).
This system contains equations coupled through their right-hand sides. Our goal is to split the solution into solving algebraic part (7) and the differential part (6) separately. To ensure that these equations are coupled, the algebraic part (7) has to be solved first. Then, its results are used in the computation of equation (6). We perform a special substitution:ρ The substitution transforms the components concentrations ρ i into modified molar densityρ. Substitution also ensures simple shape of E-L equations and enable discretization. We can also see that shape of substitution (8) is invariable to differentiation d/dr. This feature ofρ is used in differential equation. Secondly, a new artificial variable X is inserted.
To clarify the dependence in system of Eq. (7, 6) we choose . . n we then obtain a system: This system (9,10) is basis of following computation. The solution of equation (9) is possible only with the knowledge of artificial variable X in the molar density interval. For that purpose X has to be computed. The dependence of X variable enables the discretization ofρ across the molar density interval. At these discretized points we compute values of X.

Solution of the algebraic equations
As the first step we take only algebraic part (10). The artificial variable X allows to rewrite the algebraic part into a system of similar algebraic equations. As a result of adding new variable, one new equation also needs to be appended. The last equation is provided by definition (8). This completes the analytic system of equations (11). , where for j ∈ 2 . . . n we have: The system (11) needs to be solved for values of X and ρ at each point of discretization. This provides the function dependence sampled across the molar density interval. Results of each iteration X(ρ i ), ρ(ρ i ) are saved into data matrix M 1 . The secondary computational data (dX(ρ i )/dρ, dρ(ρ i )/dρ) are also saved into data matrix M 2 . The data matrices are useful in post-processing phase for computation of approximate of requested functions X and ρ. To solve this system the Newton-Raphson iterative method is used. For sufficiently fine division (500 points or more for best performance) the iteration is commenced. The Newton-Raphson method requires to compute the Jacobian of system (11) and sufficient first estimate. First guess for boundary points is provided by boundary conditions of system that we obtained from supersaturated equilibration. And points inside discretization use previous result as their first estimate. Owing to the fine discretization, the solution for each point is found in two to five Newton-Raphson iterations.
When the computation is completed, the obtained matrices M 1 and M 2 are used to compute piecewise cubic approximation of functions X(ρ), ρ 1 (ρ), and ρ 2 (ρ). Examples of these functions for saturation ratio α = 1.5 are shown in figures 1 and 2. For each point of discretization except the first we compute coefficients for cubic polynomial with function values from M 1 and with derivatives from the M 2 . The three piecewise cubic curves then provide a link between the two parts of the problem and are used for the solution of differential equation.

Solution of the differential equation
The differential part is the first equation (9). It is a second order differential equation which is now decoupled from the algebraic solution of the rest of the system. As input we use the piecewise cubic function X(ρ) computed from the solution of the algebraic system. The differential equation is accompanied with the boundary condition from the saturated system. From the mother phase end, we know the modified density of the gas phaseρ G and its derivative equal to zero. In the center of the droplet we know only that Fig. 1. Function X(ρ) computed by solving system of equations (11) for saturation ratio 1.5. / Fig. 2. Dependencies of the concentrations ρ 1 and ρ 2 on the modified densityρ computed by solving system of equations (11) for saturation ratio 1.5 the derivative of the modified density is also zero as a consequence of the spherical symmetry.
The solution method is based on the decomposition of the second order equation into system of two first order equations. The boundary problem is then solved by a shooting method, converting the boundary problem into the initial one. Our shooting parameter is the initial modified liquid densityρ 0 . The solution of the initial problem is computed with solving method ODE45, which is Runge-Kutta (4,5) method within MATLAB (R) libraries. As an option we requested to stop the computation either when the computed modified density reaches a negative value or its derivative equals zero.
The next value of shooting parameter is computed using the bisection method. At first the profile is computed for initial shooting parameter. The minimal valueρ min of the profile is found and applied for the decision of the next shooting parameterρ 0 .
We use the robust decision criteria that takes only minimal value of the profileρ min into account. If theρ min value is below theρ G level then lower shooting parameterρ 0 is needed. Similarly, whenρ min is above the gas density levelρ G , the higher shooting parameterρ 0 is chosen.
The iteration continues until specified precision is reached, or size of the interval is below numeric precision of double number format. When the iteration terminates successfully, the profile found this way is denoted as result profile and shooting parameter asρ 0,res .

Results
We implemented the splitting method and computed the density profiles for the system of carbon dioxide CO 2 and n-nonane C 9 H 20 . The phase equilibrium system conditions were set to T = 350 K and p = 10 MPa. The saturation factor alpha for this system was determined to range from α ∈ 1.15, 4 .
The discretization constant for the algebraic solution was set to 1000 points inside interval ρ G ,ρ L . This option provided enough division points during the iteration and accurately use the interpolation. Performed analysis shows that the discretization interval can be chosen even ten times smaller, still providing accurate results. This fact is illustrated in table 1 and is caused by the efficiency of the piecewise cubic interpolation.
As problematic appear areas outside of the modified density interval. Here, extrapolation is used which is sensitive to errors in interpolation coefficients. However, extrapolation is only used for initial iteration or for unphysical profiles which might occur during the iteration. The result profile does not suffer from the extrapolated data.
Table (1) illustrates also the iteration dependency. Rough discretization have the average ratio of iteration per point around 3:1 but as the refinement increases the iteration ratio changes to 2:1. This nicely correspond to section section Algebraic solution where this behaviour was explained. The last column of table represent precision of result profile initial condition computation. The initial condition is our shooting parameter and we expect it to improve as the refinement of discretization improves. The space in values is inserted for clarity reasons at position of first different digit. The problem of precision is illustrated by following figure 3 that represents the behaviour of profiles for different initial conditions (shooting parameter)ρ 0 . The figure depicts five initial conditions. But only two (red The search for optimal starting condition exhibits interesting behavior when approaching the optimal value. This behavior is captured in figure 4. Compared to previous figure, we plot only the dependency of found minimum of the density profile on initial modified density. The black full line from figure 4 is composed of two segments. A linear segment and a "logarithmic" segment. The linear segment is of no importance because it only illustrate the incrementation of the initial conditionρ 0 . The interesting situation happens for the logarithmic segment. This segment represents the search for optimal value of initial condition that is further used as result.
We know that the optimal profile is the one with minimum value equal to the modified gas densityρ G and, because of that the optimal profile is found at the intersection with the red dotted line on figure 4 representing ρ G . Before this value we find profiles similar to cyan square line from previous figure 3 that snap to thermodynamically unstable solution. And after magenta circle we see rapid descend into negative region of figure that corresponds to the green cross line from previous figure 3. This is the reason for decision criteria to evaluate the minimal values of profiles. Because of this feature we are able to distinguish the small details, where the profile changes from cyan square line into green cross line. This is the moment we are looking for and even small variation in the initial condition influence the profile dramatically.
The common problem of such behaviour is the precision of double number format. Available range of representable numbers is not sufficient for capturing optimal profile. The profile is suboptimal and display only small region where it satisfy the right boundary conditionρ opt =ρ G . The use of other number format would not solve the problem only extend the range of computable profiles. We have not performed the reimplementation into more precise number format because double is the greatest number type natively supported by MATLAB environment. The variable precision that is supported by minority of MATLAB function is also unacceptable due to is poor performance which would cause overall slowdown of computation.
For the reasons given above we chose the initial conditions in the safe region and perform computation there. We computed specified profile for three saturation factors of α ∈ {4, 1.5, 1.15} respectively to illustrate the growth of the droplet with decreasing saturation factor. With the modified densityρ used in our computation and comparison of results we also computed components densitiesρ 1 andρ 2 that correspond to CO 2 and C 9 H 20 respectively. The boundary conditions were also included. In our first system depicted in figure 5 we notice that only small droplet is formed. Due to the high saturation factor, the gas density is noticeably higher compared to figures 6 and 7. This figure also illustrates the droplet with no bulk liquid content. The droplet is so small in diameter that only phase interface is present. This fact is denoted by the density profile (black full line) that does not reach the blue dashed lineρ L .
The second system in figure 6 captures a small droplet with nearly no liquid content. This figure corresponds to the saturation factor of α = 1.5 and show interesting behavior of carbon dioxide. Carbon dioxide adsorbs at phase interface. It means that its molecules accumulate at the phase interface and density profile than have small peak within the phase interface. This behavior is noticeable also in figure 7.
The last figure shows an example of a large droplet. A significant portion of droplet is bulk liquid. If saturation factor would be lowered further, the droplet would grow until the planar phase interface is formed. This situation would happen at saturation factor α = 1 but computationally it happens sooner due the finite machine 02011-p.5 precision of computation. This error is most significant for bigger droplets with high ratio of pure liquid compared to phase interface. This problem is currently being solved with use of an analytical solution near zero radius.

Conclusion
Of the various approaches to the phase interface modeling, the gradient theory provides a good compromise between the physical detail and mathematical tractability. For this theory we developed a new solution method for multi-component mixtures and tested it for a binary system. In this work we presented the main steps of the method and offered information on variable parameters. We have also highlighted the problematic areas of computation which primarily concern very large droplets with almost "top hat" density profile.