Document [original]

This version is available at https://doi.org/10.14279/depositonce-8390

This work is licensed under a CC BY-NC-ND 4.0 License (Creative

Commons Attribution-NonCommercial-NoDerivatives 4.0 International). For

more information see https://creativecommons.org/licenses/by-nc-nd/4.0/.

Chatwell, René Spencer; Heinen, Matthias; Vrabec, Jadran (2019). Diffusion limited evaporation of a

binary liquid film. International Journal of Heat and Mass Transfer, 132, 1296–1305.

https://doi.org/10.1016/j.ijheatmasstransfer.2018.12.030

René Spencer Chatwell, Matthias Heinen, Jadran Vrabec

Diffusion limited evaporation of a binary

liquid film

Accepted manuscript (Postprint)Journal article |

Diusion limited evaporation of a binary liquid lm

René Spencer Chatwell

, Matthias Heinen

, Jadran Vrabec

∗

Thermodynamics and Process Engineering, Technical University of Berlin, 10587 Berlin,

Germany

Abstract

An analytical solution of a model uid's time behavior, known as the Ste-

fan problem, is presented. A scenario is investigated in which a planar two-

component liquid lm is continuously evaporating into a thermodynamically

non-ideal vapor phase. Evaporation is initiated and maintained by a spatial

chemical potential gradient, while its rate is limited by the components' diu-

sion uxes across the vapor-liquid interface. Local thermodynamic equilibrium

is found to be present throughout the process. In contrast to the classical ap-

proach relying on equations of state, all required non-idealities are formulated in

relation to the Gibbs energy and are determined by molecular simulations. Ini-

tially, the liquid is an equimolar mixture of two components of dierent volatil-

ity, whereas the adjacent vapor phase is dominated by a dense inert gas. To

validate the analytical model and verify all exploited assumptions, the results

are contrasted to large scale molecular dynamics simulations.

Keywords:

evaporative mass transfer, analytical modeling, large scale

molecular dynamics simulations

1. Introduction

Investigations of phase transitions for spherical and planar geometries are

methodically similar. Historically, research was focused on droplet evaporation

indicating a relationship between the current square diameter

d2(t)

and an initial

∗

vrabec@tu-berlin.de

This is the Accepted Manuscript of: Chatwell, R. S., Heinen, M., & Vrabec, J. (2019). Diffusion limited evaporation of

a binary liquid film. International Journal of Heat and Mass Transfer, 132, 1296–1305.

https://doi.org/10.1016/j.ijheatmasstransfer.2018.12.030

This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International

License, http://creativecommons.org/licenses/by-nc-nd/4.0/.

value

to exist, which for pure liquids is linear in time

d2(t) = d2

0−Kt .

(1)

For droplets immersed into a quiescent atmosphere, the evaporation rate

is exclusively a function of the thermodynamic state. The literature refers to

this equation as "d-squared law" [1]. The problem's more thorough theoretical

treatments originated from around the early twentieth century [2, 3], whereas

a rst estimation of

had already been conducted by Maxwell [4] in his 1877

work on "Diusion" and was presented in its current form by Fuchs [5]

K=8Dρv

ρlye,s −ye,∞,

(2)

where

ye,i

represent the mole fraction of evaporate at either surface or in suf-

cient distance from the droplet. This mole fraction dierence initiates and

maintains the evaporation process while a diusion coecient

in combina-

tion with the vapor-liquid density ratio

ρvρ−1

determines its rate. Maxwell's

approach was later rened, among others by Spalding [6], to account for sub-

stantial evaporation rates that are indicative of an evaporate's high vapor pres-

sure and results in its accumulation in the interface region between liquid and

vapor. Consequently, a more detailed description of the droplet's surface com-

position became indispensable and necessitated Spalding to take the quiescent

atmospheric gas' mole fraction

ya,s

into account

K=8Dρv

ρl

ln1 + ye,s −ye,∞

ya,s .

(3)

In fact, Maxwell's description (2) represents the rst-order approximation of

Spalding's formulation (3) for the case of a dominating atmospheric gas at

droplet surface, i.e.

ya,s ≈1

. Maxwell and Spalding assumed both bulk phases

to be constantly in vapor-liquid equilibrium that due to mass transport is only

sustainable via continuous evaporation.

The d-squared law in its form (1) has been experimentally validated over

a substantial range of thermodynamic states as well as for a variety of non-

sooting monocomponent droplets through a wider spectrum of chemical com-

plexity [7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. Multicomponent droplets, however,

oer a conspicuously dierent behavior owed to the continuous variation in mo-

lar composition of both liquid

x(t)

and coexisting vapor

y(t)

. Species with

genuinely dissimilar properties lead to evaporation rates that vary substantially

over time [17, 18, 19]. Complementing theoretical and experimental work, mul-

ticomponent evaporation has also been successfully addressed by atomistic sim-

ulations [20, 21].

2. Object under study

Complementary to the literature's focus on spherical droplets, this work

solves the Stefan problem [22, 23] for a multicomponent liquid lm analytically.

This analytical solution, representing the d-squared law's (1) analogue for planar

surfaces, is then validated by large scale molecular dynamics (MD) simulations.

Consequently, a two-component liquid lm's change in molar composition

x(t)

and thickness

ξ(t)

is investigated during evaporation into an inert gas dominated

vapor phase. Each phase is composed of components with deliberately specied

characteristics, cf. appendix A. Due to these characteristics, in combination

with the selected thermodynamic state, outlined in table 1, both bulk phases

are thermodynamically non-ideal. An adequate description of the establish-

ing evaporation dynamics requires all component's diusion coecients, liquid

phase activity coecients and vapor phase fugacities to be known as functions

of the present state.

Molecular simulations are employed for two dierent reasons. For one,

the sampling of all required thermodynamic observables by dedicated MD and

Monte Carlo (MC) simulations requires substantially less approximation than

modeling resting on equations of state. On the other hand, the present hy-

drodynamic formulation addresses the phases' molar composition

x(t),y(t)

di-

rectly. Large scale MD simulations oer an alternative access to determine a

phase's molar composition with a spatial and temporal resolution that is orders

of magnitude greater than what is currently achievable by experiment. To be

comparable to hydrodynamic length and time scales, both a suciently large

0z(t)

lz(t)

sz(t)

δz

I II III IV

Figure 1: Snapshot taken from one of the present large scale MD simulations depicting the

scenario in vapor-liquid equilibrium, i.e. prior to evaporation. The studied domain is de-

composed into four regions that are connected physically on microscopic length scales. The

central role is attributable to the interphase (II), characterized by a steep density decline rang-

ing from liquid (I) to vapor (III). A chemical potential gradient is established by substituting

evaporated particles that originated from the liquid phase by inert gas particles within the

control volume (IV). The scenario's symmetry ensures the net molar ux to be zero at the

origin

z= 0

, while selective thermostatisation establishes a spatially and temporally constant

temperature

T= 80

K which is accompanied by a pressure of

p= 6.34

MPa for the present

mixture.

ensemble and a substantial simulated time is necessary. Figure 1 outlines an

atomistic representation of the investigated scenario containing

2.5·105

parti-

cles with an initial liquid lm and vapor phase thickness of

zl(t0) = 15

nm and

δ=zδ(t)−zs(t) = 30

nm, respectively. Various separate simulation series were

carried out, to either sample thermodynamic observables or to pursue large scale

MD simulations, cf. Appendix A.

The present hydrodynamic formulation rests on a uid domain decompo-

sition into four separate regions, as outlined in gure 1. The liquid phase (I)

initially consists of two equimolarly mixed components (1 and 2) with dierent

volatility, that in this context should be understood as an attribute for the com-

ponent's endeavor to transition from the liquid and to remain in the coexisting

vapor phase (III). The latter was chosen to be dominated by an inert gas (com-

ponent 3), which was specied not to be condensable and be barely miscible in

the liquid lm. Both bulk phases are physically coupled by an interface region

(II) that is addressed as

interphase

in the following. A control volume (IV)

representing the invariant atmospheric conditions completes the uid domain.

2.1. Local thermodynamic equilibrium

Under vapor-liquid equilibrium the particles' entire motion is thermal and

consequently drift-free Maxwellian. The applied chemical potential gradient ini-

tiates evaporation dynamics by introducing a collective drift, i.e hydrodynamic

contribution

, to the

-th component's particle velocity

, that superimposes

to the random thermal motion

vi=ui+wi.

(4)

Local thermal equilibrium exists in every spatial domain (I) - (IV) if all parti-

cles' thermal velocities are Maxwellian distributed. The corresponding velocity

distribution function

fz,i

for the relevant mass transport direction

is most

sensibly displayed in its contracted form [24, 25]

fz,i =rmi

2πkbTexph−mivz−uz2

2kBTi.

(5)

Particle interactions identify as mechanism that drive a system towards local

equilibrium [26]. A dense vapor phase is indicative of such dominating interac-

tions, realized in the present scenario by a high system pressure

p= 6.34

MPa

ensuring the particles' thermal velocities to be Maxwellian distributed natu-

rally in each domain, even the narrow interphase. To conrm the anticipated

Maxwellian distribution (5), additional stationary large scale MD simulations

for the liquid lm being prepared and maintained at

x†= (0.48,0.52)

mol mol

−1

were carried out. The results are presented in gure (2) and fully conrm local

equilibrium to exist. Similar results were reported in the literature on interfa-

cial mass transfer, for both stationary [27, 28] and instationary [29, 30] cases.

Non-Maxwellian distributed particle velocities, in contrast, have conclusively

been demonstrated to establish during evaporation into rareed gas phases or

vacuum [31, 32, 33].

3. Hydrodynamic description

The present formalism models the vapor phase as a classical boundary layer

and predicts the lm's evaporation rate, i.e. its regression, as well as the change

i yi

ρv,i

ϑi

ρ0,i

γi

φi

mol mol

−1

mol dm

−3

- mol mol

−1

mol dm

−3

- mol cm

−2

−1

1 0.13 1.10 4.09 0.5 47.77 1.07 1.41

2 0.03 0.23 20.07 0.5 42.33 1.06 0.13

3 0.84 7.12 - 0 8.33 - -

Table 1: The system was initially prepared in a thermodynamic state with temperature

T= 80

K and pressure

p= 6.34

MPa, where both phases are under vapor-liquid equilibrium (VLE)

with the given liquid

and vapor phase mole fractions

as well as partial vapor

ρv,i

and

pure component molar densities

ρ0,i

. The thermodynamic non-idealities, i.e. the fugacity

coecient's excess contribution

ϑi

, the activity coecient

γi

and consequently the diusion

ux's entire thermodynamic contribution

ϕi

, albeit functions of composition, were considered

constant and evaluated under these VLE conditions.

of both phases' molar composition over time, while resting exclusively on particle

conservation formulated in integral form.

A domain's change in mass is determined by an imbalance of particle

uxes

jz,i

through its boundary [34] and yields for the bulk liquid's domain (I)

z∈[0, zl(t)]

, while invoking the inert gas characteristics not to be condensable

and be barely miscible in this liquid

dt ZVt

dV 



ρl,1(z, t)

ρl,2(z, t)

!+I∂Vt

dA 



jz,1(z, t)

jz,2(z, t)

=



0

.

(6)

A component's particle density decomposes into mole fraction

and overall

molar density

ρl

of the respective phase

ρl,i =xiρl

, whereas the particle ux

factorizes to

jz,i =xiρluz,i

, with the convective molar averaged velocity

uz,i

towards the interphase. Assuming spatial homogeneity of the liquid lm's molar

composition

xi6=xi(z)

and density

ρl6=ρl(z)

with the additional observation

that the particle uxes are invariant across the domain's surface leads to a

straightforward integration, where

V(t) = zl(t)A0

dt



x1(t)ρl(t)zl(t)

x2(t)ρl(t)zl(t)

+ρl(t)



x1(t)uz,1(z, t)

x2(t)uz,2(z, t)

=



0

.

(7)

Figure 2: Contrasting the results between drifting Maxwell distributed (solid line), as dened

in (5), and the sampled MD velocities for the relevant volatile (triangles) and less-volatile (cir-

cles) components does not disclose any disparities. These data were evaluated at temperature

T= 80

K and three dierent positions ranging from bulk liquid to bulk vapor. The large

dierence between both component's distribution functions is due to their dierent molar

mass. a) Values sampled in the bulk liquid phase at position

z≤zl(t)

in close vicinity to

the interphase. b) Values sampled in the interphase at position

z∈]zl(t), zs(t)[

. c) Values

sampled in the bulk vapor phase at position

z≥zs(t)

in vicinity to the interphase.

The interphase

