Mass-conserved density gradient theory model for nucleation process

Subscriber access provided by Kaohsiung Medical University

Thermodynamics, Transport, and Fluid Mechanics

Mass-conserved density gradient theory model for nucleation process Xiaoqun Mu, Florian Frank, Beatrice Riviere, Faruk Omer Alpak, and Walter G Chapman Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b03389 • Publication Date (Web): 24 Oct 2018 Downloaded from http://pubs.acs.org on October 28, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Mass-conserved density gradient theory model for nucleation process† Xiaoqun Mu,‡ Florian Frank,¶ Beatrice Riviere,§ Faruk O. Alpak ,k and Walter G. Chapman∗,‡ ‡Department of Chemical and Biomolecular Engineering, Rice University, USA ¶Department of Mathematics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Germany §Department of Computational and Applied Mathematics, Rice University, USA kShell International Exploration and Production Inc., USA E-mail: [email protected]

Abstract Modeling droplet nucleation processes requires a molecular-scale approach to describe the interfacial tension (IFT) of spherical interfaces. Density gradient theory (DGT), also referred to as square gradient theory in some publications, has been widely used to compute the IFT of many pure and mixed systems at the molecular scale. However, the application of DGT to droplet interfaces is limited by its setup in open systems, in which a stable droplet can not be achieved. In this article, we propose a mass-conserved DGT model in a closed system, i. e. with no-flux boundary conditions, where no mass exchange is allowed with the outside environment. As opposed to the traditional approach, this model enforces a canonical ensemble and guarantees energy dissipation. The proposed model has been successfully applied to systems with †

Submitted to Industrial & Engineering Chemistry Research

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32

planar interfaces as well as spherical interfaces, especially for droplet’s IFT calculation in the nucleation process. By extending the DGT model from open to closed systems, we demonstrate the potential of DGT as an inhomogeneous model for a wider range of academic and industrial applications.

Introduction The nucleation process is a commonly observed phenomenon that occurs in nature, science and industry, e. g., crystallization of salt from sea water, vapor condensation in the condenser of a distillation tower, or simply water freezing into ice. In this work, we are interested in homogeneous droplet nucleation, in which a liquid droplet forms from the supersaturated vapor phase away from a solid surface. Studying the homogeneous droplet nucleation process is of significant importance. For example, in the artificial precipitation operation, anthropogenic aerosol particles serve as seeds (nuclei) for cloud droplets to nucleate. One of the biggest challenges in simulations of climate and climate change is modeling the droplet nucleation process in clouds, which represent much of the complexity of the simulation 1,2 . In the nucleation process, the interfacial tension (IFT) is a key factor that affects the droplet’s properties such as the radius and the nucleation rate 3 . More than 200 years ago, Young 4 and Laplace 5 introduced the well-known Young–Laplace equation which correlates IFT of curved interfaces to the principal radii of curvature and the pressure difference across the interface, ∆P . For a spherical droplet, the Young–Laplace equation reduces to:

∆P =

2γ , R

(1)

where γ is the IFT on the droplet surface and R is the droplet’s radius. The Young–Laplace equation can be derived by purely thermodynamic arguments 6 , and it is valid for spatial scales that are large enough to support the sharp interface assumption. Thermodynamics approaches also approximate the effect of droplet size on vapor pressure ac2


Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cording to the Kelvin equation 7 , and the effect of droplet radius on the IFT from Tolman’s equation 8 . Unfortunately, thermodynamic approaches do not take molecular interactions into consideration, and thus lose accuracy in IFT calculations for small droplets, especially when the system’s interface thickness is comparable to the droplet radius. This is due to the fact that the system’s inhomogeneity is amplified as the droplet radius decreases, and the regional differences of molecular interactions have greater impacts on fluid properties. Therefore, models that can provide a direct link between the microscopic molecular interactions, and macroscopic thermodynamic properties are needed to compute the IFT of nano-scale droplets. Molecular dynamics simulation, which includes the forces experienced by each molecule, has been successfully implemented to model nucleation processes 3,9–12 at the expense of long computational time. Classical density functional theory (DFT), which is derived from statistical mechanics, also serves as a powerful tool in describing thermodynamic properties of inhomogeneous systems, and many studies have applied classical DFT to spherical systems 13–16 . Density gradient theory (DGT) is special form of DFT, which uses an equation of state (EoS) to describe the homogeneous part of the fluids, and a density gradient term to account for the inhomogeneity of the system. Its origin dates back to Van der Waals 17 , and the modern version became popular after Cahn and Hilliard’s reformulation of this theory 18 . With a history of over 100 years, DGT is still gaining much attention from both industry and academia, thanks to its simplicity in implementation and the accuracy in predicting interfacial properties. The conventional DGT model is based on the grand canonical ensemble system, in which chemical potential µ, volume V , and temperature T are fixed, while the system has a open boundary that allows mass exchange. There have been many successful applications of DGT in systems with planar interfaces, including alkane mixtures forming vapor–liquid equilibrium (VLE) systems 19 and water–oil mixtures forming liquid–liquid equilibrium (LLE)

3



systems 20,21 etc. In our previous work, we developed a novel stabilized DGT (SDGT) algorithm 22,23 , which allows robust and efficient solution of the DGT type of equations. In 24 , we developed a modified DGT model for surfactant molecules with more sophisticated chain structures. In the literature review, we only found a few publications attempting to apply the DGT model to spherical interfaces. However, these authors either did not conduct any IFT calculations 25 or did the IFT calculation by combining DGT with Young–Laplace equation 26 sacrificing DGT’s advantages at the nano-scale. In fact, the SDGT algorithm with current formulation of the DGT model in a grand canonical ensemble is incapable of generating a stable droplet with finite radius. This is due to the fact that the grand potential energy in an open system has a saddle point that corresponds to the droplet with the so-called critical nucleation radius Rc . When the droplet radius is smaller than Rc , the droplet keeps evaporating until no liquid phase is left; when the droplet radius is larger than Rc , it will grow until the entire system becomes liquid phase. In a recent work by Langenbach et al. 27 , conventional DGT model has been successfully applied to nucleation calculation. In our calculations however, we found it very challenging to capture the critical nucleation radius, since any small perturbation on the metastable system will lead to droplet growth or shrinkage. A closed system with no mass exchange, on the other hand, supports the stability of droplets. One example of a closed system is the canonical ensemble system, in which number of particles N , volume V , and temperature T are fixed. In this work, we develop a mass-conserved DGT model that works in canonical ensemble systems by introducing a massconservation term λi and subjecting the model equations to no-flux boundary conditions. A modified SDGT algorithm 22 is used to solve the mass-conserved DGT equations, from which the equilibrium density profiles are obtained. The mass-conserved DGT model has been applied in IFT calculations of both planar interface and spherical interface systems, the results of which are compared with those of the conventional DGT model and other theories.