z∈]zl(t), zs(t)[

balances the molar particle uxes between both

bulk phases, where

jz,i =xi(t)uz,i(zl, t)

indicates the

-th component's molar

ux exiting the bulk liquid phase and

z,i =yi(zs, t)uz,i(zs, t)

indicating the ux

entering the adjacent bulk vapor phase

−ρl(t)



x1(t)uz,1(zl, t)

x2(t)uz,2(zl, t)

+ρv



y1(zs, t)uz,1(zs, t)

y2(zs, t)uz,2(zs, t)

=



0

.

(8)

It is by bringing together (7) and (8) that the connection between the liquid

phase's mass and the uxes in the vapor phase (III)

z∈[zs(t), zδ(t)]

is achieved

dt



x1(t)ρl(t)zl(t)

x2(t)ρl(t)zl(t)

+ρl(t)



y1(zs, t)uz,1(zs, t)

y2(zs, t)uz,2(zs, t)

=



0

.

(9)

Although the inert gas is assumed not to be soluble in the liquid, its presence in

the interphase and vapor has to be accounted for. Each component's collective

drifting motion

yiui

, as dened in (4), is decomposed further into an advective

yiu

and a diusive contribution [35]

yiUi

yiui=yiu+Ui.

(10)

The evaporate's total convective ux is related to its total diusive ux, cf.

supplementary material, with all velocities being formulated in the mean molar

reference frame

i=1

yiuz,i =y1

Uz,1+y2

Uz,2,

(11)

and it is this correlation that renders the evaporation process diusion driven.

It is not necessarily trivial to assume the vapor phase mole fraction quotients

to be time invariant

bi3(z) = yi(z, t)

y3(z, t),

(12)

yet this simplication still allows to reasonably predict the vapor phase's molar

composition. To ease the notation

Lz,i =bi3Uz,i

is set for each component's

diusive ux.

3.1. Diusive motion

A general description of diusive motion in an inhomogeneous medium is

given by the Fokker-Planck diusivity law [36]. In the present case, this ux

has to be computed adjacent to the interphase, i.e at

z=zs(t)





Lz,1(zs)

Lz,2(zs)

=−∂z 



D11(y)D12(y)

D21(y)D22(y)

·



b13(z)

b23(z)

!zs

(13)

wherein the Fick diusion coecients

Dij

express the phenomenologically pos-

tulated proportionality between ux and respective mole fraction gradient [37].

For multicomponent systems, however, each diusive ux is the result of all

driving gradients. The diusion matrix, containing Fick's coecients

Dij

, is

separable into a kinetic and a thermodynamic contribution [38]





D11 D12

D21 D22 

=



B11 B12

B21 B22 



−1

·



Γ11 Γ12

Γ21 Γ22 

,

(14)

the former being the Maxwell-Stefan and the latter the thermodynamic factor

matrix, describing a mixture's departure from ideality and being composition

derivatives of the excess Gibbs energy [39]

. In this work, both matrices

were considered constant throughout the entire evaporation process. Assuming

the kinetic part to be invariant is justiable, as for the thermodynamic part,

depending strongly on composition, this simplication was conrmed by the

performed MD simulations, cf. supplementary material. Bringing together (9),

(11), (13) and (14) leads to

dt



x1(t)ρl(t)zl(t)

x2(t)ρl(t)zl(t)

−ρv



B11 B12

B21 B22 



−1

·



Γ11 Γ12

Γ21 Γ22 



×dz



b13(z)

b23(z)

zs

=



0

.

(15)

The diusion matrix couples both uxes, its coupling strength is measured by

how much the eigenvector matrix

diers from the identity matrix

, and can

spectrally be decomposed





B11 B12

B21 B22 



−1

·



Γ11 Γ12

Γ21 Γ22 

=T·



D10

0D2

·T−1.

(16)

The eigenvalues

have to be understood as eective diusion coecients that

map the weighted action of both driving forces onto each ux

T−1·dt



x1(t)ρl(t)zl(t)

x2(t)ρl(t)zl(t)

−ρv



D10

0D2

·T−1

×dz



b13(z)

b23(z)

zs

=



0

.

(17)

In the attempt to pursue rst order eects, the appearing spatial gradient

dz(bi3)

is linearized. The vapor domain

δ=zδ(t)−zs(t)

constitutes a classical boundary

layer problem with exemplary boundary conditions out of which the following

curvature-gradient correlation (18) is a statement of particle conservation, cf.

supplementary material

z=zs(t) : bi3(zs) = bi3,s , dzz(bi3)|zs=dz(bi3)|zs2,

z=zδ(t) : bi3(zδ) = bi3,∞, dz(bi3)|z∞= 0 .

(18)

It has already been established that Fick diusion originates from a spatial mole

fraction gradient comprising the stagnant gas presence, which historically arisen

has been termed mass transfer number [6]

Bi3=bi3,∞−bi3,s =yi,∞−yi,s

y3,s

(19)

The nonlinearity in the proposed boundary conditions (18) necessitates a dif-

ferent ansatz to approximate the dierential operator by a dierence quotient.

The simplest function that complies with this set of boundary conditions is a

third order polynomial, rst introduced by Pohlhausen [40] and later applied to

mass transfer problems by Spalding [41]

bi3(z)−bi3,s

Bi3

=αi3z−zs(t)

δ+βi3z−zs(t)

δ2

+γi3z−zs(t)

δ3.

(20)

A set of nonlinear algebraic equations emerges that interrelates those coecients

1 = αi3+βi3+γi3,

(21a)

0 = α2

i3−2

Bi3

βi3,

(21b)

0 = αi3+ 2βi3+ 3γi3,

(21c)

with two sets of solutions out of which only the rst yields physically sensible

boundary layer proles

αi3=2

Bi3±r1 + 3

2Bi3−1,

(22a)

βi3= 3 + 4

Bi3∓r1 + 3

2Bi3+ 1,

(22b)

γi3=2

Bi3±r1 + 3

2Bi3−1−2.

(22c)

The binomial theorem allows to approximate square roots to a convenient degree

of accuracy, given the fact the transfer number is small

Bi31

, and the

anticipated linearization is attained by truncating the series approximation to

rst order

r1 + 3

2Bi3≈1 + 3

4Bi3.

(23)

It is by considering this approximation in the coecients

αi3, βi3, γi3

that (20)

exposes the proper derivative's substitution by a dierence quotient

dz(bi3)|zs≈αi3

bi3,∞−bi3,s

δ≈2

3δ

yi,∞−yi,s

y3,s

(24)

The proportionality coecient

αi3

in front of the otherwise simple dierence

quotient is being mandated by the physically required boundary conditions.

The present scenario was set up such that an inert gas component dominates

the vapor phase, permitting to set

y3,s ≈1

even at the surface. Moreover, the

composition of the surrounding atmosphere, represented by the control volume,

was specied to be

yi,∞= 0

. Inserting the attained linearization into (17), while

assuming weak ux coupling

T≈I

, leads to

dt



x1(t)ρl(t)zl(t)

x2(t)ρl(t)zl(t)

+3ρv

2δ



D10

0D2

·



y1(zs, t)

y2(zs, t)

≈



0

.

(25)

During local equilibrium, liquid and vapor phase mole fractions

xi, yi

are entan-

gled and as such cannot be determined independently.

100

3.2. Thermodynamic non-ideality

The central element when formulating phase equilibria is the extensive Gibbs

energy

that is constructed by the weighted arithmetic mean out of every

component's chemical potential

µi

in the mixture and is reduced

g=G/(nRT)

with the universal gas constant

and the respective phase's overall mole number

g=X

xiµi(T, p, x).

(26)

The classical framework of irreversible thermodynamics oers to use its concepts

locally and to further understand the Gibbs energy as a function of both position

and time [42]. It has already been demonstrated that local diusive equilibrium

exists across the interphase, which is represented by the equality of the phase's

chemical potentials

liq

i=∂g

liq

∂xi(t)=∂g

vap

∂yi(z, t)=µ

vap

i, i ∈[1,2] .

(27)

A formulation that is true, due to the extensive Gibbs energy's scaling behavior

as a homogeneous function of rst-order in the respective mole numbers. The

ideal gas acts as the lower physical boundary of every vapor phase given a

suciently low pressure

. At elevated system pressure

, the transition to a

real gas mixing behavior leads to

vap

i=1

yiµ0

0,i(T, p0) + lnp

p0+ ln(yi) + lnϕi(T, p, y),

(28)

in which all non-ideal contributions to the chemical potential are summarized

into the fugacity coecient

ϕi(T, p, y) = ψi(T, p)ϑi(T, p, y),

(29)

that factorizes into a residual

ψi

and an excess part

ϑi

, fully covering the compo-

nent's non-ideal behavior and the deviation from ideal mixing at system pressure

, respectively. This form allows segregating the ideal from the non-ideal mixing

contribution to the vapor phase's Gibbs energy

vap

i=1

yiµ0,i(T, p) + ln(yi) + lnϑi(T, p, y).

(30)

Liquid phases should be described dierently. The components in their pure

liquid state at system pressure are chosen to be the reference, while the activity

coecient

γi

accounts for all non-ideal mixing behavior

liq

i=1

xiµ0,i(T, p) + ln(xi) + lnγi(T, p, x).

(31)

The combination of (27), (30) and (31) exposes the correlation between liquid

and vapor phase composition, given both are in local equilibrium

yi(zs, t) = xi(t)γi(T, p, x)

ϑi(T, p, y).

(32)

It is evident from the excess parts' functional dependence on their respective

phase composition that a proper substitution of the mole fractions

requires

γi

and

ϑi

to be constant, which has already been utilized in (15). In-

troducing (32) into (25), together with the diusion ux's entire thermodynamic

contribution

φi=3ρvγi

2δϑi

Di,

(33)

leads to the substitution of the vapor phase mole fractions in favor of the liquid

phase's.

3.3. Analytical solution

A set of non-linear dierential equations arises that can be solved by quadra-

ture

dt



x1(t)ρl(t)zl(t)

x2(t)ρl(t)zl(t)

+



φ10

0φ2

·



x1(t)

x2(t)

≈



0

.

(34)

Multicomponent behavior is not inevitably the mere consequence of pure com-

ponent properties. The simultaneous presence of all species generates excess

contributions, which can inuence the thermal and caloric observables dier-

ently. A multicomponent liquid may behave with respect to its thermal prop-

erties often times nearly ideal, the density is then a sum of pure component

molar densities

ρ0,i

proportional to the current molar composition

xi(t)

, cf.

supplementary material

ρl(t)−1≈

i=1

xi(t)ρ0,i(T, p)−1.

(35)

Vanishing volumetric excess does thereby in no way entail ideal caloric behavior.

105

Applying the derivative and utilizing the assumed ideal behavior (35) reveals a

seemingly proper decoupled system

ρ2

lzl



ρ−1

0,20

0ρ−1

0,1

·



˙x1

˙x2



+



ρl˙zl+φ10

0ρl˙zl+φ2

·



x2

≈



0

,

(36)

where the liquid's total molar composition enforces a constraint on the mole

fractions

x1(t) + x2(t)=1.

(37)

Two thermodynamic parameters can be identied that determine the evapo-

ration rate with the thermodynamic non-ideality quotient

being the more

decisive

=ρ0,2

ρ0,1

, λ =φ2

φ1

(38)

Multiplying the rst line of (36) with

ρ−1

0,1

and the second line with

ρ−1

0,2

allows

to fully utilize the constraint (37) and yields upon summation of the resulting

equations

˙zl(t) + ρ−1

0,2φ2x1(t)

λ−1+ 1= 0 .

(39)

Similarly, the multiplication of the rst line of (36) with

x2(t)

and the second

one with

x1(t)

yields upon subtraction

˙x1(t)−φ2

ρ0,2zl(t)x3

1(t)1−+x2

1(t)−2+x1(t)

×1−1

λ= 0 .

(40)

The lm thickness

zl(t)

has consistently been formulated as a function of time

and so has the liquid mole fraction

x1(t)

. The division of (39) by (40) dis-

closes the alternative formulation of lm thickness as a function of composition

zl(x1(t))

dzl(t)

dx1(t)=zl(t)x1(t)1−

λ−1

1−1

λx3

1(t)1−+x2

1(t)−2+x1(t).

(41)

This functional relationship will be resolved successively, where the description

of lm thickness as a function of composition

zl(x1)

is established rst and the

correlation between mole fraction

and time

subsequently. It is expedient

to factorize the denominator polynomial

dzl(t)

dx1(t)=zl(t)

1−

λ−1x1(t)1−

λ−1

x1(t)x1(t)−1x1(t)−1

1−, , λ 6= 1 ,

(42)

with the exclusion of pure component behavior, being represented by a mixture

110

of identical particles

, λ = 1

. While all denominator roots are evidently real

valued, the use of a partial fraction decomposition conveniently partitions the

occurring rational function into a primitive sum

1−

λ−1x1(t)1−

λ−1

x1(t)x1(t)−1x1(t)−1

1−=λ

1−λx1(t)

λ−1x1(t)−1

x1(t)−1

1−

(43)

A variable separation of (42) in combination with the given decomposition (43)

renders the emerging equation integrable

115

Zξ(t)

ξ(t0)

dξ1

ξ(t)=Zx1(t)

x1(t0)

dx1 λ

1−λx1(t)+1

λ−1x1(t)−1

x1(t)−1

1−!,

(44)

with a non-dimensional lm thickness

ξ(t) = zl(t)zl(t0)−1

that allows to formu-

late a simpler initial condition

ξ(t0)=1 , x1(t0) = x1,0.

(45)

An integration reveals the anticipated description of lm thickness as a function

of composition and unfolds its physical interpretation as being the product of

an exponentially weighted composition and liquid density ratio

ξ(x1) = x1

x1,0λ

1−λx2,0

x2

1−λρl(t0)

ρl(t).

(46)

To lighten the notation and simplify further calculations, all invariant terms are

grouped into one constant

ξ(x1,0)

and the liquid density is formulated in terms

of molar composition

ξ(x1) = ξ(x1,0)xλ

1−x1

1−λx1−1

1−.

(47)

The description of liquid composition over time

x1(t)

cannot easily be con-

structed from already formulated functions, however, the inverse solution of

time as a function of composition

t(x1)

is readily attainable

dξ(x1(t))

dt =dξ(x1(t))

dx1(t)

dt ,

(48)

where the derivative on the left hand side is given in (39) and the rst derivative

on the right hand side can be computed from (47)

dξ

dx1

=ξ(x1,0)x

2λ−1

1−λ

11−x1λ−2

1−λx1+λ

1−λx1−1

1−+xλ

1−x1

1−λ.

(49)

Equations (39) and (49) contain negligible terms. A reduction in complexity is

easily obtainable by exploiting the strongly dierent volatilities

λ1

of the

two components in the liquid, which physically determines the more volatile

component to undergo phase transition preferentially, a characteristic that was

built into the particle model used for all simulations

120

d¯

dx1

≈ξ(x1,0)x1+λx1−1

1−

x11−x12+1

1−x1,

(50a)

d¯

dt ≈−φ1

ρ0,2zl,0x1+λ.

(50b)

The required relationship between composition and time is then determined

by an additional separation of variables ansatz in (48) resulting in a quotient

function

dt =Zx1

x1,0

dx1dξ

dx1

/dξ

dt ,

(51)

that in combination with the attained approximations (50a) and (50b) generates

an analytically solvable integral equation

dt =Zx1

x1,0

dx1 −ρ0,2ξ0zl,0

φ1x1+λx1−1

1−

x11−x12x1+λ

1−x1x1+λ!,

(52)

where the emerging rational function similarly partitions into a primitive sum,

given all denominator roots being real valued

x1+λx1−1

1−

x11−x12x1+λ=a1(, λ)

+a2(, λ)

1−x1

+a3(, λ)

1−x12+a4(, λ)

x1+λ,

(53a)

1−x1x1+λ=a5(, λ)

1−x1

+a6(, λ)

x1+λ.

(53b)

All occurring integration constants show the already proclaimed dependence on

125

the thermodynamic parameters

, λ

and do thereby signicantly contribute to

the lm's evaporation rate

a1(, λ) = 1

−1, a2(, λ) = λ2++λ+λ2

−1+λ2,

a3(, λ) = λ +

−1+λ, a4(, λ) = λ −−λ

+λ2,

a5(, λ) = a6(, λ) = 1

+λ.

(54)

As a result, the functional relationship between composition and time is decom-

posed into four linear terms that are straightforwardly integrable to yield

t(x1) = ρ0,2ξ0zl,0

φ1 a1lnx1,0

x1+a7lnx2

x2,0

+a3

x2−x2,0

x2x2,0

+a8lnx1,0+λ

x1+λ!,

(55)

with

a7=a2+a5

and

a8=a4+a6

The solution of a liquid lm's Stefan problem, i.e. the d-squared law's (1)

planar surface analogue for both cases, the two-component and pure liquid phase

is readily obtained. The former is given by the combination of (47) and (55)

relating lm thickness

ξ(x1)

to time

t(x1)

via the liquid's mole fraction

The latter is a special case for

, λ = 1

, that is not entailed in (47) and (55),

yet is alternatively obtained by integration of (39) with the addition of (32) to

substitute the vapor phase's non-ideality

ϑi

with the evaporate's mole fraction

ye,s

near the surface

zl(t) = zl,0−3Dρv

2δρl

ye,st .

(56)

Comparing pure component liquid lms to droplets oers a salient dierence

in the order to which both cases evaporate with respect to their characteristic

length

zl(t)

and

d(t)

. In this rst-order approximation (56) the liquid lm evap-

orates linearly in time, where the spherical droplet's analogue for the same level

of approximation, as given by Maxwell's equation (2), regresses quadratically

d2(t) = d2

0−8Dρv

ρl

ye,st .

(57)

A high evaporation rate, i.e. short liquid phase lifetime, is facilitated by a high

130

diusivity or low liquid density, respectively. For planar surfaces, the boundary

layer thickness

remains in the description as a nite size eect, which is in full

accordance with classical boundary layer theory. Small values of

correspond

to steep chemical potential gradients that inevitably lead to faster evaporation

rates.

135

4. Results

Under vapor-liquid equilibrium (VLE), the mass uxes between both bulk

phases balance and consequently no net mass transfer occurs. During evapora-

tion, the liquid lm is constantly forfeiting mass in the attempt to dispose of

the developed chemical potential gradient and restore a VLE. The establishing

140

evaporation dynamics aects each bulk phase's molar composition dierently.

Due to the volatile component's characteristic to evaporate preferentially, the

less-volatile component successively enriches within the liquid, as depicted in g-

ure 3 for the interphase's vicinity at

z=zl(t)

. The process causes the lm not

only to regress but also increases its overall density, since the less-volatile com-

145

ponent has a higher saturated liquid density, cf. table 1 and gure 4. The bulk

vapor phase's molar composition in vicinity to the interphase does marginally

change during evaporation. While the volatile component depletes quite rapidly,

the inert gas accumulates at the surface. The less-volatile's mole fraction, in

contrast, remains almost stagnant, cf. gure 5 evaluated at

z=zs(t)

150

The combination of (47) and (55) allows to compute a two-component lm's

dimensionless thickness

ξ(t)

as function of time, if the phases' non-idealities

ϑi, γi

and the components' Fick diusion coecients

Dij

are known. Similarly,

(56) determines a monocomponent lm's regression rate, given the evaporate's

diusion coecient

and molar composition adjacent to the interphase

ye,s

155

are known. All required thermodynamic data were sampled by molecular sim-

ulations, cf. Appendix A, and are listed in table 1, Appendix B and the

enclosed supplementary material. Three dierently composed liquid lms were

investigated. Ranging from pure volatile over an initially equimolarly composed

two-component mixture to pure less-volatile. The monocomponent lms' evap-

160

oration rates expectedly envelop those of all possible binary mixtures, cf. gure

6. The equimolarly prepared lm regresses initially quite rapidly and then,

due to a limitation embedded in (55) that predicts the volatile component's

full evaporation to be reached only asymptotically

limx1→0t(x1) = ∞

, transi-

tions into an apparently stagnant and nite thickness for small values of

, i.e.

165

x1≤0.01

mol mol

−1

. In contrast, the large scale MD simulations indicate a

transition into a constant rate that closely resembles, yet not equals, the pure

less-volatile component's evaporation rate

K(T, p, x)≈K2(T, p)

, since residual

volatile particles still remain in the vapor phase for an extended period of time.

Figure 3: Time evolution of the molar fractions of the volatile (red) and less-volatile (blue)

components in the liquid phase. The dashed lines represent the analytical solution and the

solid lines are large scale MD simulation data, evaluated on the liquid side of the interphase

zl(t)

. The liquid is successively enriched with the less-volatile component.

5. Conclusion

170

The Stefan problem for a multicomponent planar lm was solved. An analyt-

ical model was presented, describing evaporation dynamics of a liquid lm into

a dense non-ideal vapor phase. The outlined formalism requires specied bulk

phase non-idealities

γi, ϑi, ψi

, as well as Fick diusion coecients

Dij

at the

interphase's vicinity, i.e. at

z=zs(t)

. Consequently, its solutions (47), (55) and

175

(56) allow to determine the time evolution of both phases' molar compositions

x(t),y(t)

as well as the lm's thickness

ξ(t)

quantitatively. The components'

deliberately specied characteristics, i.e. strongly dierent volatilities

λ1

weak diusion ux coupling

T≈I

and the inert gas behavior, facilitate expe-

dient simplications in the derivation of the model and its solutions.

180

The limiting quantity of a two-component lm's evaporation rate

ξ(t)

identied as the volatile component's diusion ux

Lz,1

, wherein all driving

inuences, i.e. the activity coecient

γ1

, the fugacity coecient's excess contri-

bution

ϑ1

, the vapor phase thickness

and the eective diusivity

, appear

Figure 4: Time evolution of the partial density proles of both components in the liquid mix-

ture as sampled by the present large scale MD simulations. Two opposing density gradients

emerge, resulting in two counter-oriented currents, where the volatile component (red,1) prop-

agates towards the interphase and preferentially evaporates, while the less-volatile component

(blue,2) evaporates at a lower rate and partly regresses towards the liquid lm's center. A local

volatile component's enrichment within the interphase is observable. The respective overall

density is depicted in black. The liquid lm's density increase over time is a mixture eect

due to the enrichment of component 2 that has a higher saturated liquid density

(ρl,1< ρl,2)

Figure 5: Time evolution of the vapor phase mole fraction on the vapor side of the interphase

zs(t)

, where the dominating inert gas' behavior (cyan) is prominent.

Figure 6: Time evolution of dimensionless lm thickness

ξ(t) = zl(t)/zl(t0)

, where the dotted

lines represent the analytical solution while the solid lines the large scale MD simulations.

The evaporation rate is given by the slope of these curves, where the equimolar liquid case

(yellow) is enveloped by the pure volatile's (red) and pure less-volatile's (blue) rate

K1(T, p)>

K(T, p, x)> K2(T, p)

. The strong dierence between the evaporation rates of the two pure

uids

K1, K2

is intentional and a consequence of the chosen intermolecular force eld model.

linearly and consequently aect the lm's behavior equally.

185

The analytical model was intentionally contrasted to large scale MD sim-

ulations serving three purposes. First, the phases' molar compositions were

assessed in an unconstrained manner. Second, the existence of local thermo-

dynamic equilibrium was veried unambiguously in every domain and third it

was demonstrated that the hydrodynamic formalism can justiably be applied

190

to such a nanoscale scenario. The MD simulations fully conrm that neglect-

ing the lm's internal dynamics, by assuming spatial homogeneity

xi6=xi(z)

as in (6) and (7), does not compromise the accurate prediction of its composi-

tion

x(t)

and thickness

ξ(t)

over time. Consequently, the hydrodynamic eld

description is applicable to system sizes of the order of 50 nm, given the ther-

195

mophysical properties are known a priori. Atomistic simulations prove to be an

adequate methodology determining all required properties without constraints.

The symbiosis of both approaches oers a promising route for further studies

on evaporation.

Supplementary material

200

All calculated and sampled data are disclosed in this manuscript's supple-

mentary material.

Acknowledgements

The present work contributes to the Collaborative Research Center (SFB) 75

of Deutsche Forschungsgemeinschaft (DFG) and was funded under the grant VR

205

6/9-2. All computations were performed either on the HPC cluster

OCuLUS

the Paderborn Center for Parallel Computing (PC

) or on the Cray XC40 system

Hazel Hen

at the High Performance Computing Center Stuttgart (HLRS) with

resources allocated according to grant MMHBF2. This work was carried out

under the auspices of the Boltzmann-Zuse society.

210

0z(t2)

lz(t2)

sz(t2)

δz

I II III IV

0z(t1)

lz(t1)

sz(t1)

δz

I II III IV

Figure A.7: Snapshot taken from the present large scale MD simulations depicting two

successive instants of time

t2> t1

during evaporation. While the liquid lm (I) regresses, the

control volume (IV) was simultaneously expanded to keep the vapor phase's (III) thickness

constant

δ6=δ(t)

Appendix A. Molecular dynamics simulations

The model uid's thermodynamic behavior is determined by deliberately

specied Lennard-Jones parameters as outlined in table A.2. The component's

distinct volatilities were realized here by dierent energy parameters

i

. The in-

ert gas characteristic to be barely miscible in the liquid was modeled primarily

215

by the Berthelot parameter

ζi36= 1

All particles in the investigated scenario interact via the full Lennard-Jones

potential with

3.5σ

cut-o radius and analytic long range corrections [43]. The

chosen Lennard-Jones parameters approximately describe a liquid lm com-

posed of either nitrogen + oxygen or argon + krypton. The large scale MD

220

simulations were carried out in a cuboid volume with dimensions

Lx=Ly= 15

nm and

Lz= 90

nm, an integrator time step of

∆t= 5

fs and periodic boundary

conditions in

and

directions.

The present large scale MD simulations are, technically speaking, nonequi-

librium molecular dynamics (NEMD) simulations, which does not inevitably

225

lead to non-Maxwellian distributed particle velocities, but allows possible local

component

i σi

/nm

i

ζ1i

ζ2i

volatile 1 0.3 75 2 - 1

less-volatile 2 0.3 100 20 1 -

inert gas 3 0.3 10 30 0.3 0.3

Table A.2: The thermodynamic behavior of the uid sampled by simulation was determined

by this set of Lennard-Jones parameters that were chosen to produce a strongly non-ideal

vapor, yet retain a fairly ideal liquid phase. The length parameter

determines a particle's

eective diameter while the energy parameter



is a measure for a particle's dispersive inter-

action strength. Additionally, a modication of the Berthelot combining rule via

ζki

allows

to alter the interaction energy between unlike particle species. The inert gas was in a highly

supercritical state and consists of mainly repulsive interacting heavy particles. While all simu-

lations sampling thermodynamic observables were carried out with the open source molecular

simulation tool

2 [44], the large scale MD calculations were performed with the massively

parallel open source code

mardyn

[45].

equilibria to arise naturally, as demonstrated in section 2.1. In order to inves-

tigate evaporation dynamics that are driven exclusively by a spatial chemical

potential gradient the system was selectively thermostated. Consequently, the

less-volatile particles in the liquid

z≤zl(t)

and the inert gas particles in the

230

vapor

z∈[zs(t), zδ(t)]

were thermostated in

and

directions only, avoiding

articial interference with the investigated mass transport in the relevant

di-

rection.

The necessary concentration gradient to initiate and maintain evaporation

dynamics was realized by substituting all particles of components 1 and 2 by

235

inert gas particles (component 3) within the control volume, i.e. maintaining

a molar composition of

y3= 1

z > zδ(t)

. During evaporation, the control

volume was expanded to the same extent as the liquid lm regressed, keeping

the vapor phase thickness

δ6=δ(t)

constant, cf. gure A.7.

The symmetry of the uid system was exploited to generate considerably

240

smoother sampling proles by calculating the arithmetic mean between both

sides' density and molar composition proles.

Appendix B. Diusion properties

A thermodynamically non-ideal system can display a conspicuous disparity

between multicomponent and collective pure component behavior, as measured

by the Gibbs energy's excess contribution

ge(T, p, y) =

i=1

yiln(ϑi).

(B.1)

The fugacity coecient's excess contribution

ϑi

can alternatively be calculated

from the species' chemical potentials

µi

(Appendix C) and in contrast to the

conventional approach, utilizing equations of state, were computed here via

dedicated MC simulations

ln(ϑi) = µi−ln(yi)−µ0,i .

(B.2)

Non-ideal behavior, however, is not merely described by the excess Gibbs en-

ergy's numerical value but also by its derivatives [46], where the thermodynamic

245

factor is determined by the curvature

Γij =δij +yi ∂ge

∂yj∂yi

−∂ge

∂yi∂y3

−

k=1

yk∂ge

∂yk∂yj

−∂ge

∂yk∂y3!.

(B.3)

A discrete data set of

values for a varying molar composition

y= (y1, y2, y3)

was produced by comprehensive MC simulations (see supplementary material)

and the computation of the necessary derivatives in B.3 was made possible

by selecting an appropriate

model. An appropriate ansatz is a Margules

250

type polynomial that is the result of an empirically motivated Wohl's expansion

[47] truncated to desired order. As data generated by simulations were made

suciently precise and plentifully available, an advanced fourth-order ansatz

[48] became feasible

ge=y1y2y3A12 +A21 +A32 −C1y1−C2y2−C3y3

+y1y2A21y1+A12y2−E12y1y2

+y1y3A31y1+A13y3−E13y1y3

+y2y3A32y2+A23y3−E23y2y3,

(B.4)

where the rst term models the genuine ternary contribution with parameters

and the following three terms describe the respective binary subsystems' in-

uence with parameters

Aij, Eij

. All parameters

(Aij, Eij, Ci)

were determined

with a standard nonlinear least squares t to the respective excess Gibbs energy

data obtained from MC simulations.

The Fick diusion matrix

can be decomposed into a kinetic

B−1

and a

thermodynamic

contribution, which are functions of the vapor phase's molar

composition. It has been argued above that both matrices were considered con-

stant and will only be evaluated at equilibrium condition

y∗= (0.13,0.03,0.84)

mol mol

−1

Γ(y∗) = 



0.82 −0.18

−0.03 1.0

.

(B.5)

Determining the kinetic contribution

B−1

, however, is a two-step procedure.

255

The Maxwell-Stefan diusivities Ð

have to be computed rst, which for this

ternary mixture are readily accessible via MD simulations

13 = 171.0·10−9

−1,

23 = 58.0·10−9

−1,

12 = 34.5·10−9

−1,

(B.6)

and subsequently the matrix

with its diagonal

Bii

and o-diagonal elements

Bij

has to be calculated [39]

Bii =yi

k=16=i

(B.7a)

Bij =−yi1

−1

i3, i, j ∈[1,2] .

(B.7b)

The Maxwell-Stefan diusivities Ð

quantify the relative motion between the

th and

-th component in the ternary mixture and are incorporated reciprocally

in the description. Consequently, it is this matrix's inverse that is part of the

thermodynamic contribution to relate the driving thermodynamic force, i.e. the

concentration gradient to its corresponding diusion ux, and is given as

B−1(y∗) = 



155.78 2.62

25.02 53.69 

·10−9

−1.

(B.8)

The Fick diusion matrix

D=B−1·Γ

couples the simultaneous inuence

of both gradients to each molar diusion ux

Lz,i

, as specied in (13). A

spectral decomposition facilitates the decoupling by mapping both gradients'

inuence onto an eective diusion coecient

that is given as the system's

-th eigenvalue

D(y∗) = 



120.97 0

0 55.88 

·10−9

−1,

(B.9)

where the deviation of the eigenvector matrix

from the identity matrix

is a measure of coupling strength and the arguably present mild coupling is

neglected by setting

T=



1 0.35

0.26 1 

≈I.

(B.10)

The simplication introduced in (50a) and (50b) was based on the components'

strongly dierent volatilities

λ1

and the exclusion of a pseudo mixture

6= 1

λ=φ2

φ1

= 0.093 ,  =ρ0,2

ρ0,1

= 1.18 .

(B.11)

Appendix C. Thermodynamic non-ideality

260

The vapor phase's Gibbs energy could be decomposed into the sum of a

variety of physically interpretable contributions

vap

(T, p, y) =

i=1

0,i(T, p) + g

res

0,i (T, p) + yiln(yi) + ge

i(T, p, y),

(C.1)

with the rst term being the ideal gas contribution of the pure component, given

at system temperature and pressure, and the second its residual contribution.

The third term describes the ideal mixing term and the fourth its non-ideal

mixing correction. The Gibbs energy can also be decomposed into the weighted

sum of its components' chemical potentials, as stated in (26)

vap

(T, p, y) =

i=1

yiµi(T, p, y).

(C.2)

The ideal gas contribution is readily understood as the ideal gas chemical po-

tential

0,i =yiµ0

0,i

. The residual and excess contribution's mathematical form

has, historically arisen, been dened to be a weighted logarithmic function akin

to the ideal gas mixing entropy

vap

(T, p, y) =

i=1

yiµ0

0,i(T, p) + lnψi(T, p)+ ln(yi)

+ lnϑi(T, p, y).

(C.3)

Generally, all pure component contributions were grouped together to ease no-

tation

µ0,i(T, p) = µ0

0,i(T, p) + ln(ψi(T, p)) ,

(C.4)

which consequently leads to the vapor phase's Gibbs energy

vap

(T, p, y) =

i=1

yiµ0,i(T, p) + ln(yi) + ln(ϑi(T, p, y)).

(C.5)

Bringing (C.2) and (C.5) together leads to description (B.2) for the case of

known chemical potentials

µi

, which were sampled by here by dedicated MC

simulations

ln(ϑi) = µi−ln(yi)−µ0,i .

(C.6)

Alternatively, if the mixture's Gibbs excess energy is known, e.g. from a

model, the component's excess contribution to the fugacity coecient can be

computed [46]

ln(ϑi) = ge+∂ge

∂yi

−

k=1

∂ge

∂yk

(C.7)

References

265

[1] C. Law, Recent advances in droplet vaporization and combustion, Energy

Combust. Sci. 8 (1982) 171201.

[2] H. Morse, On evaporation from the surface of a solid sphere. Preliminary

note, Proc. Amer. Acad. Sci. 45 (1910) 363367.

[3] N. Fuchs, Über die Verdampfungsgeschwindigkeit kleiner Tröpfchen in einer

270

Gasatmosphäre, Phys. Z. Sowjetunion 6 (1934) 224243.

[4] J. C. Maxwell, in: W. D. Niven (Ed.), The Scientic Papers of James Clark

Maxwell, Vol. 2, Cambridge Libarary Collection, 2011, pp. 625646.

[5] N. Fuchs, Evaporation and Droplet Growth in Gaseous Media, Pergamon

Press, 1959.

275

[6] D. Spalding, The combustion of liquid fuels, Fourth Symp. of Comb. 4

(1953) 847864.

[7] S. Kumagai, T. Sakai, S. Okajima, Combustion of free fuel droplets in a

freely falling chamber, Proc. Combust. Inst. 13 (1971) 779785.

[8] S. Okajima, S. Kumagai, Further investigations of combustion of free

280

droplets in a freely falling chamber including moving droplets, Proc. Com-

bust. Inst. 15 (1975) 401407.

[9] C. T. Avedisian, J. C. Yang, C. H. Wang, On low-gravity droplet combus-

tion, Proc. R. Soc. Lond. 420 (1988) 183200.

[10] I. Gökalp, C. Chauveau, J. Richard, M. Kramer, W. Leuke, Observations

285

on the low temperature vaporization and envelope or wake ame burning

of n-heptane droplets at reduced gravity during parabolic ights, Proc.

Combust. Inst. 22 (1989) 20272035.

[11] C. Chauveau, I. Gökalp, Droplet combustion under microgravity condi-

tions, in: B. Kaldreich (Ed.), VIIIth European Symposium on Materials

290

and Fluid Sciences in Microgravity, Vol. 2, European Space Agency, Paris,

1992, pp. 467472.

[12] C. Chauveau, X. Chesnau, I. Gökalp, Burning characteristics of n-heptane

droplets under dierent regimes, in: 31st Aerospace Science Meeting, 1993,

pp. 824832.

295

[13] H. Nomura, Y. Ujiie, H. J. Rath, J. Sato, M. Kono, Experimental study

on high-pressure droplet evaporation using microgravity conditions, Proc.

Combust. Inst. 26 (1996) 12671273.

[14] V. Nayagam, J. B. Haggard, R. O. Colantonio, A. J. Marchese, F. L. Dryer,

B. L. Zhang, F. A. Williams, Microgravity n-heptane droplet combustion

300

in oxygen-helium mixtures at atmospheric pressure, AIAA J. 36 (1998)

13691378.

[15] C. Morin, C. Chaveau, I. Gökalp, Droplet vaporization characteristics of

vegetable oil derived biofuels at high temperatures, Exp. Therm Fluid. Sci.

21 (2000) 4150.

305

[16] H. Ghassemi, S. W. Baek, Q. S. Khan, Experimental study on evaporation

of kerosine droplets at elevated pressures and temperature, Combust. Sci.

Technol. 178 (2006) 16691684.

[17] V. Glushkov, O. Todes, V. Fedoseev, Evaporation kinetics of drops of binary

mixtures of organic liquids, Kolloidn Zh. 32 (1967) 678683.

310

[18] E. Ravindran, J. Davis, Multicomponent evaporation of single aerosol

droplets, J. Coll. Int. Sci. 85 (1982) 278288.

[19] G. Brenn, L. Deviprasath, F. Durst, C. Fink, Evaporation of acoustically

levitated multi-component liquid droplets, Int. J. Heat Mass Transf. 50

(2007) 50735086.

315

[20] S. Sumardiono, J. Fischer, Molecular dynamics simulations of mixture

droplet evaporation, in: J. Pagliarini, S. Rainieri (Eds.), Eurotherm Sem-

inar 77, Heat and mass transfer in food processing, Edizione ETS Parma,

Parma, 2005, p. 323.

[21] A. Frezzotti, L. Gibelli, D. A. Lockerby, J. E. Sprittles, Mean-eld kinetic

320

approach to evaporation of a binary liquid into vacuum, Phys. Rev. Fluids

3 (2018) 054001.

[22] J. Stefan, Über die Theorie der Eisbildung, Monatsh. f. Mathematik und

Physik 1 (1890) 16.

[23] J. Stefan, Über die Verdampfung und die Auösung als Vorgänge der Dif-

325

fusion, Ann. Phys. 277 (1890) 725747.

[24] W. Baumjohann, R. Treusmann, Basic Space Plasma and Physics, Imperial

College Press, 2012.

[25] A. Lofti, J. Vrabec, J. Fischer, Evaporation from a free liquid surface, Int.

J. Heat and Mass Transf. 73 (2014) 303317.

330

[26] R. Schunk, A. Nagy, Ionospheres: Physics, Plasma Physics, and Chemistry,

Cambridge University Press, 2009.

[27] B. Hafskjold, S. Kjelstrup, Criteria for local equilibrium in a system with

transport of heat and mass, J. Stat. Phys. 78 (1995) 463494.

[28] I. Inazoli, S. Kjelstrup, D. Bedeaux, J. Simon, Thermodynamic properties

335

of a liquidvapor interface in a two-component system, Chem. Eng. Sci. 65

(2010) 41054116.

[29] J. Xu, S. Kjielstrup, D. Bedeaux, A. Røsjorde, L. Rekvig, Verication of

Onsager's reciprocal relations for evaporation and condensation using non-

equilibrium molecular dynamics, J. Colloid Interface Sci. 299 (2006) 452

340

463.

[30] S. Kjelstrup, D. Bedeaux, I. Inazoli, J.-M. Simon, Criteria for validity of

thermodynamic equations from non-equilibrium molecular dynamics simu-

lations, Energy 33 (2008) 11851196.

[31] T. Ishiyama, T. Yano, S. Fujikawa, Molecular dynamics study of kinetic

345

boundary condition at an interface between argon vapor and its condensed

phase, Phys. Fluids 16 (2004) 28992906.

[32] A. Frezzotti, L. Gibelli, S. Lorenzani, Mean eld kinetic theory description

of evaporation of a uid into vacuum, Phys. Fluids 17 (2005) 012102.

[33] A. Frezzotti, Boundary conditions at the vapor-liquid interface, Phys. Flu-

350

ids 23 (2011) 030609.

[34] L. D. Landau, E. Lifshitz, Course of Theoretical Physics. Fluid Mechanics,

Pergamon Press, 1987.

[35] S. de Groot, P. Mazur, Non-Equilibrium Thermodynamics, Dover Press,

1984.

355

[36] B. van Milligen, P. Bons, B. Carreras, R. Sánchez, On the applicability of

Fick's law to diusion in inhomogeneous systems, Eur. J. Phys. 26 (2005)

913925.

[37] A. Fick, Über Diusion, Ann. Phys. 170 (1855) 5986.

[38] C. Truesdell, Mechanical basis of diusion, J. Chem. Phys. 37 (1962) 2336

360

2344.

[39] R. Taylor, R. Krishna, Multicomponent Mass Transfer, Wiley, 1993.

[40] K. Pohlhausen, Zur näherungsweisen Integration der Dierentialgleichung

der laminaren Grenzschicht, ZAMM 1 (1921) 252290.

[41] D. Spalding, Mass transfer in laminar ow, Proc. R. Soc. A 221 (1954)

365

7899.

[42] G. Lebon, D. Jou, J. Casas-Vázquez, Understanding Non-Equilibrium

Thermodynamics, Springer, 2008.

[43] J. Jane£ek, Long range corrections in inhomogeneous simulations, J. Phys.

Chem. B 110 (2006) 62646269.

370

[44] C. W. Glass, S. Reiser, G. Rutkai, S. Deublein, A. Köster, G. Guevara-

Carrion, A. Wafai, M. Horsch, M. Bernreuther, T. Windmann, H. Hasse,

J. Vrabec, ms2: A molecular simulation tool for thermodynamic properties,

new version release, Comp. Phys. Commun. 185 (2014) 33023306.

[45] C. Niethammer, S. Becker, M. Bernreuther, M. Buchholz, W. Eckhardt,

375

A. Heinecke, S. Werth, H. Bungartz, C. Glass, H. Hasse, J. Vrabec,

M. Horsch, ls1 mardyn: The massively parallel molecular dynamics code

for large systems, J. Chem. Theory Comput. 10 (2014) 44554464.

[46] R. Taylor, H. Koojiman, Composition derivatives of activity coecient

models (for the estimation of thermodynamic factors in diusion), Chem.

380

Eng. Commun. 102 (1991) 87106.

[47] J. Prausnitz, R. Lichtenthaler, E. de Azvedo, Molecular Thermodynamics

of Fluid-Phase Equilibria, 3rd Edition, Prentice Hall, 1999.

[48] E. Hála, J. Pick, V. Fried, O. Vilím, Vapor-Liquid Equilibrium, Pergamon

Press, 1967.

385

Diffusion limited evaporation of a multicomponent liquid film

René Spencer Chatwell, Matthias Heinen, and Jadran Vrabec1, a)

Thermodynamics and Energy Technology, University of Paderborn, 33098 Paderborn,

Germany

This supplementary material discloses all calculated

thermodynamic data for the present analytical solution

together with the gemodel relying on data sampled

by MD or MC simulations. Values from the analytical

solution are indicated by the subscript ana, values from

the gemodel are denoted by mod and simulated data are

specified by sim. In the following tables, all quantities

are understood as functions of either liquid [xi] = mol

mol−1or vapor phase [yi] = mol mol−1mole fractions, if

not otherwise specified. The Gibbs excess energy geand

all chemical potentials µiare given in non-dimensional

form, i.e. ge=Ge/nRTand µi= ¯µi/RT

All time-independent thermodynamic properties were

determined by MC simulations in the isothermal-isobaric

ensemble containing N= 2000 particles. Transport

properties were determined by MD simulations in the

canonical ensemble containing N= 4914 particles and

an integrator time step of ∆t= 0.88 fs. Values were

obtained after about 9·108time steps, which represents

a sampling time of 0.8 µs.

containing N= 4914 particles, with values being ob-

tained after about 9·108time steps, which represents a

sampling time of 0.8 µs.

The selected gemodel (fourth order Margules type

polynomial) was fitted with a standard nonlinear least

squares algorithm to minimize the root mean square error