4


Page 4 of 32

Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The mass-conserved DGT model DGT equations in closed systems For an inhomogeneous fluid consisting of N components occupying a volume V at fixed temperature T , the Helmholtz free energy A is approximated using the following expression: Z A(ρ) = A0 (ρ) + A∇ (ρ) =

a0 ρ(r) dr +

V

Z V

N 1 X vij ∇ρi (r) · ∇ρj (r)dr , 2 i,j=1

(2)

for i = 1, . . . , N with A0 being the homogeneous Helmholz free energy, which is given by a bulk EoS. The vector of molar densities ρ indicates the dependency of all molar densities ρi . The gradient energy A∇ includes the (constant) influence parameters vij . In a canonical ensemble system, the Helmholtz free energy A is the characteristic energy that is minimized at equilibrium. The minimization of A leads to N

X δA = µ0i (ρ) − vij ∇2 ρj , δρi j=1

(3)

where δA/δρi is the functional derivative of the Helmholtz free energy A with respect to molar density ρi . The term µ0i (ρ) := ∂a0 (ρ)/∂ρi is the homogeneous chemical potential of component i given by a bulk EoS. In this work, the perturbed chain statistical association fluid theory (PC-SAFT) EoS 28,29 is employed for µ0i (ρ). With fixed number of particles in the system, the following constraint holds for each component i: Z ρi (r) dr = Ni ,

(4)

V

where Ni is the total molar amount of component i in the system. Since the volume V is assumed closed, there is no mass flux across the boundary ∂V , i. e., the homogeneous

5



Page 6 of 32

Neumann boundary condition holds:

∇ρi · n = 0

(5)

on the boundary ∂V with n denoting the unit normal outward of the domain V . In previous works of classical DFT calculations 15,16 , the density profile ρi = ρi (r) with mass-conservation is computed by

µex i (r) RT

Ni exp ex 0 , µi (r ) 0 exp dr RT V

ρi (r) = R

(6)

where µex i is the excess chemical potential. Eqn. 6 implies mass conservation (4), but it does not enforce a closed boundary of the system, cf. (5). In other words, the system is still a grand canonical ensemble with an open boundary, while the mass entering the system equals the mass existing. In the mass-conserved DGT model, we are interested in finding density distributions ρi in V that minimize the free energy A, cf. (2) while enforcing the mass conservation (4) and no-flux boundary conditions (5). In order to embed these conditions into an energyminimization process, we introduce the pseudo-time variable t and the mass-conservation term λi (t). Subjecting the Euler–Lagrange equation to a time-evolution term ∂ρi /∂t 22 and λi (t), the density profiles of interest are hence those at the stationary state of equation 1 ∂ρi (t, r) δA(t) + + λi (t) = 0 for i = 1, 2, . . . , N . M ∂t δρi

(7)

In Eqn. (7), the quantity M > 0 is a phenomenological mobility coefficient, which is assumed constant and equal for every component. The mass-conservation term λi = λi (t) is assumed to depend on time only (not on space) and is motivated by Rubinstein and Sternberg 30 who used a similar expression in the context of the Allen–Cahn equation 31 .

6


Page 7 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Mass-conservation term λi In Eqn. 7, the term λi is supposed to enforce mass conservation of the system. The closedform expression of λi can be derived analytically: With the mass conservation assumption (4), the total mass of each component i does not change in time, i. e. d dNi = dt dt

Z

Z ρi (t, r) dr =

V

V

∂ρi (t, r) =0. ∂t

(8)

Substituting Eqn. 7 into Eqn. 8 yields: Z M V

δA(t, r) + λi (t) dr = 0 δρi

(9)

for all time points t. Since λi is only a function of time, it follows Z

1 λi (t) = − |V |

V

δA(t, r) dr , δρi

(10)

where |V | is the measure of the volume V . Eqn. 10 is the general expression of λi for any formulations of inhomogeneous Helmholtz free energy. Within the DGT framework, this integral can be expanded and simplified: 1 λi (t) = − |V |

Z

1 |V |

Z

(3)

=−

N X 0 µi ρ(t, r) dr − vij

V

Z

µ0i ρ(t, r) dr −

V

N X j=1

∇2 ρj (t, r) dr

V

j=1

!

!

Z

∇ρj (t, r) · n dr ,

vij

(11)

∂V

where Gauss’ theorem is used to transform the volume integral to the surface integral. Due to the boundary condition (5), the surface integral in Eqn. 11 vanishes, hence: 1 λi (t) = − |V |

Z

µ0i ρ(t, r) dr .

V

7


(12)


Page 8 of 32

Using λi (t) in Eqn. 7 leads to the mass-conserved DGT model (18a) for canonical ensemble systems.

Physical justification of λi To justify the mass-conservation term λi physically, we start from the concept of diffusion due to the differences in chemical potentials. Consider an inhomogeneous system at time t with a density distribution of ρ(t, r) in the domain as shown in Figure 1, where r is the position vector. In the diffusion process, a reference cell at position r0 can only have mass exchange with its neighboring cells, i. e., an adjacent cell at some position r1 . The density of component i at a fixed time point t + ∆t in the cell at r0 equals the current density plus the amount of mass diffused into the cell during a small time interval ∆t: h i ρi (t + ∆t, r0 ) = ρi (t, r0 ) + ∆tM µi (ρi (t, r1 )) − µi (ρi (t, r0 )) .