between simulation data and the model around the equi-

librium vapor phase composition y∗= (0.13,0.03,0.84)

and yielded for the parameters Aij, Eij and Ci





–A12 A13

A21 –A23

A31 A32 –



=



–−39.9 1.7

4.2–3.3

3.8 0.7–



(1)





–E12 E13

– – E23

– – –



=



–−0.01 1.63

– – −1.82

– – –



,(2)







=



−50.6

−126.3

−40.9



.(3)

To describe the liquid phase’s non-idealities with greater

precision, an additional fit around x∗= (0.495,0.505,0)

was performed to yield different parameters for the same

a)jadran.vrab[email protected]

FIG. 1. Gibbs ternary plot illustrating compositions where

thermodynamic data were generated by MC simulations. Un-

der equilibrium, the vapor phase’s composition is given as

y= (0.13,0.027,0.843) mol mol−1(red bullet) while the liq-

uid phase’s is given as x= (0.495,0.505,0) mol mol−1(green

bullet). The dashed lines show the vapor and liquid binodals.

FIG. 2. Segment of the Gibbs ternary plot depicting differ-

ent paths for the vapor phase composition during evaporation

from simulation (blue) and analytic solution (violet).

model applied to the binary subsystem only





A21

A12

E12



=



0.26

0.25

−0.01



.(4)

Vapor phase convective motion

The linear momentum of a multicomponent mixture is

constructed out of the vector sum of its constituent’s

linear momenta

ρvuz=

∑

i=1

ρu,ivz,i .(5)

The i-th component’s partial density is given by ρv,i =

yiρv. The third component is seen as a stagnant gas

uz,3≈0, i.e. no collective motion in a spatial direc-

tion is measurable, although random thermal motion is

present. The molar averaged velocity uzin the relevant

zdirection is then given as a proportionate sum of every

component’s convective velocity vz,i

uz=

∑

i=1

yiuz,i ≈

∑

i=1

yiuz,i .(6)

The convective motion uz,i can be decomposed into an

advective and a diffusive contribution uz,i =uz+Uz,i

∑

i=1

yiuz,i =

∑

i=1

yi(uz+Uz,i).(7)

Substituting the molar averaged motion uzwith the re-

sult already obtained from (6) yields

∑

i=1

yiuz,i =

∑

i=1

yi(

∑

i=1

yiuz,i +Uz,i),(8)

which after a simple algebraic manipulation connects

each component’s convective yivz,i to its diffusion flux

yiUz,i

∑

i=1

yiuz,i(1−

∑

i=1

yi)=

∑

i=1

yiUz,i .(9)

In a last step, the mixture’s total molar composition con-

straint y1+y2+y3= 1 is utilized to retain the description

used in equation (11) of the manuscript

∑

i=1

yiuz,i =

∑

i=1

(yi

y3)Uz,i =

∑

i=1

bi3Uz,i .(10)

Boundary condition

In equation (18) of the manuscript, a boundary

condition was used that needs a more detailed clarifi-

cation concerning its limitations. Akin to the particle

density conservation equation

∂t(ρv,i(z, t))+∂z(ρv,i(z, t)vz,i(z, t))= 0 ,(11)

a conservation equation for the weighted mole fractions

bi3can be constructed

∂t(bi3(z))+∂z(bi3(z)vz,i(z, t))= 0 ,(12)

where utilizing assumption (12) of the manuscript re-

duces the complexity significantly. It should be noted

that the parameters bi3are indeed functions of time,

yet their variation in time was considered negligible. A

formulation for this derivative at the interphase z=zs

needs to be found

∂z(bi3(z)vz,i(z, t))



zs

= 0 .(13)

Decomposing the convective velocity vector vz,i =uz+

Uz,i into its advective and diffusive parts again yields

∂z(bi3(z)(uz+Uz,i))



zs

= 0 .(14)

The diffusion flux has already been identified to be de-

scribed by Fick’s diffusion law that with the effective

diffusion coefficient Direads

bi3Uzi=−Didz(bi3(z)),(15)

after neglecting bi3dz(uz)

zsfor being small compared

to the other terms

uzdz(bi3(z))



zs

− Didzz(bi3(z))



zs

= 0 .(16)

In a last step, the connection between this vapor phase

argumentation and the liquid phase’s surface regression

velocity needs to be found. The latter is identified as the

negative convective motion of evaporate at the interphase

and is given

uz,liq(zs) = −uz(zs) =

∑

i=1

Didz(bi3(z))



zs

.(17)

After inserting (17) into (16) and expanding the sum, a

further assumption has to be made

(D1dz(b13(z))

zs+D2dz(b23(z))

zs)dz(bi3(z))

zs

−Didzz(bi3(z))

zs= 0 .(18)

All crossover terms are assumed to be negligible over the

non-crossover terms

Djdz(bj3(z))

zsdz(bi3(z))

zs≪ Di(dz(bi3(z))

zs)2

,∀j=i .

(19)

One can finally conclude the problem’s boundary condi-

tion at film’s surface to be

(dz(bi3(z))

zs)2

=dzz(bi3(z))

zs.(20)

TABLE I. Pure component simulation results that give the behavior of each component’s chemical potential µ0,i and residual

contribution ψias a function of pressure at constant temperature T= 80 K. The simulations yielded a uniform molar density

of ρ= 0.0151 mol dm−3for each component prepared in the state T= 80 K, p= 0.01 MPa that is precisely the ideal gas

equation’s result and determines the ideal gas chemical potentials to be µ0

0=−(7.86,7.86,7.85). Since the residual contribution

ψi, according to (27), measures a component’s departure from the ideal gas behavior, the components 1 and 2 are liquid while

component 3 is gaseous for the system being prepared in the state T= 80 K, p= 6.34 MPa.

p/MPa µ0,1µ0,2µ0,3ψ1ψ2ψ3

0.01 -7.86 -7.86 -7.85 1 1 1

0.10 -5.56 -5.57 -5.55 0.990 0.983 1.002

1.00 -3.37 -4.80 -3.23 0.885 0.247 1.021

2.00 -3.05 -4.68 -2.51 0.764 0.120 1.044

3.00 -3.01 -4.65 -2.09 0.425 0.082 1.067

3.25 -3.00 -2.00 0.397 1.073

3.50 -2.99 -1.92 0.373 1.079

3.75 -2.98 -1.85 0.350 1.085

4.00 -2.97 4.62 -1.78 0.332 0.064 1.091

4.25 -2.96 -1.71 0.315 1.097

4.50 -2.95 -1.65 0.300 1.103

4.75 -2.94 -1.59 0.287 1.109

6.00 -2.93 -4.59 -1.53 0.275 0.053 1.115

5.25 -2.92 -1.48 0.265 1.122

5.50 -2.91 -1.42 0.255 1.128

5.75 -2.91 -1.37 0.246 1.134

6.00 -2.90 -1.33 0.238 1.141

6.25 -2.89 -1.28 0.230 1.147

6.34 -2.89 -4.56 -1.26 0.228 0.043 1.150

FIG. 3. Graphical representation of the individual pure component’s chemical potentials µ0,i as a function of pressure at

constant temperature. The three curves, as they should, approach identical values for the ideal gas case.

TABLE II. The binary subsystem that is composed of volatile (component 1) and less-volatile (component 2) particles remains

liquid for all possible compositions. The fourth order Margules type polynomial covers also the respective infinite dilution

activity coefficients γ∞

i, as depicted in Fiq.(4).

x1x2y3ge

sim ge

mod γ1,sim γ2,sim γ1,mod γ2,mod

0 1 0 0 0 - 1 1.285 1

0.025 0.975 0 0.01 0.01 1.269 1.004 1.270 1.000

0.050 0.950 0 0.02 0.01 1.258 1.007 1.256 1.001

0.075 0.925 0 0.02 0.02 1.238 1.001 1.242 1.001

0.100 0.900 0 0.02 0.02 1.222 1.001 1.228 1.002

0.125 0.875 0 0.02 0.03 1.204 0.998 1.215 1.004

0.150 0.850 0 0.03 0.03 1.194 1.003 1.202 1.005

0.175 0.825 0 0.04 0.04 1.194 1.014 1.190 1.007

0.200 0.800 0 0.03 0.04 1.168 1.003 1.178 1.010

0.225 0.775 0 0.05 0.04 1.116 1.014 1.166 1.012

0.250 0.750 0 0.05 0.05 1.158 1.021 1.155 1.016

0.275 0.725 0 0.05 0.05 1.144 1.020 1.145 1.019

0.300 0.700 0 0.06 0.05 1.137 1.027 1.135 1.023

0.325 0.675 0 0.05 0.06 1.120 1.023 1.125 1.027

0.350 0.650 0 0.06 0.06 1.118 1.034 1.115 1.031

0.375 0.625 0 0.06 0.06 1.104 1.032 1.106 1.036

0.400 0.600 0 0.06 0.06 1.099 1.043 1.098 1.041

0.425 0.575 0 0.06 0.06 1.089 1.045 1.090 1.046

0.450 0.550 0 0.06 0.06 1.081 1.050 1.082 1.052

0.475 0.525 0 0.06 0.06 1.075 1.058 1.074 1.059

0.495 0.505 0 0.06 0.06 1.070 1.061 1.069 1.064

0.525 0.475 0 0.06 0.06 1.061 1.071 1.061 1.072

0.550 0.450 0 0.06 0.06 1.056 1.079 1.054 1.080

0.575 0.425 0 0.07 0.06 1.052 1.090 1.048 1.088

0.600 0.400 0 0.06 0.06 1.041 1.091 1.043 1.096

0.625 0.375 0 0.06 0.06 1.041 1.105 1.037 1.105

0.650 0.350 0 0.06 0.06 1.034 1.113 1.033 1.114

0.675 0.325 0 0.06 0.06 1.029 1.123 1.028 1.124

0.700 0.300 0 0.05 0.05 1.022 1.128 1.024 1.134

0.725 0.275 0 0.05 0.05 1.020 1.142 1.020 1.144

0.750 0.250 0 0.05 0.05 1.017 1.153 1.016 1.155

0.775 0.225 0 0.04 0.05 1.013 1.163 1.013 1.167

0.800 0.200 0 0.04 0.04 1.012 1.179 1.011 1.179