(13)

Outside environment Target system

r

1

r

0

Figure 1: In the diffusion process, mass is exchanged between adjacent cells, i. e., at r0 and r1 . Eqn. 13 represents the change of density in time due to the driving force of chemical potential differences. In the DGT model, we focus on the equilibrium state, not on the process that the system takes to reach equilibrium. Therefore, we relax the constraint that mass exchange only occurs with neighbors, and let each cell have mass exchange with the 8


Page 9 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


outside reservoir directly, which has a constant chemical potential µi,bulk . Eqn. 13, which describes diffusion, then becomes: h i ρi (t + ∆t, r0 ) = ρi (t, r0 ) + ∆tM µi,bulk − µi (ρi (t, r0 )) .

(14)

By letting ∆t → 0, Eqn. 14 becomes 1 ∂ρi (t, r0 ) + µi (ρi (t, r0 )) − µi,bulk = 0 , M ∂t

(15)

which is exactly the familiar objective function of the stabilized DGT equation in grand canonical ensemble systems 22 . Now consider a closed system, i. e., mass exchange with the outside environment is not possible. Instead, we let each cell have mass exchange with any other cells inside the system regardless of distance to achieve an overall equilibrium. This can also be interpreted as a process of redistributing molecules inside the closed system as long as chemical potential differences exist. In this scenario, µi,bulk in Eqn. 15 is replaced by the system’s average chemical potential at each time point to describe the new diffusion rule: 1 1 ∂ρi (t, r0 ) + µi (ρi (t, r0 )) − M ∂t |V |

Z

µi ρi (t, r 0 ))dr 0 = 0 ,

(16)

V

in which the last term represents the mass conservation term λi (t). It is updated over time to reflect the current difference between the local chemical potential and the system average chemical potential, and equilibrium is achieved when no difference exists. The mass conservation is validated by integrating Eqn. 16 over V : 1 M

Z V

∂ρi (t, r0 ) dr0 = ∂t

Z " V

1 |V |

#

Z

0

0

µi ρi (t, r ))dr − µi ρi (t, r)) dr0 = 0 , V

i. e., there is no change of mass in time.

9


(17)


Page 10 of 32

Mathematical model Model equations The mass-conserved DGT model is represented by the following system of partial differential equations, whose stationary state solutions are the molar density distributions that satisfy the energy-minimum equation (7) and the mass-balance constraint (4) for given molar amounts Ni : For component i = 1, . . . , N , we seek density distributions ρi = ρi (t, r) as a function of time t ∈ J and space r ∈ V such that: N X 1 ∂ρi (t, r) 1 0 + µi ρ(t, r) − vij ∇2 ρj (t, r) − M ∂t |V | j=1

Z

µ0i ρ(t, r 0 ) dr 0 = 0 ,

(18a)

V

subjected to the boundary condition

∇ρi (t, r) · n = 0 for t ∈ J, r ∈ ∂V ,

(18b)

ρi (0, r) = ρ0i (r) for r ∈ V ,

(18c)

and initial condition

where the time interval J := (0, ts ) is assumed large enough such that the system reaches a stationary state at time ts , and the initial data ρ0i = ρ0i (r) is chosen to satisfy Eqn. 4. System (18) is an extension of the system used by the SDGT algorithm 22 , which introduced an artificial time dependence term in order to approach the stationary state gradually instead of attempting to solve the stationary state equations directly. In this work, we made several modifications to the SDGT algorithm. Firstly, system (18) is formulated for the NVT setting, and thus prescribes Neumann boundary conditions (18b) instead of having Dirichlet boundary conditions. Secondly, the time evolution term in Eqn. 18a is physically consistent by having the mobility coefficient M . Furthermore, −µi,bulk is replaced by the

10


Page 11 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


term λi (t) = − |V1 |

R V

µ0i ρ(t, r 0 ) dr 0 , which guarantees mass-conservation (4) of each com-

ponent over time.

Energy dissipation In this section, we prove that solutions of system (18) satisfy the following energy dissipation law: N Z 1 X ∂ρi (t, r) 2 dA(t) =− dr ≤ 0 for t ∈ J , dt M i=1 V ∂t

(19)

i. e, the Helmholtz free energy A dissipates over time in the NVT setting. Proof. We suppress the argument (t, r) of ρi for compactness. Multiplying Eqn. 18a by ∂ρi /∂t, integrating over V , and summing up over i yields:

−

N Z X i=1

V

N

X 1 ∂ρi 2 dr = M ∂t i=1

Z

∂ρi µ0i (ρ) ∂t

V

N X

dr −

i,j=1

Z

∂ρi dr ∂t V Z Z N X ∂ρi 1 µi ρ) − dr |V | V V ∂t i=1 ∇2 ρj

vij

(20)

=: Γ1 + Γ2 + Γ3 . Using the definition of µi , term Γ1 is equivalent to Z Z N Z X ∂a0 (ρ) ∂ρi ∂a0 (ρ) d dr = dr = a0 (ρ) dr . ∂ρi ∂t ∂t dt V V i=1 V For term Γ2 , we use the identity Z N d 1X vij ∇ρi · ∇ρj dr dt 2 i,j=1 V

! =

=

N X

Z ∇ρj · ∇

vij

i,j=1

V

N X

Z

(18b)

∂V

= −

N X

Z vij

i,j=1

11

i

∂t

dr

∂ρi dr − ∇ρj · n ∂t

vij

i,j=1

∂ρ

V

∇2 ρ j

∂ρi dr . ∂t


Z V

∂ρi ∇ ρj dr ∂t 2

(21)


Page 12 of 32

Term Γ3 vanishes due to mass conservation (8). Adding Γ1 , Γ2 , and Γ3 leads to dA(t)/dt according to Eqn. 2. Therefore, the energy dissipation law (19) holds.