0.825 0.175 0 0.04 0.04 1.009 1.192 1.008 1.192

0.850 0.150 0 0.03 0.03 1.005 1.203 1.006 1.205

0.875 0.125 0 0.03 0.03 1.004 1.218 1.004 1.218

0.900 0.100 0 0.02 0.02 1.003 1.233 1.003 1.233

0.925 0.075 0 0.02 0.02 1.002 1.250 1.001 1.247

0.950 0.050 0 0.01 0.01 1.001 1.266 1.001 1.263

0.975 0.025 0 0.01 0.01 1.000 1.282 1.000 1.279

1 0 0 0 0 1 - 1 1.295

TABLE III. Data for the binary subsystem that is composed of volatile (component 1) and inert gas (component 3) particles.

Thermodynamic stability analysis relates the thermodynamic factor’s algebraic sign to either a stable (Γ>0) or an unstable

(Γ<0) state. This binary mixture must exhibit a vapor-liquid phase separation between 0.45 < y1<0.92. Both the Gibbs

excess energy and the fugacity coefficient ϑ1are well described in the vapor phase by the selected gemodel.

y1y2y3ge

sim ge

mod ϑ1,sim ϑ1,mod Γ1

0 0 1 0 0 - 5 1

0.010 0 0.990 0.02 0.02 5.10 5.16 0.98

0.020 0 0.980 0.03 0.03 5.02 5.07 0.96

0.025 0 0.975 0.04 0.04 4.98 5.03 0.96

0.030 0 0.970 0.05 0.05 4.94 4.98 0.95

0.040 0 0.960 0.06 0.06 4.87 4.90 0.93

0.050 0 0.950 0.08 0.08 4.79 4.82 0.92

0.060 0 0.940 0.10 0.10 4.72 4.74 0.91

0.070 0 0.930 0.11 0.11 4.65 4.67 0.89

0.080 0 0.920 0.13 0.13 4.59 4.60 0.88

0.090 0 0.910 0.14 0.14 4.52 4.54 0.87

0.100 0 0.900 0.16 0.16 4.45 4.47 0.86

0.150 0 0.850 0.23 0.23 4.15 4.18 0.81

0.200 0 0.800 0.30 0.30 3.88 3.93 0.75

0.250 0 0.750 0.36 0.36 3.63 3.68 0.67

0.275 0 0.725 0.39 0.39 3.52 3.56 0.63

0.300 0 0.700 0.42 0.42 3.40 3.44 0.58

0.350 0 0.650 0.48 0.48 3.20 3.20 0.46

0.400 0 0.600 0.53 0.53 3.00 2.95 0.32

0.450 0 0.550 0.57 0.57 2.81 2.70 0.16

0.500 0 0.500 0.61 0.60 2.63 2.45 -0.01

0.550 0 0.450 0.63 0.63 2.31 2.21 -0.18

0.600 0 0.400 0.64 0.64 2.01 1.98 -0.35

0.639 0 0.362 0.64 0.64 1.86 1.81 -0.46

0.650 0 0.350 0.63 0.64 1.83 1.76 -0.49

0.700 0 0.300 0.62 0.62 1.64 1.57 -0.60

0.750 0 0.250 0.59 0.58 1.50 1.41 -0.64

0.800 0 0.200 0.54 0.52 1.33 1.27 -0.60

0.850 0 0.150 0.46 0.44 1.23 1.15 -0.44

0.900 0 0.100 0.37 0.32 1.12 1.07 -0.15

0.910 0 0.090 0.35 0.30 1.10 1.06 -0.07

0.920 0 0.080 0.32 0.27 1.08 1.05 0.02

0.930 0 0.070 0.29 0.24 1.06 1.04 0.11

0.940 0 0.060 0.26 0.21 1.04 1.03 0.21

0.950 0 0.050 0.23 0.18 1.02 1.02 0.32

0.960 0 0.040 0.18 0.15 1.01 1.01 0.44

0.970 0 0.030 0.14 0.11 1.01 1.01 0.57

0.980 0 0.020 0.09 0.08 1.00 1.00 0.70

0.990 0 0.010 0.05 0.04 1.00 1.00 0.85

1 0 0 0 0 1 1 1

TABLE IV. Data for the binary subsystem that is composed of less-volatile and inert gas particles. This subsystem exhibits a

wider phase transition in terms of molar variation than the system in table (III). This subsystem must exhibit a vapor-liquid

phase separation in the range 0.175 < y2<0.8.

y1y2y3ge

sim ge

mod ϑ2,sim ϑ2,mod Γ2

0 0 1 0 0 - 35.34 1

0 0.010 0.990 0.03 0.04 26.33 32.93 0.95

0 0.020 0.980 0.07 0.07 25.71 30.74 0.90

0 0.031 0.969 0.10 0.11 25.06 28.55 0.83

0 0.040 0.960 0.13 0.14 24.54 26.91 0.78

0 0.043 0.957 0.14 0.15 24.31 26.50 0.77

0 0.050 0.950 0.16 0.17 23.98 25.23 0.72

0 0.100 0.900 0.32 0.33 21.40 18.69 0.41

0 0.150 0.850 0.47 0.47 19.13 14.28 0.10

0 0.175 0.825 0.54 0.53 18.07 12.62 -0.06

0 0.200 0.800 0.61 0.59 16.61 11.23 -0.21

0 0.255 0.775 0.62 0.65 9.23 10.05 -0.36

0 0.250 0.750 0.66 0.70 7.65 9.04 -0.49

0 0.300 0.700 0.74 0.78 6.09 7.42 -0.73

0 0.350 0.650 0.81 0.85 5.27 6.20 -0.91

0 0.400 0.600 0.87 0.89 4.53 5.25 -1.02

0 0.450 0.550 0.90 0.91 3.88 4.50 -1.07

0 0.500 0.500 0.93 0.91 3.40 3.89 -1.05

0 0.550 0.450 0.94 0.89 2.95 3.38 -0.96

0 0.600 0.400 0.93 0.85 2.56 2.95 -0.81

0 0.650 0.350 0.91 0.79 2.19 2.58 -0.60

0 0.700 0.300 0.88 0.72 1.89 2.26 -0.35

0 0.750 0.250 0.82 0.63 1.63 1.98 -0.07

0 0.800 0.200 0.75 0.52 1.44 1.73 0.22

0 0.850 0.150 0.64 0.40 1.30 1.51 0.49

0 0.900 0.100 0.51 0.28 1.13 1.32 0.74

0 0.950 0.005 0.31 0.14 1.04 1.15 0.92

0 1 0 0 0 1 1 1

FIG. 4. Binary liquid mixture’s (components 1 and 2) activity coefficients as given in table (III). The lines represent the

Margules polynomial model while the bullets represent simulation data. The activity coefficients at infinite dilution γ∞

(squares) are well covered by the chosen Margules type polynomial.

TABLE V. Data for the ternary system, where the first component’s mole fraction is varied and the second’s is kept constant

at y2= 0.025 mol mol−1.

y1y2y3ge

sim ge

mod ϑ1,sim ϑ2,sim ϑ1,mod ϑ2,mod Γ11 Γ12 Γ21 Γ22

0 0.025 0.975 0.08 0.08 - 25.41 4.96 25.48 1 0 -0.06 0.95

0.025 0.025 0.950 0.12 0.12 4.75 24.25 4.74 24.25 0.95 -0.05 -0.05 0.98

0.050 0.025 0.925 0.16 0.16 4.58 23.16 4.55 23.22 0.92 -0.09 -0.04 1.00

0.070 0.025 0.905 0.19 0.20 4.45 22.33 4.41 22.51 0.89 -0.11 -0.04 1.01

0.075 0.025 0.900 0.20 0.20 4.42 22.13 4.38 22.34 0.88 -0.11 -0.04 1.01

0.100 0.025 0.875 0.23 0.24 4.26 21.16 4.23 21.58 0.85 -0.15 -0.04 1.01

0.110 0.025 0.865 0.25 0.26 4.20 20.79 4.17 21.29 0.84 -0.16 -0.03 1.01

0.120 0.025 0.855 0.26 0.27 4.15 20.42 4.11 21.02 0.83 -0.17 -0.03 1.01

0.125 0.025 0.850 0.27 0.28 4.12 20.24 4.09 20.88 0.82 -0.18 -0.03 1.00

0.130 0.025 0.845 0.28 0.29 4.09 20.07 4.06 20.75 0.82 -0.18 -0.03 1.00

0.140 0.025 0.835 0.29 0.30 4.03 19.72 4.01 20.49 0.81 -0.19 -0.03 0.99

0.150 0.025 0.825 0.30 0.31 3.98 19.38 3.96 20.24 0.79 -0.20 -0.03 0.98

0.160 0.025 0.815 0.32 0.33 3.93 19.04 3.91 19.99 0.78 -0.21 -0.03 0.97

0.170 0.025 0.805 0.33 0.34 3.87 18.71 3.86 19.74 0.77 -0.23 -0.03 0.96

0.175 0.025 0.800 0.34 0.35 3.85 18.55 3.83 19.62 0.76 -0.23 -0.03 0.95

0.180 0.025 0.795 0.34 0.36 3.82 18.39 3.81 19.49 0.76 -0.24 -0.03 0.94

0.190 0.025 0.785 0.36 0.37 3.77 18.08 3.76 19.24 0.74 -0.25 -0.03 0.92

0.200 0.025 0.775 0.37 0.38 3.72 17.76 3.71 18.99 0.73 -0.26 -0.03 0.91

0.225 0.025 0.750 0.40 0.42 3.60 17.01 3.59 18.36 0.69 -0.30 -0.04 0.85

0.250 0.025 0.725 0.43 0.45 3.49 16.29 3.48 17.69 0.65 -0.35 -0.04 0.79

0.270 0.025 0.705 0.45 0.47 3.40 15.74 3.39 17.13 0.62 -0.38 -0.04 0.73

0.275 0.025 0.700 0.46 0.48 3.38 15.60 3.37 16.98 0.61 -0.39 -0.04 0.71

0.300 0.025 0.675 0.49 0.51 3.27 14.94 3.25 16.23 0.56 -0.45 -0.05 0.63

0.310 0.025 0.665 0.50 0.52 3.23 14.68 3.21 15.91 0.54 -0.47 -0.05 0.59

0.320 0.025 0.655 0.51 0.53 3.19 14.42 3.16 15.59 0.52 -0.49 -0.05 0.55

0.325 0.025 0.650 0.51 0.53 3.17 14.29 3.14 15.42 0.51 -0.51 -0.06 0.54

0.340 0.025 0.635 0.54 0.55 3.11 13.91 3.07 14.91 0.47 -0.54 -0.06 0.47

0.350 0.025 0.625 0.55 0.56 3.07 13.66 3.02 14.56 0.45 -0.57 -0.06 0.43

0.360 0.025 0.615 0.56 0.57 3.03 13.42 2.98 14.21 0.42 -0.60 -0.06 0.39

0.370 0.025 0.605 0.56 0.58 2.99 13.17 2.93 13.84 0.40 -0.60 -0.07 0.34

0.380 0.025 0.595 0.58 0.59 2.95 12.93 2.88 13.47 0.37 -0.65 -0.07 0.29

0.390 0.025 0.585 0.59 0.60 2.92 12.69 2.84 13.09 0.34 -0.68 -0.07 0.24

0.400 0.025 0.575 0.61 0.60 2.88 12.45 2.79 12.71 0.32 -0.71 -0.08 0.19

0.425 0.025 0.550 0.64 0.62 2.79 11.85 2.67 11.73 0.25 -0.77 -0.09 0.06

0.475 0.025 0.500 0.65 0.66 2.59 10.44 2.44 9.72 0.10 -0.90 -0.10 -0.24

0.525 0.025 0.450 0.66 0.68 2.16 6.10 2.22 7.74 -0.06 -1.00 -0.12 -0.58

0.575 0.025 0.400 0.65 0.69 1.96 4.94 2.00 5.91 -0.20 -1.03 -0.15 -0.95

0.675 0.025 0.300 0.63 0.67 1.63 3.50 1.62 2.99 -0.40 -0.77 -0.19 -1.81

0.725 0.025 0.250 0.60 0.63 1.48 2.93 1.47 1.98 -0.42 -0.40 -0.22 -2.29

0.825 0.025 0.150 0.48 0.49 1.22 2.09 1.23 0.75 -0.17 1.11 -0.27 -3.35

0.875 0.025 0.100 0.38 0.37 1.11 1.76 1.16 0.43 0.14 2.36 -0.29 -3.94

0.885 0.025 0.090 0.36 0.35 1.09 1.70 1.15 0.38 0.23 2.66 -0.29 -4.06

0.895 0.025 0.080 0.33 0.32 1.07 1.64 1.14 0.34 0.32 2.97 -0.30 -4.18

0.905 0.025 0.070 0.30 0.29 1.05 1.59 1.13 0.30 0.42 3.31 -0.30 -4.30

0.915 0.025 0.060 0.27 0.26 1.04 1.54 1.12 0.26 0.52 3.66 -0.31 -4.43

0.925 0.025 0.050 0.24 0.23 1.02 1.48 1.12 0.23 0.63 4.03 -0.31 -4.55

0.935 0.025 0.040 0.19 0.20 1.01 1.44 1.11 0.21 0.76 4.42 -0.32 -4.68

0.955 0.025 0.020 0.10 0.13 1.00 1.35 1.11 0.16 1.03 5.27 -0.33 -4.94

0.965 0.025 0.010 0.05 0.09 1.00 1.31 1.11 0.14 1.18 5.72 -0.33 -5.07