Planar interface case Consider a planar interface system in a rectangular cuboid domain V that is aligned with Cartesian axes. Let the density vary along the z direction only i. e. it is homogeneous in x and in y direction. Assume that the domain size (0, D) in z direction is large enough for the interface boundaries to locate in the two bulk phases. Eqn. 18 then simplifies to: For i = 1, . . . , N , we seek density distributions ρi = ρi (t, z) as a function of time t ∈ J and space z ∈ (0, D) such that N X ∂ 2 ρj (t, z) 1 1 ∂ρi (t, z) 0 + µi ρ(t, z) − vij − 2 M ∂t ∂z D j=1

Z

D

µ0i ρ(t, z) dr = 0

(22a)

0

with the boundary conditions for each component i, ∂ρi (t, 0) = 0 and ∂z

∂ρi (t, D) =0 ∂z

(22b)

for t ∈ J, where the initial condition ρi (0, r) is chosen to satisfy the mass constraint Z As

D

ρi (0, z) dz = Ni ,

(22c)

0

where As is the x-y surface area.

Spherical interface case Without loss of generality, assume that V is a sphere with radius R having its origin at the center of the droplet. The radius R is assumed large enough such that ∂V is entirely in one bulk phase. In this setting, the density ρi is only a function of the radial distance r to the origin. The problem given in (18) then reads: 12


Page 13 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For i = 1, . . . , N , we seek density distributions ρi = ρi (t, r) as a function of time t ∈ J and space r ∈ (0, R) such that N X 1 ∂ρi (t, r) 0 + µi ρ(t, r) − vij M ∂t j=1

∂ 2 ρj (t, r) 2 ∂ρj (t, r) + ∂r2 r ∂r

!

3 − 3 R

Z

R

r2 µ0i ρ(t, r) dr = 0

0

(23a) with the boundary conditions for each component i, ∂ρi (t, 0) = 0 and ∂r

∂ρi (t, R) = 0. ∂r

(23b)

The initial conditions ρi (0, r) are chosen to satisfy (4), which simplifies to Z 4π

R

r2 ρi (0, r) dr = Ni .

(23c)

0

Discretization Space discretization The systems of time-dependent boundary value problems (22) and (23) are solved with a time-marching scheme until a stationary state is reached. On every time level, a finite difference scheme is used to discretize equations in space. By interpolating S points inbetween the boundaries, 0 =: z0 < z1 < . . . < zS+1 := D denote equidistant nodes on a grid with mesh size ∆r :=

D S+1

for the planar interface case. By suppressing the time

argument here and using the abbreviation ρi,k := ρi (zk ), We approximate the first and second derivatives by central difference quotients: dρi (zk ) ρi,k+1 − ρi,k−1 ≈ , dz 2∆z

d2 ρi (rk ) ρi,k−1 − 2ρi,k + ρi,k+1 ≈ , 2 dz ∆z 2

13


(24a)


Page 14 of 32

At the boundary, the first derivative dρi /dz vanishes due to the homogeneous Neumann boundary condition, and the second derivatives are approximated by d2 ρi (z0 ) 2ρi,1 − 2ρi,0 ≈ dz 2 ∆z 2

d2 ρi (zS+1 ) 2ρi,S − 2ρi,S+1 ≈ . dz 2 ∆z 2

and

(24b)

Same discretization is applied to the spherical interface case with equidistant nodes 0 =: r0 < r1 < . . . < rS+1 := R.

Time discretization Let ∆tn denote the time step size between time level tn+1 and tn , which may vary from step to step according to a time step adaptation scheme 22 . The time derivative is approximated by the backward difference quotient (suppressing the space index): ρn+1 − ρni ∂ρi (tn+1 ) ≈ i ∂t ∆tn

(25)

A convex–concave splitting 22,23 is applied to the nonlinear part of the equation, which is the expression of the homogeneous chemical potential:

µ0i ρ(tn+1 ) ≈ µ0,+ ρn+1 + µ0,− ρn , i i

(26)

where µ0,+ is associated with the convex part of the energy and µ0,− with the concave part. i i

Fully discretized equations Discretizing Eqn. 22a and Eqn. 23a in time and space, while applying the convex–concave splitting (26) yields the following fully discrete problem. For every time level n + 1, we seek ρn+1 i,k such that for k = 1, . . . , S, i = 1, . . . , N , the residual equations equal zero:

14


Page 15 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


N n+1 n+1 n X ρn+1 ρn+1 j,k−1 − 2ρj,k + ρj,k+1 i,k − ρi,k 0,− 0,+ n+1 n vij + µi ρk − + µi ρk M ∆tn ∆z 2 j=1 Z 1 D 0,+ n+1 n dz = 0 ρ + µ0,− − µi ρ k k i D 0

(27)

subjected to the boundary conditions n+1 ρn+1 j,0 − ρj,1 = 0 and ∆z

n+1 ρn+1 j,S − ρj,S+1 = 0. ∆z

for the planar interface case and N n+1 n+1 n X ρn+1 ρn+1 j,k−1 − 2ρj,k + ρj,k+1 i,k − ρi,k 0,+ 0,− n+1 n + µi (ρk ) + µi (ρk ) − vij M ∆tn ∆r2 j=1 ! Z R n+1 n+1 3 2 ρj,k+1 − ρj,k−1 0,− n+1 n − 3 r2 µ0,+ + (ρ ) + µ (ρ ) dr = 0 , k i i k rk 2∆r R 0

(28)

subjected to the boundary conditions n+1 ρn+1 j,0 − ρj,1 = 0 and ∆r

n+1 ρn+1 j,S − ρj,S+1 = 0. ∆r

for the spherical interface case. The fully discretized equations (27) and (28) are linearized and solved using Newton’s method to obtain the density distribution ρn+1 on each time level. More details about using Newton’s method for DGT equations can be found in our algorithm paper 22 . The stationary state of the system is reached when either of the following stopping P criteria is satisfied: Sk=1 |ρnk − ρn+1 | ≤ εd , or |(An − An+1 )/An | ≤ εe , or ∆tn ≥ εt . The first k stopping criteria means that the density difference between the two time levels is smaller than the tolerance εd , the second one means that the relative change of the Helmholtz free energy between two time levels is smaller than the tolerance εe , and the third one means that