0.975 0.025 0 0.01 0.05 1.00 1.28 1.11 0.12 1.33 6.20 -0.33 -5.20

TABLE VI. Continuation of table (V) for y2= 0.050 mol mol−1.

y1y2y3ge

sim ge

mod ϑ1,sim ϑ2,sim ϑ1,mod ϑ2,mod Γ11 Γ12 Γ21 Γ22

0 0.050 0.950 0.16 0.16 - 23.98 4.78 23.87 1 0 -0.04 0.84

0.025 0.050 0.925 0.20 0.21 4.54 22.90 4.58 23.45 0.95 -0.02 -0.05 0.90

0.050 0.050 0.900 0.24 0.25 4.38 21.88 4.41 22.93 0.91 -0.06 -0.06 0.93

0.070 0.050 0.880 0.27 0.28 4.26 21.11 4.28 22.42 0.88 -0.07 -0.07 0.93

0.075 0.050 0.875 0.27 0.29 4.23 20.92 4.25 22.28 0.88 -0.07 -0.07 0.93

0.100 0.050 0.850 0.31 0.33 4.09 20.02 4.10 21.51 0.84 -0.17 -0.09 0.92

0.110 0.050 0.840 0.32 0.34 4.03 19.67 4.05 21.16 0.83 -0.20 -0.09 0.91

0.120 0.050 0.830 0.34 0.36 3.98 19.33 3.99 20.78 0.82 -0.24 -0.10 0.90

0.125 0.050 0.825 0.34 0.37 3.95 19.16 3.96 20.59 0.81 -0.25 -0.10 0.89

0.130 0.050 0.820 0.35 0.37 3.92 18.99 3.94 20.39 0.80 -0.27 -0.11 0.88

0.140 0.050 0.810 0.36 0.39 3.87 18.66 3.88 19.97 0.79 -0.31 -0.12 0.86

0.150 0.050 0.800 0.38 0.40 3.82 18.34 3.83 19.52 0.78 -0.34 -0.12 0.84

0.160 0.050 0.790 0.39 0.42 3.77 18.03 3.78 19.06 0.76 -0.38 -0.13 0.82

0.170 0.050 0.780 0.40 0.43 3.72 17.72 3.73 18.58 0.75 -0.42 -0.14 0.79

0.175 0.050 0.775 0.41 0.44 3.69 17.56 3.71 18.33 0.74 -0.45 -0.15 0.77

0.180 0.050 0.770 0.41 0.44 3.67 17.41 3.68 18.07 0.73 -0.47 -0.15 0.75

0.190 0.050 0.760 0.43 0.46 3.76 16.82 3.64 17.55 0.72 -0.51 -0.16 0.72

0.200 0.050 0.750 0.44 0.47 3.57 16.82 3.59 17.01 0.70 -0.56 -0.17 0.68

0.225 0.050 0.725 0.47 0.50 3.36 16.10 3.47 15.60 0.66 -0.69 -0.19 0.57

0.250 0.050 0.700 0.50 0.53 3.35 15.41 3.35 14.12 0.62 -0.82 -0.22 0.44

0.270 0.050 0.680 0.52 0.56 3.27 14.88 3.26 12.91 0.58 -0.93 -0.24 0.33

0.275 0.050 0.675 0.52 0.56 3.25 14.75 3.24 12.60 0.58 -0.96 -0.25 0.29

0.300 0.050 0.650 0.55 0.59 3.14 14.11 3.12 11.09 0.53 -1.10 -0.28 0.13

0.310 0.050 0.640 0.56 0.60 3.10 13.85 3.08 10.49 0.50 -1.16 -0.29 0.06

0.320 0.050 0.630 0.57 0.61 3.06 13.60 3.03 9.90 0.48 -1.21 -0.30 -0.02

0.325 0.050 0.625 0.58 0.62 3.04 13.48 3.01 9.61 0.47 -1.24 -0.31 -0.06

0.340 0.050 0.610 0.59 0.63 2.99 13.11 2.94 8.75 0.44 -1.33 -0.33 -0.18

0.350 0.050 0.600 0.60 0.64 2.95 12.86 2.89 8.19 0.42 -1.38 -0.34 -0.26

0.360 0.050 0.590 0.61 0.65 2.91 12.62 2.85 7.65 0.39 -1.43 -0.35 -0.35

0.370 0.050 0.580 0.62 0.66 2.87 12.38 2.80 7.12 0.37 -1.49 -0.37 -0.44

0.380 0.050 0.570 0.63 0.66 2.83 12.13 2.76 6.62 0.35 -1.54 -0.38 -0.53

0.390 0.050 0.560 0.63 0.67 2.80 11.88 2.71 6.13 0.32 -1.59 -0.39 -0.63

0.400 0.050 0.550 0.64 0.68 2.76 11.62 2.67 5.66 0.30 -1.63 -0.41 -0.72

0.450 0.050 0.500 0.66 0.71 2.35 7.13 2.45 3.64 0.18 -1.83 -0.48 -1.26

0.500 0.050 0.450 0.67 0.73 2.11 5.51 2.24 2.18 0.06 -1.92 -0.56 -1.86

0.550 0.050 0.400 0.68 0.73 1.92 4.59 2.04 1.20 -0.04 -1.86 -0.64 -2.54

0.600 0.050 0.350 0.67 0.72 1.77 3.93 1.86 0.61 -0.10 -1.61 -0.72 -3.28

0.650 0.050 0.300 0.65 0.69 1.62 3.35 1.70 0.29 -0.12 -1.11 -0.80 -4.08

0.700 0.050 0.250 0.61 0.64 1.46 2.81 1.57 0.12 -0.06 -0.28 -0.88 -4.95

0.750 0.050 0.200 0.56 0.56 1.33 2.39 1.46 0.05 0.10 0.93 -0.96 -5.89

0.800 0.050 0.150 0.49 0.46 1.21 2.03 1.39 0.02 0.38 2.60 -1.04 -6.88

0.850 0.050 0.100 0.39 0.32 1.11 1.72 1.35 0.01 0.81 4.82 -1.12 -7.94

0.870 0.050 0.080 0.34 0.26 1.07 1.61 1.35 0.00 1.03 5.87 -1.15 -8.38

0.880 0.050 0.070 0.31 0.23 1.05 1.56 1.35 0.00 1.16 6.44 -1.17 -8.60

0.890 0.050 0.060 0.28 0.19 1.03 1.51 1.35 0.00 1.29 7.04 -1.18 -8.82

0.900 0.050 0.050 0.25 0.16 1.02 1.46 1.35 0.00 1.43 7.66 -1.20 -9.05

0.910 0.050 0.040 0.20 0.12 1.01 1.41 1.36 0.00 1.58 8.31 -1.21 -9.28

0.920 0.050 0.030 0.16 0.08 1.00 1.37 1.37 0.00 1.74 8.99 -1.23 -9.51

0.930 0.050 0.020 0.11 0.04 1.00 1.33 1.38 0.00 1.90 9.70 -1.24 -9.75

0.940 0.050 0.010 0.06 0 1.00 1.29 1.39 0.00 2.08 10.45 -1.26 -9.98

0.950 0.050 0 0.01 -0.05 1.00 1.27 1.41 0.00 2.27 11.22 -1.27 -10.22

TABLE VII. Continuation of table (VI) for y2= 0.075 mol mol−1.

y1y2y3ge

sim ge

mod ϑ1,sim ϑ2,sim ϑ1,mod ϑ2,mod Γ11 Γ12 Γ21 Γ22

0.025 0.075 0.900 0.28 0.29 4.35 21.64 4.53 21.98 0.94 0 -0.03 0.76

0.050 0.075 0.875 0.31 0.33 4.19 20.69 4.32 21.80 0.89 -0.04 -0.07 0.78

0.070 0.075 0.855 0.34 0.36 4.08 19.97 4.16 21.41 0.85 -0.10 -0.11 0.78

0.075 0.075 0.850 0.35 0.37 4.05 19.80 4.12 21.27 0.84 -0.11 -0.12 0.78

0.100 0.075 0.825 0.38 0.41 3.92 18.94 3.95 20.41 0.80 -0.22 -0.17 0.75

0.110 0.075 0.815 0.40 0.43 3.86 18.62 3.88 19.98 0.78 -0.26 -0.19 0.73

0.120 0.075 0.805 0.41 0.44 3.81 18.29 3.81 19.50 0.76 -0.32 -0.21 0.70

0.125 0.075 0.800 0.41 0.45 3.79 18.14 3.78 19.24 0.76 -0.34 -0.22 0.69

0.130 0.075 0.795 0.42 0.46 3.76 17.98 3.75 18.97 0.75 -0.37 -0.24 0.67

0.140 0.075 0.785 0.43 0.47 3.71 17.67 3.69 18.40 0.73 -0.43 -0.26 0.64

0.150 0.075 0.775 0.45 0.49 3.66 17.36 3.62 17.80 0.71 -0.50 -0.28 0.60

0.160 0.075 0.765 0.46 0.50 3.62 17.06 3.57 17.15 0.70 -0.57 -0.31 0.56

0.170 0.075 0.755 0.47 0.52 3.57 16.77 3.51 16.48 0.68 -0.64 -0.33 0.51

0.175 0.075 0.750 0.48 0.52 3.55 16.62 3.48 16.14 0.67 -0.68 -0.34 0.49

0.180 0.075 0.745 0.48 0.53 3.52 16.48 3.45 15.78 0.67 -0.71 -0.35 0.46

0.190 0.075 0.735 0.50 0.54 3.48 16.19 3.39 15.06 0.65 -0.79 -0.38 0.41

0.200 0.075 0.725 0.51 0.56 3.43 15.91 3.34 14.33 0.63 -0.87 -0.40 0.35

0.225 0.075 0.700 0.54 0.59 3.32 15.23 3.21 12.46 0.59 -1.08 -0.47 0.18

0.250 0.075 0.675 0.56 0.61 3.22 14.56 3.08 10.59 0.55 -1.30 -0.54 -0.01

0.270 0.075 0.655 0.58 0.64 3.14 14.05 2.98 9.15 0.51 -1.48 -0.59 -0.19

0.275 0.075 0.650 0.59 0.64 3.12 13.92 2.96 8.80 0.51 -1.52 -0.61 -0.23

0.300 0.075 0.625 0.61 0.66 3.02 13.29 2.84 7.14 0.46 -1.74 -0.68 -0.48

0.310 0.075 0.615 0.62 0.67 2.98 13.04 2.80 6.52 0.44 -1.83 -0.71 -0.58

0.320 0.075 0.605 0.63 0.68 2.94 12.78 2.75 5.93 0.43 -1.92 -0.73 -0.69

0.325 0.075 0.600 0.64 0.68 2.92 12.66 2.73 5.65 0.42 -1.96 -0.75 -0.75

0.330 0.075 0.595 0.64 0.69 2.90 12.53 2.71 5.38 0.41 -2.00 -0.76 -0.81

0.340 0.075 0.585 0.65 0.70 2.86 12.21 2.66 4.85 0.39 -2.09 -0.79 -0.92

0.350 0.075 0.575 0.66 0.70 2.81 11.62 2.62 4.36 0.38 -2.17 -0.82 -1.04

0.360 0.075 0.565 0.66 0.71 2.74 10.87 2.58 3.91 0.36 -2.24 -0.85 -1.17

0.370 0.075 0.555 0.67 0.72 2.65 9.67 2.54 3.48 0.34 -2.32 -0.89 -1.30

0.375 0.075 0.550 0.67 0.72 2.59 8.93 2.51 3.28 0.33 -2.35 -0.90 -1.36

0.380 0.075 0.545 0.67 0.72 2.55 8.54 2.49 3.09 0.33 -2.39 -0.92 -1.43

0.390 0.075 0.535 0.67 0.73 2.47 7.71 2.45 2.74 0.31 -2.46 -0.95 -1.57

0.400 0.075 0.525 0.67 0.73 2.38 6.95 2.41 2.41 0.29 -2.52 -0.98 -1.71

0.425 0.075 0.500 0.68 0.74 2.25 6.06 2.32 1.72 0.26 -2.65 -1.06 -2.08

0.475 0.075 0.450 0.69 0.74 2.06 5.05 2.14 0.81 0.20 -2.79 -1.22 -2.89

0.575 0.075 0.350 0.68 0.70 1.74 3.73 1.84 0.13 0.21 -2.36 -1.55 -4.79

0.625 0.075 0.300 0.66 0.65 1.60 3.23 1.73 0.04 0.30 -1.64 -1.72 -5.88

0.675 0.075 0.250 0.63 0.58 1.45 2.73 1.65 0.01 0.48 -0.47 -1.89 -7.05

0.725 0.075 0.200 0.58 0.48 1.32 2.33 1.60 0.00 0.78 1.22 -2.06 -8.31

0.775 0.075 0.150 0.50 0.34 1.21 1.99 1.60 0.00 1.21 3.52 -2.22 -9.66

0.795 0.075 0.130 0.47 0.28 1.17 1.86 1.61 0.00 1.43 4.62 -2.29 -10.21

0.825 0.075 0.100 0.40 0.18 1.10 1.68 1.64 0.00 1.82 6.55 -2.38 -11.08

0.835 0.075 0.090 0.37 0.14 1.08 1.62 1.66 0.00 1.97 7.25 -2.42 -11.38

0.845 0.075 0.080 0.35 0.10 1.07 1.58 1.68 0.00 2.12 7.98 -2.45 -11.67