15



Page 16 of 32

the time step length ∆tn is larger than the tolerance εt . We use εd = 1e−3, εe = 1e−6 and εt = 1e4 in this work. The obtained equilibrium density distribution is an approximation of the solution of the original mass-conserved DGT model.

Applications of the mass-conserved DGT model Results for planar interface cases Before applying the mass-conserved DGT model in a closed system, we first conduct conventional DGT calculations in open systems, of which the grand potential energy is minimized at equilibrium. Figure 2(a) shows the density profile convergence history of pure hexane at 293.15 K in an open system. PC-SAFT parameters and influence parameter of hexane are summarized in Table. 1. Starting from the linear density distribution (red solid line), the density profile converges to its equilibrium state (blue solid line) in 110 time steps using the SDGT algorithm. At equilibrium, the system’s total molar amount per unit area Ni /As is calculated by integrating the equilibrium density profile along the r direction, which equals to 1.98960503e–5 mol/m2 on a 5 nm domain. Table 1: Model parameters of Hexane. Component

Mi [g/mol]

mi [1]

σi [Å]

i /kB [K]

vi · 1020 [J · m5 /mol2 ]

n-Hexane

86.177

3.0576

3.7983

236.77

35.575

Same amount of hexane is used as the mass input to the mass-conserved DGT model. The density profiles are obtained by solving Eqn. 27 on each time level, and the convergence history is displayed in Figure 2(b). From a linear density distribution (red solid line), the system evolves to its equilibrium state (blue solid line) in 130 time steps. At equilibrium, the system’s total molar amount per unit area is 1.98960502e–5 mol/m2 , which is less than 1e–13 mol/m2 mass loss compared with the initial condition. This means that the mass has been conserved accurately during the calculation. 16


8000

7000

7000

6000

6000 3

Density [mol/m ]

8000

5000 4000 3000 2000 1000 0 0

5000 4000 3000 2000

Initial density profile Equilibrium density profile

1

2 3 Distance [nm]

4

Initial density profile Equilibrium density profile

1000 0 0

5

(a) Open system.

1

2 3 Distance [nm]

4

5

(b) Closed system.

Figure 2: Density profile convergence history of hexane at 293.15 K in (a) an open system and (b) a closed system. 8000