0.855 0.075 0.070 0.32 0.06 1.05 1.52 1.70 0.00 2.28 8.75 -2.48 -11.97

0.865 0.075 0.060 0.29 0.01 1.03 1.48 1.73 0.00 2.46 9.55 -2.51 -12.28

0.875 0.075 0.050 0.25 -0.03 1.01 1.43 1.76 0.00 2.64 10.39 -2.54 -12.58

0.905 0.075 0.020 0.12 -0.17 1.00 1.31 1.87 0.00 3.24 13.14 -2.63 -13.52

0.915 0.075 0.100 0.07 -0.22 1.00 1.28 1.91 0.00 3.46 14.13 -2.66 -13.84

0.925 0.075 0 0.02 -0.28 1.00 1.25 1.97 0.00 3.69 15.16 -2.69 -14.16

TABLE VIII. Spatially resolved vapor phase molar composition bi3data. The normalized form (bi3−bi3,s)/Bi3can be computed

with the following Spalding transfer numbers B13 = -0.154 and B23 = -0.003.

(z−zs)/δ b13(z)b23(z) (b13 −b13,s)/B13 (b23 −b23,s)/B23

0 0.155 0.032 0 0

0.01 0.152 0.031 0.016 0.025

0.03 0.147 0.030 0.048 0.076

0.05 0.142 0.028 0.079 0.125

0.07 0.138 0.027 0.111 0.174

0.09 0.133 0.025 0.142 0.222

0.11 0.128 0.024 0.173 0.269

0.13 0.123 0.022 0.204 0.314

0.15 0.118 0.021 0.234 0.358

0.17 0.114 0.020 0.264 0.401

0.19 0.109 0.018 0.294 0.443

0.21 0.105 0.017 0.323 0.483

0.23 0.110 0.016 0.352 0.521

0.25 0.096 0.015 0.381 0.558

0.27 0.091 0.013 0.409 0.594

0.29 0.087 0.012 0.437 0.628

0.31 0.083 0.011 0.465 0.660

0.33 0.079 0.010 0.492 0.691

0.35 0.071 0.009 0.518 0.720

0.37 0.067 0.008 0.544 0.747

0.39 0.063 0.008 0.570 0.772

0.41 0.061 0.007 0.595 0.796

0.43 0.059 0.006 0.619 0.819

0.45 0.056 0.005 0.643 0.840

0.47 0.052 0.005 0.666 0.859

0.49 0.048 0.004 0.689 0.876

0.51 0.045 0.004 0.711 0.892

0.55 0.038 0.003 0.753 0.920

0.57 0.035 0.002 0.773 0.932

0.59 0.032 0.002 0.792 0.943

0.61 0.029 0.002 0.811 0.952

0.63 0.027 0.001 0.828 0.961

0.65 0.024 0.001 0.846 0.968

0.67 0.021 0.001 0.862 0.974

0.69 0.019 0.001 0.877 0.980

0.71 0.017 0.005 0.892 0.984

0.75 0.013 0.000 0.919 0.991

0.77 0.011 0.000 0.931 0.993

0.79 0.009 0.000 0.942 0.995

0.81 0.007 0.000 0.952 0.997

0.83 0.006 0.000 0.961 0.998

0.85 0.005 0.000 0.970 0.999

0.87 0.004 0.000 0.977 0.999

0.89 0.003 0.000 0.984 1.000

0.91 0.002 0.000 0.989 1.000

0.93 0.001 0.000 0.993 1.000

0.95 0.005 0.000 0.997 1.000

0.97 0.002 0.000 0.999 1.000

0.99 0.000 0.000 1.000 1.000

1 0 0 1 1

FIG. 5. The normalized vapor phase molar composition was assumed to be a classical boundary layer problem and hence can

be considered analogous to the velocity distribution in a compressible boundary layer.

FIG. 6. Course of the weighted i-th component mole fraction bi3(z) = yi(z, t)/y3(z, t)as a function of position. Both values

approach bi3= 0 for z=zδdue to the enforced boundary condition.

TABLE IX. Data for the liquid xand vapor yphase molar compositions, as well as the film thickness’ time evolution. Since

the disparity between the analytical solution and the NEMD simulations ζsim < ζana became substantial, disclosing further

data after t > 213.50 ns was considered obsolete.

t/ns x1,sim x2,sim x1,ana x2,ana y1s,sim y2s,sim y3s,sim y1s,ana y2s,ana y3s,ana ζsim ζana

0 0.495 0.505 0.495 0.505 0.130 0.027 0.843 0.130 0.027 0.843 1 1

1.25 0.495 0.508 0.495 0.505 0.098 0.023 0.878 0.129 0.027 0.844 0.99 1.00

3.25 0.495 0.517 0.485 0.515 0.078 0.021 0.901 0.127 0.027 0.846 0.97 0.97

5.50 0.47 0.53 0.47 0.53 0.072 0.020 0.908 0.123 0.028 0.849 0.94 0.94

7.50 0.46 0.54 0.46 0.54 0.070 0.020 0.910 0.120 0.029 0.851 0.92 0.92

9.50 0.45 0.55 0.45 0.55 0.068 0.019 0.912 0.118 0.029 0.853 0.90 0.89

11.50 0.44 0.56 0.44 0.56 0.067 0.020 0.913 0.115 0.030 0.855 0.88 0.87

13.50 0.43 0.57 0.43 0.57 0.065 0.020 0.915 0.112 0.030 0.857 0.85 0.85

15.50 0.42 0.58 0.42 0.58 0.063 0.021 0.916 0.110 0.031 0.860 0.83 0.83

17.25 0.41 0.59 0.41 0.59 0.062 0.022 0.917 0.107 0.031 0.862 0.81 0.81

19.25 0.49 0.60 0.40 0.60 0.059 0.022 0.919 0.105 0.032 0.864 0.79 0.80

21.25 0.39 0.61 0.39 0.61 0.057 0.021 0.922 0.102 0.032 0.866 0.78 0.78

23.00 0.38 0.62 0.38 0.62 0.056 0.021 0.923 0.099 0.033 0.868 0.76 0.76

25.00 0.37 0.63 0.37 0.63 0.054 0.021 0.925 0.097 0.033 0.870 0.74 0.74

26.75 0.36 0.64 0.36 0.64 0.052 0.021 0.927 0.094 0.034 0.872 0.73 0.73

28.50 0.35 0.65 0.35 0.65 0.050 0.021 0.928 0.092 0.034 0.874 0.71 0.71

30.50 0.34 0.66 0.34 0.66 0.049 0.022 0.929 0.089 0.035 0.876 0.69 0.70

32.25 0.33 0.67 0.33 0.67 0.049 0.021 0.930 0.086 0.035 0.878 0.68 0.68

34.00 0.32 0.68 0.32 0.68 0.047 0.021 0.931 0.084 0.036 0.880 0.66 0.67

36.00 0.31 0.69 0.31 0.69 0.046 0.022 0.932 0.081 0.036 0.882 0.65 0.65

37.75 0.30 0.70 0.30 0.70 0.044 0.023 0.933 0.078 0.037 0.885 0.63 0.64

39.50 0.29 0.71 0.29 0.71 0.042 0.023 0.935 0.076 0.038 0.887 0.62 0.63

41.50 0.28 0.72 0.28 0.72 0.041 0.023 0.936 0.073 0.038 0.889 0.61 0.61

43.25 0.27 0.73 0.27 0.73 0.039 0.023 0.938 0.071 0.039 0.891 0.59 0.60

45.25 0.26 0.74 0.26 0.74 0.037 0.024 0.939 0.068 0.039 0.893 0.58 0.59

47.00 0.25 0.75 0.25 0.75 0.036 0.024 0.939 0.065 0.040 0.895 0.57 0.58

49.00 0.23 0.77 0.24 0.76 0.035 0.024 0.941 0.063 0.040 0.897 0.56 0.57

51.00 0.22 0.78 0.23 0.77 0.033 0.024 0.944 0.060 0.041 0.899 0.55 0.55

53.00 0.21 0.79 0.22 0.78 0.031 0.023 0.946 0.058 0.041 0.901 0.53 0.54

55.00 0.20 0.80 0.21 0.79 0.029 0.023 0.948 0.055 0.042 0.903 0.52 0.53

57.00 0.19 0.81 0.20 0.80 0.028 0.024 0.948 0.052 0.042 0.905 0.51 0.52

59.00 0.18 0.82 0.19 0.81 0.027 0.024 0.949 0.050 0.043 0.908 0.50 0.51

61.25 0.17 0.83 0.18 0.82 0.025 0.024 0.951 0.047 0.043 0.910 0.49 0.50

63.50 0.16 0.84 0.17 0.83 0.023 0.024 0.952 0.044 0.044 0.912 0.48 0.49

65.75 0.15 0.85 0.16 0.84 0.022 0.025 0.953 0.042 0.044 0.914 0.47 0.48

68.25 0.14 0.86 0.15 0.85 0.020 0.025 0.955 0.039 0.045 0.916 0.46 0.47

70.75 0.13 0.87 0.14 0.86 0.019 0.025 0.956 0.037 0.045 0.918 0.45 0.46

73.25 0.12 0.88 0.13 0.87 0.017 0.024 0.958 0.034 0.046 0.920 0.44 0.45

76.00 0.11 0.89 0.12 0.88 0.016 0.025 0.959 0.031 0.047 0.922 0.43 0.44

79.00 0.10 0.90 0.11 0.89 0.014 0.026 0.960 0.029 0.047 0.924 0.42 0.43

82.00 0.09 0.91 0.10 0.90 0.013 0.025 0.962 0.026 0.048 0.926 0.41 0.42

85.50 0.08 0.92 0.09 0.91 0.011 0.025 0.964 0.024 0.048 0.928 0.40 0.41

89.25 0.07 0.93 0.08 0.92 0.010 0.026 0.964 0.021 0.049 0.930 0.39 0.40

93.25 0.06 0.94 0.07 0.93 0.008 0.025 0.966 0.018 0.049 0.933 0.38 0.39

98.00 0.05 0.95 0.06 0.94 0.007 0.026 0.967 0.016 0.050 0.935 0.37 0.38

103.25 0.04 0.96 0.05 0.95 0.005 0.027 0.968 0.013 0.050 0.937 0.35 0.36

109.75 0.03 0.97 0.04 0.96 0.004 0.025 0.972 0.010 0.051 0.939 0.34 0.35

117.75 0.02 0.98 0.03 0.97 0.003 0.027 0.970 0.008 0.051 0.941 0.33 0.34

128.75 0.01 0.99 0.02 0.98 0.002 0.028 0.971 0.005 0.052 0.943 0.31 0.32

147.25 0.01 0.99 0.01 0.99 0.001 0.027 0.972 0.003 0.052 0.945 0.28 0.29

150.25 0.006 0.994 0.009 0.991 0.001 0.027 0.972 0.002 0.052 0.945 0.27 0.29

t/ns x1,sim x2,sim x1,ana x2,ana y1,sim y2,sim y3,sim y1,ana y2,ana y3,ana ζsim ζana

153.25 0.006 0.994 0.008 0.992 0.000 0.027 0.973 0.002 0.052 0.945 0.27 0.29

156.75 0.005 0.995 0.007 0.993 0.000 0.027 0.973 0.002 0.052 0.946 0.26 0.28

165.75 0.004 0.996 0.005 0.995 0.000 0.027 0.973 0.001 0.053 0.946 0.25 0.27

171.50 0.005 0.995 0.004 0.996 0.000 0.027 0.973 0.001 0.053 0.946 0.24 0.26

189.50 0.004 0.996 0.002 0.998 0.000 0.027 0.973 0.001 0.053 0.947 0.22 0.25

207.75 0.003 0.997 0.001 0.999 0.000 0.027 0.973 0.000 0.053 0.947 0.189 0.229

210.50 0.0038 0.996 0.0009 0.9991 0.000 0.028 0.972 0.000 0.053 0.947 0.185 0.226

213.50 0.0032 0.997 0.0008 0.9992 0.000 0.026 0.974 0.000 0.053 0.947 0.181 0.223

TABLE X. Comparing the deviation between the analytical solution liquid density time evolution to the NEMD simulation’s,

justifies the assumption that has been made in (33).

t/ns ρl,ana ρl,sim dev. in %

0.31 46.097 46.642 1.17

2.56 46.172 46.694 1.12

4.77 46.246 46.757 1.09

9.07 46.395 46.903 1.08

21.33 46.842 47.348 1.07

25.27 46.991 47.495 1.06

29.16 47.140 47.643 1.06

38.86 47.513 48.018 1.05

48.76 47.885 48.375 1.01

50.81 47.960 48.448 1.01

54.99 48.109 48.577 0.96

57.14 48.183 48.663 0.99

59.34 48.258 48.719 0.95

61.60 48.332 48.776 0.91

63.91 48.407 48.859 0.92

66.30 48.481 48.913 0.88

68.77 48.556 48.974 0.85

71.34 48.630 49.053 0.86

74.02 48.705 49.104 0.81

86.31 49.003 49.346 0.70

89.95 49.078 49.407 0.67

93.93 49.152 49.448 0.60

103.33 49.301 49.561 0.52

109.09 49.376 49.610 0.47

116.00 49.450 49.636 0.37

124.71 49.525 49.672 0.30

136.73 49.599 49.699 0.20

156.85 49.674 49.723 0.10

163.25 49.689 49.730 0.08

171.47 49.703 49.720 0.03

176.67 49.711 49.718 0.01

183.01 49.718 49.741 0.04

191.17 49.726 49.740 0.03

202.65 49.733 49.740 0.01