Density [mol/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Density [mol/m3]

Page 17 of 32

6000

4000

Open system Closed system

2000

0 0

1

2

3

4

5

Distance [nm] Figure 3: Comparison of hexane’s equilibrium density profiles (293.15 K) of DGT models in open and closed systems. By comparing the two density profile convergence histories in Figure 2, we demonstrate the differences of density profile behaviors in µVT and NVT systems. The conventional DGT model in a µVT system fixes boundary densities, which guarantees constant bulk chemical potentials during the calculation. At the same time, the system is allowed to have mass exchange with the environment, which leads to the change of the overall mass (area under the density profile) from time to time. The mass-conserved DGT model in a

17



NVT system, on the other hand, retains the area under the density profile as a constant on every time step, while the density values on each position, including boundaries, are subjected to change. Since the mass-conserved DGT model used the same amount of mass from equilibrium state of the open system case, the two models produce identical equilibrium density profiles as shown in Figure 3. When applying the conventional DGT model in open systems, the phase equilibrium calculation is required. The obtained bulk densities are used to set up boundary conditions for DGT equations. For the mass-conserved DGT model, no phase equilibrium calculation is required, and the bulk properties as well as the interfacial properties can be obtained from equilibrium density profiles. The work-flows of the two models are compared and summarized in Figure 4. To verify the accuracy of the bulk densities calculations from the mass-conserved DGT model, we compare the results with those from the phase equilibrium calculations at different temperatures in Figure 5. The comparison shows that the massconserved DGT model reproduces the bulk densities accurately without explicitly performing bulk phase equilibrium calculations. Conventional DGT in open systems:

Mass-conserved DGT in closed systems:

Inputs: T (for pure system) T, P, zi (for mixed system)

Phase equilibrium calculation for bulk properties

Inputs: T, V, Ni

DGT calculation

DGT calculation Bulk properties and interfacial properties

Interfacial properties

Figure 4: Comparison of workflows of the conventional DGT model and the mass-conserved DGT model. From the phase envelope in Figure 5, it can be concluded that the vapor–liquid phase 18


Page 18 of 32

Page 19 of 32

550 VLE calculation Mass conserved DGT

500

Temperature [K]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


450 400 350 300 250 200 0

2000

4000

6000

8000

10000

3

Density [mol/m ]

Figure 5: The phase envelope of hexane generated by the phase equilibrium calculation (solid line) and the mass-conserved DGT model (square markers). Markers on the dash line indicate the average input densities as illustrated in Figure 6. splitting occurs when the average density ρavg falls inside the phase equilibrium curve: ρV < ρavg < ρL , where ρavg = Ni /(As D). Otherwise, only one vapor or liquid phase would exist. For example, the equilibrium bulk densities of hexane at 293.15 K are ρV = 6.7 and ρL = 7588.1 mol/m3 . Figure 6 shows the equilibrium density profiles with different input average densities, as those marked on the dash line of Figure 5. For ρavg = 3000 (black line), 4000 (blue line), and 5000 (red line) mol/m3 , the input average densities fall in the spinodal region and vapor–liquid phase split takes place. The volume of the liquid phase increases as more hexane is added to the system, while the interface shape and bulk densities do not change. If ρavg < ρV or ρavg > ρL , such as when ρavg = 5 (magenta line) or 10000 (green line) mol/m3 , no phase splitting occurs, and thus no interface appears.

Results for spherical interface cases As discussed in the introduction, it is technically impossible to obtain a stable droplet in an open system using the SDGT algorithm due to the metastable saddle point on the grand potential energy surface. In this section, we study the droplet formation by applying the 19



12000 10000

Density [mol/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

8000 6000

Average density [mol/m3]

10000 5000 4000 3000 5

4000 2000 0 0

1

2

3

4

5

Distance [nm] Figure 6: Equilibrium density profiles of hexane at 293.15 K with different mass in the system. mass-conserved DGT model in a closed spherical system. Figure 7 shows the formation of the hexane droplet in a system of radius R = 10 nm. A step function is used as initial density distribution (red solid line), which is equivalent to having a 5 nm homogeneous liquid droplet surrounded by homogeneous vapor phase. The total molar amount is 1.59542e–21 mol at initial time. The dash lines depict the convergence history of the density profiles. After 210 time steps the stopping criteria is met which indicates a thermodynamically stable state. The mass measured at the equilibrium is 1.59536e–21 mol, which indicates the mass loss of less than 6e–26 mol. Compared with the planar interface cases, finer grids are required here since numerical errors are amplified at nodes that are close to R. From the equilibrium density profile (blue solid line), it can be observed that the hexane forms a droplet with a radius of about 4 nm, and it has an interface thickness of about 2 nm, which is comparable to the droplet radius. At such a small scale with a diffusive interface, the traditional thermodynamic approaches such as Young–Laplace equation loose their accuracy, while DGT is more suitable for this scenario. The mass-conserved DGT calculation is also carried out at different temperatures using

20


Page 21 of 32

8000 Initial density profile Equilibrium density profile

Density [mol/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6000

4000

2000

0 0

2

4

6

8

10

Distance [nm]

Figure 7: Density profile convergence history of hexane droplet at 293.15 K. the same initial setup. The equilibrium density profiles are shown in Figure 8. As the temperature increases, the radius of the liquid droplet decreases and the interface becomes more diffusive. At the same time, the liquid density decreases while the vapor density increases at a higher temperature. This is due to the fact that molecules are energized at higher temperatures, and more liquid vaporizes into the vapor phase. These temperaturerelated behaviors are consistent with density functional theory calculations 16 . In a system with a planar interface, pressures are equal everywhere at equilibrium. Otherwise the pressure difference will lead to fluid flow until the pressure balance is achieved. In a system with spherical interface, however, there is a pressure difference across the interface region. Equilibrium liquid and vapor pressures are calculated on different radius droplets, and the results are displayed in Figure 9. The solid line shows the pressure as a function of the droplet radius at 293.15 K, and the dash line marks the pressure of a planar interface system at the same temperature, which is 0.1618 bar as computed by PC-SAFT EoS. For a droplet with small radius, both the liquid and vapor pressures are much higher than the planar interface pressure. As the droplet radius grows, the liquid and vapor pressures decrease, and will eventually drop to the same level as the planar interface pressure at infinitely large

21



10000 253 K 273 K 293 K 313 K 333 K 373 K 393 K

8000

Density [mol/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

6000

4000

2000

0 0

2

4

6

8

10

Distance [nm]

Figure 8: Equilibrium density profiles of hexane droplet at different temperatures. radius (infinitely small curvature). The reason that the pressure difference exists is because the IFT on the curved surface tends to contract the droplet, and therefore, the inner pressure increases to balance the force. With the equilibrium density profile, the IFT of a droplet can be calculated using the following equation: γ=

Ω + P A VA + P B VB , As

where Ω is the grand potential energy of the system, PA VA and PB VB are products of the bulk pressure with the bulk volume in the vapor and liquid phases respectively. For a spherical interface, a dividing surface is required to determine VA , VB and As . In this work, we use the Gibbs dividing surface 6 , which is the position where the surface excess number of particles Γs is equal to zero. For a spherical system, the zero adsorption is defined as: Z Γs = 4π

Re

r

2

Z

R

ρ(r) − ρA dr + 4π

r2 ρ(r) − ρB dr = 0 ,

Re

0

22


(29)

Page 23 of 32

from which, the position of the Gibbs dividing surface Re can be determined by Re3

=

3

RR 0

r2 ρ(r) − ρB dr . ρA − ρB

(30)

0.6

250

Spherical interface Planar interface

Spherical interface Planar interface

0.5

Pressure [bar]

200

Pressure [bar]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


150

100

0.4

0.3

0.2

50

0.1

0 0

50

100

0

150

50

100

150

Radius [nm]

Radius [nm]

(a) Liquid pressure

(b) Vapor pressure

Figure 9: Liquid and vapor phase pressures at 293.15 K and different radius calculated by the mass-conserved DGT model. Figure. 10 shows pressure as a function of bulk density at different radius. The liquid pressure is calculated by taking the liquid density at the center of the droplet into the PCSAFT EoS, and the vapor pressure is calculated based on the vapor density away from the interface region. As discussed in Figure 9, pressure increases as a result of decreasing radius. Higher pressures in/outside the droplet lead to denser liquid/vapor phases compared with their saturated densities at the same temperature. For more compressible vapor phase, there is an almost 5 times density difference due to the formation of the droplet. As the radius keeps decreasing, the droplet would eventually vanish, and there is only a homogeneous supersaturated vapor phase left as marked by the green dot in Figure 9b. The IFT of the hexane droplet at 293.15 K with different radii are calculated and displayed as the blue solid line in Figure 11. The black dotted line marks the IFT on a planar interface. From infinitely small radius, the system’s IFT increases rapidly to a maximum value. After that, the IFT starts to decrease, and will eventually drop to the same value as the planar interface IFT at an infinitely large radius. 23



0.6

250 Spherical interface Planar interface

Pressure [bar]

150

100

0.4 0.3 0.2

50

0 7500

Spherical interface Planar interface Homogeneous vapor

0.5

200

Pressure [bar]

7600

7700

7800

0.1

7900

5

10

15

20

25

30

Density [mol/m3]

3

Density [mol/m ]

(a) Liquid phase

(b) Vapor phase

Figure 10: Bulk pressure as a function of bulk density for liquid and vapor phases calculated by the mass-conserved DGT model at 293.15 K and different radius. 19.2 19

IFT [mN/m]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 32

Spherical interface Planar interface Tolman equation

18.8 18.6 18.4 18.2 18 0

50

100

150

Radius [nm]

Figure 11: IFT results of hexane droplet at 293.15 K using the mass-conserved DGT model (blue solid line) and Tolman equation (red dash line). We recall below the Tolman equation that has been widely used as a macroscopic thermodynamic approach to describe the relationship between the droplet IFT and droplet radius: γ(r) 1 = , γ∞ 1 + 2δ/r

(31)

in which, γ∞ is the IFT on the planar interface, and δ is the so-called Tolman length, which

24


Page 25 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


by definition is: δ = Re − Rs ,

(32)

where Rs is called the surface of tension, at which the formal derivative of the IFT to the radius equals to zero: dγ = 0. dr r=Rs

(33)

According to the Tolman equation, the droplet’s IFT is a monotonic function of the droplet radius, and the sign of the Tolman length δ has a significant effect on the behavior of the IFT function: with growing radius, the IFT increases if Tolman length is positive, and the IFT decreases if Tolman length is negative. We determine the value of δ at the point where the IFT reaches its maximum in DGT calculations, so that the monotonicity can be best captured 3 . For hexane at 293.15 K, Re = 2.9585 nm at IFT maximum, and Rs = 3 nm. Therefore, δ = −0.0415 nm using Eqn. 32. The red dash line in Figure 11 shows the Tolman equation results. It can be seen that the two approaches have excellent agreement at radius r > 10 nm. For smaller droplets, however, the Tolman equation results deviate from DGT results. The DGT model successfully captured the sharp drop of IFT at very small radii that has been confirmed by both MD simulation 32 and DFT calculations 15,33 . In general, the Tolman equation is a fast and reasonably accurate approach to determine the γ − r relationship for large droplet, and DGT can serve as a reliable helper tool to calculate the Tolman length. For very small droplets, DGT is significantly a more suitable approach for the calculation.

Conclusion In this work, we extended the conventional DGT model from open systems to closed systems, in which the number of particles is fixed. By implementing a mass conservation term λi and no-flux boundary condition, we were able to enforce a rigorous canonical ensemble system

25



in simulations. The mass-conserved DGT model has been proved to conserve mass while dissipating the Helmholtz free energy over time. With this new model, we have successfully calculated the interfacial properties in closed systems, in particular a spherical system in which nucleation process takes place. From the calculation results, we: (I) validated the mass conservation feature of the new model; (II) demonstrated that no explicit phase equilibrium calculations were needed in the new model, while it still delivered accurate equilibrium bulk density predictions; (III) showed that the pressure in the droplet is much higher than that in the planar interface system at the same condition, and will reduce to the same level at sufficiently large droplet radius; (IV) showed that the droplet IFT would increase sharply and then decrease as the radius grows from infinitely small to infinitely large, which indicated a negative Tolman length. For future work, we plan to compare the mass-conserved DGT model with the DFT model and the molecular dynamics simulation on droplet’s IFT calculations. Another topic of interest is to combine the mass-conserved DGT model with the modified DGT model for surfactant molecules 24 , so as to study micelle formation.

Notes The authors would like to thank Shell International Exploration & Production Inc. for their financial support. The authors would also like to thank Dr. Kai Langenbach for valuable discussions.

The first two authors contributed equally. In particular, Mu performed the model implementation and Frank derived the mathematical model (18).

26


Page 26 of 32

Page 27 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Solomon, S. Climate Change 2007 – The Physical Science Basis: Working Group I Contribution to the Fourth Assessment Report of the IPCC ; Cambridge University Press, 2007; Vol. 4. (2) Ghan, S. J.; Abdul-Razzak, H.; Nenes, A.; Ming, Y.; Liu, X.; Ovchinnikov, M.; Shipway, B.; Meskhidze, N.; Xu, J.; Shi, X. Droplet nucleation: Physically-based parameterizations and comparative evaluation. Journal of Advances in Modeling Earth Systems 2011, 3 . (3) Lau, G. V.; Ford, I. J.; Hunt, P. A.; Müller, E. A.; Jackson, G. Surface thermodynamics of planar, cylindrical, and spherical vapour-liquid interfaces of water. The Journal of Chemical Physics 2015, 142, 114701. (4) Young, T. III. An essay on the cohesion of fluids. Philosophical Transactions of the Royal Society of London 1805, (5) Laplace, P.-S. Traité de mécanique céleste; supplément au dixième livre, sur lâ action capillaire. Courcier, Paris 1806, (6) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Courier Corporation, 2013. (7) Thomson, S. P. The Life of William Thomson, Baron Kelvin of Largs; Macmillan, 1910. (8) Tolman, R. C. The effect of droplet size on surface tension. The Journal of Chemical Physics 1949, 17, 333–337. (9) Yasuoka, K.; Matsumoto, M. Molecular dynamics of homogeneous nucleation in the vapor phase. II. Water. The Journal of Chemical Physics 1998, 109, 8463–8470.

27



(10) Horsch, M.; Vrabec, J.; Bernreuther, M.; Grottel, S.; Reina, G.; Wix, A.; Schaber, K.; Hasse, H. Homogeneous nucleation in supersaturated vapors of methane, ethane, and carbon dioxide predicted by brute force molecular dynamics. The Journal of chemical physics 2008, 128, 164510. (11) Horsch, M.; Vrabec, J.; Hasse, H. Modification of the classical nucleation theory based on molecular simulation data for surface tension, critical nucleus size, and nucleation rate. Physical Review E 2008, 78, 011603. (12) Horsch, M.; Hasse, H.; Shchekin, A. K.; Agarwal, A.; Eckelsbach, S.; Vrabec, J.; Müller, E. A.; Jackson, G. Excess equimolar radius of liquid drops. Physical Review E 2012, 85, 031605. (13) Zeng, X. C.; Oxtoby, D. W. Binary homogeneous nucleation theory for the gas–liquid transition: A nonclassical approach. The Journal of Chemical Physics 1991, 95, 5940– 5947. (14) Talanquer, V.; Oxtoby, D. W. Density functional analysis of phenomenological theories of gas-liquid nucleation. The Journal of Physical Chemistry 1995, 99, 2865–2874. (15) Malijevsk` y, A.; Jackson, G. A perspective on the interfacial properties of nanoscopic liquid drops. Journal of Physics: Condensed Matter 2012, 24, 464121. (16) Li, Z.; Wu, J. Toward a quantitative theory of ultrasmall liquid droplets and vaporliquid nucleation. Industrial & Engineering Chemistry Research 2007, 47, 4988–4995. (17) Van der Waals, J. D. The equation of state for gases and liquids. Nobel lectures in Physics 1910, 1, 254–265. (18) Cahn, J. W.; Hilliard, J. E. Free energy of a nonuniform system. I. Interfacial free energy. The Journal of Chemical Physics 1958, 28, 258–267.

28


Page 28 of 32

Page 29 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(19) Miqueu, C.; Mendiboure, B.; Graciaa, C.; Lachaise, J. Modelling of the surface tension of binary and ternary mixtures with the gradient theory of fluid interfaces. Fluid Phase Equilibria 2004, 218, 189–203. (20) Miqueu, C.; Míguez, J. M.; Piñeiro, M. M.; Lafitte, T.; ; Mendiboure, B. Simultaneous application of the gradient theory and Monte Carlo molecular simulation for the investigation of methane/water interfacial properties. The Journal of Physical Chemistry B 2011, 115, 9618–9625. (21) Niño-Amézquita, O. G.; Enders, S. Phase equilibrium and interfacial properties of water+ methane mixtures. Fluid Phase Equilibria 2016, 407, 143–151. (22) Mu, X.; Frank, F.; Alpak, F. O.; Chapman, W. G. Stabilized density gradient theory algorithm for modeling interfacial properties of pure and mixed systems. Fluid Phase Equilibria 2017, 435, 118–130. (23) Qiao, Z.; Sun, S. Two-phase fluid simulation using a diffuse interface model with Peng– Robinson equation of state. SIAM Journal on Scientific Computing 2014, 36, B708– B728. (24) Mu, X.; Xi, S.; Alpak, F. O.; Chapman, W. G. Modified Density Gradient Theory for Surfactant Molecules Applied to Oil/Water Interfaces. Industrial & Engineering Chemistry Research 2018, 57 (22), 7643–7654. (25) Planková, B.; Hrub` y, J.; Vinš, V. Homogeneous droplet nucleation modeled using the gradient theory combined with the PC-SAFT equation of state. EPJ Web of Conferences. 2013; p 01076. (26) Vinš, V.; Planková, B.; Hrub` y, J. Surface tension of binary mixtures including polar components modeled by the density gradient theory combined with the PC-SAFT equation of state. International Journal of Thermophysics 2013, 34, 792–812.

29



(27) Langenbach, K.; Heilig, M.; Horsch, M.; Hasse, H. Study of homogeneous bubble nucleation in liquid carbon dioxide by a hybrid approach combining molecular dynamics simulation and density gradient theory. The Journal of chemical physics 2018, 148, 124702. (28) Gross, J.; Sadowski, G. Perturbed-chain SAFT: An equation of state based on a perturbation theory for chain molecules. Industrial & Engineering Chemistry Research 2001, 40, 1244–1260. (29) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New reference equation of state for associating liquids. Industrial & Engineering Chemistry Research 1990, 29, 1709–1721. (30) Rubinstein, J.; Sternberg, P. Nonlocal reaction–diffusion equations and nucleation. IMA Journal of Applied Mathematics 1992, 48, 249–264. (31) Allen, S. M.; Cahn, J. W. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metallurgica 1979, 27, 1085–1095. (32) Koga, K.; Zeng, X. C.; Shchekin, A. K. Validity of Tolman?s equation: How large should a droplet be? The Journal of chemical physics 1998, 109, 4063–4070. (33) Talanquer, V.; Oxtoby, D. W. Density functional analysis of phenomenological theories of gas-liquid nucleation. The Journal of Physical Chemistry 1995, 99, 2865–2874.

30


Page 30 of 32

Page 31 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nomenclature a0

Homogeneous Helmholtz free energy density (J m−3 ).

A

Helmholtz free energy (J).

δ

Tolman length (m).

γ

Interfacial tension (J m−2 ).

J

Time interval, J := (0, ts ) (s).

λi

Mass-conservation enforcing term (J mol−1 ).

µ0i

Homogeneous chemical potential (J mol−1 ).

M

Mobility (assumed constant) (mol2 J−1 m−3 s−1 ).

Ni

Total molar amount of ith component (mol).

n

Outward unit normal on ∂V (1).

PA

Bulk pressure of phase A (Pa).

ρi

Molar density of ith component (mol m−3 ).

Rc

Critical nucleation radius (m).

Re

Radius of the Gibbs dividing surface (m).

Rs

Radius of the surface of tension (m).

S

Number of interpolation points, (1).

t

Time variable, t ∈ J (s).

tn

Time level n (s).

ts

Stationary state time level (s).

τn

Time step size (s).

vij

Influence parameter (constant, J m5 mol−2 ).

V

Physical domain, filled with mixture (m3 ).

∂V

Boundary of domain V (m2 ).

r

Space variable in 1D, r ∈ [0, D] (m).

r

Space variable, r ∈ V (m3 ).

DFT

Density functional theory.

DGT

Density gradient theory.

EoS

Equation of state.

IFT

Interfacial tension.

LLE

Liquid–liquid equilibrium.

SDGT Stabilized density functional theory. VLE

Vapor-liquid equilibrium.

31



Page 32 of 32

Vapor phase Diffusive interface

Liquid droplet

Figure 12: Table of Content graphics.

32


Mass-conserved density gradient theory model for nucleation process

Recommend Documents