Document [original]

Droplet coalescence by molecular dynamics

and phase-field modeling

Cite as: Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131

Submitted: 22 January 2022 .Accepted: 10 March 2022 .

Published Online: 4 April 2022

Matthias Heinen,

Marco Hoffmann,

Felix Diewald,

Steffen Seckler,

Kai Langenbach,

and Jadran Vrabec

1,a)

AFFILIATIONS

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

Institute of Applied Mechanics, Technical University of Kaiserslautern, 67653 Kaiserslautern, Germany

Department of Optimization, Fraunhofer ITWM, 67663 Kaiserslautern, Germany

Scientiﬁc Computing in Computer Science, Technical University of Munich, 85748 Garching, Germany

Thermal Process Engineering, Innsbruck University, 6020 Innsbruck, Austria

Author to whom correspondence should be addressed: [email protected]

ABSTRACT

Coalescence of argon droplets with a radius of 25, 50, and 100 nm is studied with computational methods. Molecular dynamics (MD)

simulations are carried out to generate reference data. Moreover, a phase-ﬁeld model resting on a Helmholtz energy equation of state is

devised and evaluated by computational ﬂuid dynamics (CFD) simulations. Exactly the same scenarios in terms of geometry, ﬂuid, and state

are considered with these approaches. The MD and CFD simulation results show an excellent agreement over the entire coalescence process,

including the decay of the inertia-induced oscillation of the merged droplet. Theoretical knowledge about the asymptotic behavior of coales-

cence process regimes is conﬁrmed. All considered scenarios cross from the inertially limited viscous regime over to the inertial regime

because of the low shear viscosity of argon. The particularly rapid dynamics during the initial stages of the coalescence process in the thermal

regime is also captured by the phase-ﬁeld model, where a closer look at the liquid density reveals that metastable states associated with nega-

tive pressure are attained in the emerging liquid bridge between the coalescing droplets. This demonstrates that this model is even capable of

adequately handling the onset of coalescence. To speed up CFD simulations, the phase-ﬁeld model is transferred to coarser grids through an

interface widening approach that retains the thermodynamic properties including the surface tension.

Published under an exclusive license by AIP Publishing. https://doi.org/10.1063/5.0086131

I. INTRODUCTION

Droplet coalescence is ubiquitous in nature, e.g., during the forma-

tion of raindrops, and in technical applications, like viscous sintering

spray cooling.

Recently, Jin et al.

presented a general overview, discus-

sing ways to externally induce coalescence, which is a key factor for

designing lab-on-a-chip devices. Hence, there is an ongoing effort to

expand its understanding and to develop numerical models for the

design and optimization of applications. Therefore, coalescence has been

investigated by numerous theoretical,

4–9

experimental,

10–24

and numeri-

cal

8,9,21–36

approaches, and also review articles

37,38

have appeared.

The coalescence process between two droplets is characterized

by the formation of a bridge at their initial contact point that

develops into a meniscus at its outer boundary. According to the

Young–Laplace law, a strong pressure gradient emerges at this menis-

cus due to an initially strong negative curvature so that liquid ﬂows

in-between the droplets driven by the surface tension c, causing the

bridge to grow. The growth rate of the bridge also depends on the

velocity with which the liquid can ﬂow from the droplet interior into

the bridge. For the asymptotic limit of Stokes ﬂow, Hopper

provided

an analytical solution for the two-dimensional case of viscous cylin-

ders. Later, Eggers et al.

showed the applicability of Hopper’s solution

to the case of spherical droplets. However, they remarked that the

Stokes approximation is only valid up to a Reynolds number Re 1

so that for low shear viscosity g, the ﬂow is expected to become domi-

nated by inertial rather than viscous forces, once the transition regime

has been passed. The scaling laws derived from this analysis describe

some aspects of droplet coalescence dynamics. The process begins in a

viscous regime, where the time evolution of the bridge radius r

terms of the initial droplet radius Rapproximately satisﬁes rb=R

s=svln ðs=svÞfor rb=R<0:03 with the time scale sv¼gR=c.

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-1

Published under an exclusive license by AIP Publishing

Physics of Fluids ARTICLE scitation.org/journal/phf

In the case of sufﬁciently low shear viscosity, the ﬂow characteristics

will transition to an inertial (or inviscid) regime at later stages, where

rb=Rðs=siÞ1=2applies with the time scale si¼ðqR3=cÞ1=2. Here,

s¼tt0measures the time elapsed since the droplet contact at t

and qis the mass density.

While the 1/2 power law was conﬁrmed by many experimen-

tal

13–17

and numerical

9,29–32

approaches, inconsistencies were revealed

with respect to viscous scaling. Among others, Aarts et al.

found

their experimental data to obey linear scaling rb=Rs=sv.Also,

Paulsen et al.

—who employed an electrical method to probe the very

early stages of droplet coalescence (down to 10 ns after contact) in the

laboratory—could not corroborate the logarithmic correction

and

attributed the discrepancy to an inappropriate length scale of the

Reynolds number. They suggested the bridge height r2

b=Rinstead of

thebridgeradiusr

as the dominant length scale. Using this alternative

length scale deﬁnition, they managed to collapse their experimental

data covering a variation of shear viscosity over two orders of magni-

tude onto a single curve, pointing to the universality of that scaling

law. For this purpose, the data were rescaled by the crossover time s

and crossover radius r

obtained from the condition Re ¼1.

Subsequently, Paulsen et al.

identiﬁed the inertially limited viscous

(ILV) regime, upstream of the Stokes regime. They argued that the lat-

ter can only be attained once the inward force of the bridge generated

by the surface tension is sufﬁciently large to induce an asymptotic

acceleration of the coalescing droplet’s centers of mass. Consequently,

this condition leads to a threshold criterion for entering the Stokes

regime. Paulsen et al.

presented this in the form of a revised phase

diagram for coalescence, by plotting the Ohnesorge number

Oh ¼g=ﬃﬃﬃﬃﬃﬃﬃﬃ

qcR

pover the reduced bridge radius rb=R. This phase dia-

gram indicates that coalescence always begins in the ILV regime.

Depending on Oh, there is a subsequent crossover to either the Stokes

regime (Oh >1) or the inertial regime (Oh <1).

Reviewing the development of theory on coalescence dynamics,

it is noticeable that it has beneﬁted greatly from numerical investiga-

tions. For instance, Paulsen et al.

were able to corroborate their dis-

covery of the ILV regime by streamline analyses from computational

ﬂuid dynamics (CFD) simulations, showing qualitatively different ﬂow

patterns for all three regimes identiﬁed at that time (ILV, Stokes, and

inertial). Recently, Perumanath et al.

showed with molecular dynam-

ics (MD) simulation that another so-called thermal regime upstream

of the classical hydrodynamic regimes exists, where the bridge growth

dynamics is dominated by “collective molecular jumps.” MD simula-

tion is particularly suited for this regime, since the length and time

scales that this method can resolve are on the order of picometers and

picoseconds. Coalescence processes with their rapid dynamics, espe-

cially at the very beginning, can thus be studied with atomistic detail,

cf. Fig. 1 (Multimedia view). Another advantage is that the interface

and its properties, such as surface tension, evolve naturally during MD

simulation, without interventions or assumptions about its nature.

Residing on a lower level of abstraction, the sole input for MD simula-

tions is the molecular interaction model. It entails, however, a large

computational effort; therefore, such studies

32,39–43

have been limited

to nanoscopic droplets so far.

To extend the accessible length and time scales to those of hydro-

dynamics, Hoogerbrugge and Koelman

devised the dissipative parti-

cle dynamics (DPD) simulation technique. Compared with MD, the

computational effort is reduced by combining several molecules into

coarse-grained DPD particles. The challenge is then to appropriately

calibrate their interaction potential, e.g., to experimental data or MD

simulation results. Recently, a ﬁrst work

appeared corroborating the

universality of coalescence scaling laws by means of DPD simulations.

With traditional CFD methods, the Navier–Stokes (NS) equations

are solved directly,

21,24–26

which was done for coalescing liquid–

liquid systems

25,26

and liquid droplets surrounded by a gas phase.

21,24,33

Next to classical CFD simulations, several publications

27–29

used the

lattice Boltzmann (LB) method to study coalescence. In both NS- and

LB-based simulations, the interface requires speciﬁc treatment due to its

small width and surface tension, which is commonly not directly

included in the pressure tensor. Generally, there are two ways to model

the interface. A sharp interface deﬁnes a discontinuity of the thermody-

namic properties, while a diffuse interface imposes a continuous transi-

tion over a ﬁnite width. The ﬁrst additionally requires an interface

24,25

or boundary tracking method

to apply the surface tension. The phase-

ﬁeld method,

47–50

which is a diffuse interface method, is based on the

ideas of Cahn and Hillard,

who proposed a free energy functional that

depends on a scalar state variable and its spatial derivative. This scalar

property, referred to as order parameter in phase-ﬁeld models, can be,

e.g., the mole fraction

25–27

for binary liquid systems or the density

for

vapor–liquid systems.

The CFD method presented in this work is based on a phase-

ﬁeld model and the thermodynamic properties of the ﬂuid are

described by a Helmholtz energy functional that consists of a density

and a density gradient-dependent part. There are physically sound

expressions for the Helmholtz energy, i.e., equations of state (EoS), like

the perturbed chain statistical associating ﬂuid theory (PC-SAFT),

perturbed truncated and shifted (PeTS),

or the co-oriented ﬂuid

functional equation for electrostatic interactions (COFFEE),

to prop-

erly describe pure ﬂuids and mixtures. Combined with classical density

gradient theory,

51,55,56

the gradient term of phase-ﬁeld models can be

understood as a resistance to the formation of an interface and is

related to the second moment of the direct correlation function.

Using the PeTS EoS,

some of the present authors have shown that

vapor–liquid interfaces of pure ﬂuids

and mixtures

in static

58,59

FIG. 1. Snapshot of the initial conﬁguration of the R¼100 nm scenario consisting

of about 200 10

molecules, rendered by the open source visualization frame-

work MegaMol.

Molecules constituting the two droplets were colored green and

red, respectively, to follow their propagation. In order to obtain a clear view on the

interface, molecules of the vapor phase were rendered with low opacity. The magni-

ﬁed view on the contact region illustrates the small initial distance of a single mole-

cule diameter. For the droplets with R¼25 nm, the coalescence process was

rendered as a sequence. Multimedia view: https://doi.org/10.1063/5.0086131.1

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-2

Published under an exclusive license by AIP Publishing

dynamic

scenarios with

and without

58,59

the presence of solid walls

can be described almost solely based on this Helmholtz energy EoS.

The only adjustable parameter of the density gradient term was shown

to be state-independent.

For dynamic scenarios, additional models

for the transport properties are necessary.

Such phase-ﬁeld models

have the advantage that all thermodynamic properties are inherently

included. However, the length and time scales that need to be consid-

ered become small, demanding for a ﬁne discretization in space and

time, if no additional measures are taken. The present phase-ﬁeld

model does not require any a priori assumptions about the ﬂuid, its

phase equilibrium, or the investigated process. The only requirements

are a Helmholtz energy EoS that appropriately covers the metastable

and unstable regions by having a single van der Waals loop, as well as

expressions for the shear viscosity and the surface tension.

In a concerted effort, the phase-ﬁeld model and MD simulation

were applied in this work to coalescence, considering argon droplets in

three sizes with an initial radius ranging from 25 to 100 nm sur-

rounded by coexisting vapor.

II. MOLECULAR DYNAMICS

Resting on statistical mechanics, MD simulation is a physically

sound approach to study phenomena like droplet coalescence with an

atomistic level of detail. Except for the force ﬁeld describing the inter-

molecular interactions, practically no further assumptions have to be

made. Consequently, the results from this approach are well suited to

serve as a benchmark to evaluate macroscopic solutions that involve

much more severe assumptions. This opportunity was exploited here

to assess the phase-ﬁeld model described in Sec. III and ﬁndings from

the literature. Since the applicability of continuum assumptions

reaches its limit on the nanoscale, comparatively large droplets were

prepared to allow for a direct comparison, considering exactly the

same scenario in terms of geometry, ﬂuid, and state.

The truncated and shifted Lennard-Jones (LJTS) potential was

chosen for the intermolecular interactions. Since it has a small cutoff

radius of 2.5 molecule diameters rand no long range corrections have

to be taken into account, it is computationally rather cheap. Moreover,

the PeTS EoS

is available for the ﬂuid deﬁned by this potential that is

suitable for phase-ﬁeld modeling. The LJTS potential can be parameter-

ized to adequately represent the thermophysical properties of the noble

gases and methane.

In fact, the energy parameter =kB¼137:9K,the

size parameter r¼3:3916 ˚

A, and the molar mass M¼39:948 g mol1

were speciﬁed such that the present simulations mimic argon.

Three droplet sizes were considered, having an initial radius of

R¼25, 50, and 100 nm. To create an initial conﬁguration of the coa-

lescence process, the droplets were ﬁrst prepared in a vapor–liquid

equilibrium simulation at the desired temperature T¼110 K, which

corresponds to 73% of the critical temperature. The resulting satu-

rated liquid and vapor densities that slightly depend on the droplet

size were in good agreement with a correlation for exactly that inho-

mogeneous ﬂuid.

Subsequently, two copies of equilibrated droplets

were arranged along the xaxis of a cuboid volume with a distance of a

single molecule diameter, surrounded by vapor in equilibrium with

the liquid, cf. Fig. 1. In total, the generated systems consisted of about

3, 25, and 200 10

molecules. All MD simulations were conducted

withtheopensourcecodels1mardyn

that is well suited for mas-

sively parallel execution.

To capture the rapid dynamics of the coalescence process, the

spatially and temporally resolved density distribution was sampled

during simulation. For this purpose, the computational domain was

divided into equidistant cylinder shells around the xaxis, exploiting

the system’s rotational symmetry, and equidistant slices along the x

axis, both with a width of 3r=4. Within these sampling bins in the

form of cylinder shell elements, the mass density

q¼M

hNi

VðrÞ(1)

was averaged over a period of 5000 MD integration time steps with

Dt¼2 fs, thus amounting to 10 ps. V(r)isthebinvolumeand

depends on the shell radius r, i.e., distance to the xaxis; hNiis the

average number of molecules counted within the according bin; and

is the Avogadro constant. Because of the largely varying cylinder

shell volumes, the statistical quality of the data sampled near the xaxis

is lower than elsewhere.

The computational effort of the conducted MD simulations is

summarized in Table I.

III. PHASE-FIELD MODEL

To investigate droplet coalescence, the NS equations were solved

on the basis of a phase-ﬁeld model. The mass and momentum balan-

ces in the computational domain Bare

qþqrv¼0;(2)

v¼rr:(3)

Here, qðx;tÞand vðx;tÞdenote the density and velocity as a function

of position xand time t, the material time derivative is _

ðÞ¼ dðÞ=dt

and the symbol rstands for the nabla operator. The temperature in

the domain was held constant, and bulk viscous effects and volume

forces, like gravity, were neglected. The Cauchy stress tensor ris given

r¼2grsv1

3trðrsvÞ1



þfql

½

1jrqrq;(4)

with rsðÞ ¼ ðrðÞþðrðÞÞTÞ=2, trace trðÞ,identitymatrix1,

dyadic product , chemical potential l, and inﬂuence parameter j

that was adjusted to the surface tension.

This constitutive relation

includes a viscous part as well as pressure and surface tension resulting

from the phase-ﬁeld model. The Helmholtz energy density freads

fðq;rqÞ¼qaþj1

2jrqj2;(5)

TABLE I. Computational effort for the MD simulations in terms of core-h spent on

the national supercomputer HPE Apollo (Hawk) at the High Performance Computing

Center Stuttgart (HLRS), depending on the initial droplet radius R, volume of the

computational domain V(width height depth), number of molecules, number of

time steps N

, and the resulting physical time period t.

R(nm) V(nm

)N(–) N

(–) t(ns) core-h (h)

25 150 100 100 3:21065:010610.0 1:9105

50 300 200 200 2:510717:210634.3 2:7106

100 600 400 400 2:010810:210620.3 8:8106

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-3

Published under an exclusive license by AIP Publishing

where the Helmholtz energy awas speciﬁed by the PeTS EoS,

which

adequately covers metastable and unstable states and appropriately

describes the transition between the coexisting bulk phases. An expres-

sion for the shear viscosity gthat is based on MD simulation data

was taken from Ref. 60.

In conjunction with the NS equations, the chemical potential,

given by

l¼@fðq;rqÞ

@qr@fðq;rqÞ

@rq;(6)

was treated as a variable and solved at every time step.

Unlike MD, CFD simulations can be applied to problems exceed-

ing the nanoscale because the computational effort is substantially

lower. This is evident from the core-h spent on the supercomputers,

which is three orders of magnitude larger for MD simulations, cf.

Tables I and II. For the ﬁnite element (FE) implementation, the inter-

face width nonetheless represents a limiting factor for the mesh size

because it has to be resolved by at least four grid elements to ensure a

stable computation. To overcome this scale limitation, the interface

was artiﬁcially widened following Ref. 65.

For this purpose, a scaled formulation of the Helmholtz energy

density f

wasintroducedasfollows:

fscðq;rqÞ¼aðqaðqÞlsqþps

|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}

xðqÞ

Þþlsqpsþb1

2jrqj2:(7)

Here, p

denotes the vapor pressure and l

the chemical potential at

saturation. The two constants aand bscale the grand potential density

xðqÞand the density gradient-dependent part of the Helmholtz

energy density, respectively. They are deﬁned by the speciﬁed interface

width L

. When scaling the Helmholtz energy density, the thermody-

namic properties of the system, except for the interface width, must

not change. The ﬁrst derivative of the scaled Helmholtz energy with

respect to the density under saturation conditions reads

@fscðq;rqÞðÞ

@qq2fq0;q00g¼a@qaðqÞðÞ

@qq2fq0;q00gls

|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}

¼0

þls:(8)

This is identical with the original Helmholtz energy formulation (6)

and retains the consistency of the saturated densities. Furthermore, the

second derivative of the scaled Helmholtz energy with respect to the

density, given by

@2fscðq;rqÞðÞ

@q2¼a@2qaðqÞðÞ

@q2;(9)

is equal to the second derivative of the original Helmholtz energy for-

mulation, except for the constant a. Consequently, the spinodals,

where the second derivative of the Helmholtz energy density with

respect to the density is zero, remain unaltered. The effects of scaling

on the Helmholtz energy and grand potential density are visualized in

Fig. 2.Fora¼1, the curves are identical with the original Helmholtz

TABLE II. Computational effort for the CFD simulations in terms of core-h spent on

the cluster Elwetritsch at Regionales Hochschulrechenzentrum Kaiserslautern

(RHRK), depending on the initial droplet radius R, interface widening factor n, num-

ber of Q2P1 elements N

used for the discretization of the computational domain,

time step width Dt, number of time steps N

, and the resulting physical time period t.

The entry of core-h for R¼100 nm, n¼4 was omitted because that measurement

failed.

R(nm) n(–) Nel (–) Dt(ps) N

(–) t(ns) core-h (h)

25 1 336 144 1 400 0.4 9.7

25 2 168 72 1 400 0.4 4.2

25 4 84 36 1 400 0.4 2.4

25 1 336 144 10 1040 10.4 33.1

50 1 672 288 1 400 0.4 38.0

50 2 336 144 1 400 0.4 9.6

50 4 168 72 1 400 0.4 4.1

50 1 672 288 20 1270 25.4 398.5

100 1 1344 576 1 400 0.4 149.7

100 2 672 288 1 400 0.4 37.7

100 4 336 144 1 400 0.4 9.1

100 4 336 144 40 2040 81.6 

FIG. 2. Scaled Helmholtz energy density (top) and grand potential density (bottom)

over mass density qat a temperature of T¼110 K for different values of a. The

two minima of the grand potential density are located at the saturated densities q0

and q00.

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-4

Published under an exclusive license by AIP Publishing

energy formulation. Consistent saturated densities and an unaltered

unstable region width are retained for varying values of a.

The scaling of the interface width was carried out under the

assumption of a planar interface. However, nanoscopic droplets have

slightly deviating interface properties due to their strong curvature and

also the coexisting densities do not exactly correspond to those of a pla-

nar interface.

Since the purpose of interface widening is to deal with

larger scenarios, these differences were neglected. The scaled width L

of a planar interface that is perpendicular to the xaxis was deﬁned as

dxx¼xi¼q0q00

Lsc

with xi¼xðqxi ¼ðq0þq00Þ=2Þ;(10)

with x

being the position where the density takes the average of the

saturated liquid and vapor densities q0and q00.Withthisdeﬁnition,a

tangent is laid out on the density proﬁle at x

. The interface width is

then the distance of this tangent along the xaxis between its two inter-

section points with the saturated densities q0and q00, respectively.

In addition, the surface tension taken from Ref. 62 had to be

retained. The necessary requirement is

c¼X

A;(11)

where Ais the interface area and Xthe grand potential. For a planar

interface under equilibrium that is perpendicular to the xaxis,

X¼Aðþ1

1

axðqÞþ1

2bjrqj2



dx:(12)

Following the approach presented in Ref. 51 yields

axðqÞ¼1

2bjrqj2:(13)

This allows to rewrite the integral over xin Eq. (12) into an integral

over qusing

dx¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2axðqÞ

sdq;(14)

leading to

ðq0

q00 ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2abxðqÞ

pdq¼c:(15)

With Eq. (14), the deﬁnition of the interface width L

can be rewritten

dxx¼xi¼q0q00

Lsc ¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2axðq¼qxiÞ

s:(16)

Solving Eqs. (15) and (16) yields

a¼q0q00

2ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

xq¼qxi

ðÞ

pðq0

q00 ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

xðqÞ

pdq

Lsc

;(17)

b¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

xq¼qxi

ðÞ

q0q00

ðÞ

ðq0

q00 ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

xðqÞ

pdq

cLsc:(18)

The interface width can be widened by specifying a value for L

which allows to reduce the number of elements needed to discretize

the computational domain. Instead, or even additionally, an adaptive

grid could have been used to resolve the interface and reduce the num-

ber of elements even further, but this was beyond the scope of this

work.

The computational domain was formed utilizing a cylindrical

coordinate system with a radial distance r,angleu,andlengthx.Due

to the radial symmetry of the investigated coalescence scenario, all vec-

torial quantities in the direction of uand all derivatives with respect to

the angle uvanish, i.e.,

vu¼0;@v

@u¼0;@l

@u¼0;@q

@u¼0;(19)

thus, two-dimensional discretization sufﬁces. The boundary condi-

tions for the computational domain Bread

v¼0on@Bn@Br¼0;(20)

vr¼0on@Br¼0;(21)

rqn¼0on@B:(22)

Here, v

denotes the velocity perpendicular to the symmetry axis

@Br¼0and nthe outer normal to the boundary of the computational

domain @B,cf.Fig. 3. Due to these boundary conditions, the mass in

the computational domain remains constant. The temperature was

throughout speciﬁed to be constant at T¼110 K. The setup was ini-

tialized by a radial density proﬁle approximated by

FIG. 3. Bottom: setup of the computational domain for the coalescence scenario.

The magniﬁed view on the droplets’ contact region emphasizes the initial interface

width L0¼L, where Lis the interface width that was found under equilibrium condi-

tions. Top: detailed views with grids indicating the discretization used for the FE

scheme. In simulations with a widened interface Lsc ¼nL, for n¼2, half as many

but twice as large grid elements were used compared to the simulations without

employing the interface widening approach (n¼1).

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-5

Published under an exclusive license by AIP Publishing

qðrÞ¼q0q00

2tanh 2ðRrÞ



þ1



þq00;(23)

where the saturated liquid and vapor densities q0and q00 were taken

from Ref. 62. Following Ref. 58,avalueofL02:5rwas used for the

initial interface width. In simulations employing the interface widen-

ing approach, a scaled interface width Lsc ¼nLset to a multiple of the

true interface width Lunder equilibrium conditions was used instead.

In the present work, simulations with n¼1, 2, and 4 were performed.

An FE implementation was used as a numerical solution scheme;

details are provided in Ref. 65. For the coalescence scenario consider-

ing two droplets with an initial radius of R¼25 nm, a rectangle of

336 144 Q2P1 elements was used to discretize a domain of about

120 50 nm2. The larger scenarios were scaled proportionally in

terms of the domain’s extent and element count. Consequently, these

were doubled twice for R¼50 nm and R¼100 nm. For simulations

employing the interface widening approach, however, the number of

elements was scaled down by the factor nin both dimensions, cf. Fig.

3, yielding four (n¼2) and 16 times (n¼4) fewer elements and, thus,

reduced the computational effort roughly by a factor n2,assumma-

rized in Table II.Itshouldbenoted,however,thatthecomputational

effort for different time increments varied signiﬁcantly, depending on

theconvergenceoftheNewtonmethod.AsintheMDsimulations,

the initial distance between the two droplets was a single molecule

diameter r.Att¼0, the velocity v¼0 was set in the entire domain.

Figure 3 depicts the initialized domain for the R¼25 nm scenario.

The time discretization was chosen such that the numerical solution

efﬁciency was optimal, while ensuring convergence of the Newton

method. Because of the rapidly growing bridge in the early stages, a

time step of Dt¼1 ps was used up to time t¼400 ps. For the later

stages, time steps of Dt¼10;20, and 40 ps were used for the scenar-

ios with R¼25, 50, and 100 nm, respectively.

IV. RESULTS

For a qualitative comparison of the simulation results, it is useful

to superimpose the evolution of the coalescing droplets’ interface con-

tour obtained from MD and CFD simulations. Naively comparing

states with equal elapsed time shas the disadvantage that differences

in terms of droplet dynamics accumulate over time. Moreover, the

instant of contact s¼0 initiating the coalescence process cannot be

unambiguously determined; so, a perfectly synchronous time axis can-

not be guaranteed. Therefore, it is advantageous to compare states

with equal bridge radius r

,asshowninFig. 4.

The superimposed interface contours coincide almost perfectly

for all droplet sizes and states, ranging from rb=R0:3to1.3.This

indicates that the phase-ﬁeld model can reproduce the physical behav-

ior very well, since several conditions have to be satisﬁed for an accu-

rate evolution of the contour. First, the appropriate growth rate of the

bridge is crucial. Its growth, driven by the surface tension, induces a

stress in the nearby liquid, which can only be relaxed by setting the ini-

tially resting bulk liquid into motion. The bridge generates an inward

force that pulls the droplets together. The upper bound of this force is

given by the surface tension around the circumference of the bridge at

its minimum radius.

To better understand the sequence of cause and effect, Figs. 5

and 6show the evolution of the density and velocity ﬁelds resulting

from simulations with droplets of initial radius R¼25 nm and

R¼50 nm (the velocity ﬁeld was only sampled with CFD simulation).

To reveal density variations in the droplets, the color map in Fig. 5 was

scaled to a narrow range around the saturated liquid density. Since this

FIG. 4. Comparison of droplet interface

contours sampled with MD (symbols) and

CFD (solid lines), evaluated at time instan-

ces with equal bridge radius r

. Top: view

at the bridge at rb=R0:3, 0.5, and 0.7.

Bottom: global view at rb=R0:7, 1.1,

and 1.3. From left to right, columns show

results of simulations with initial droplet

radius R¼25, 50, and 100 nm. Results

for R¼100 nm were not sampled until

rb=R1:3 to limit the computational

effort.

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-6

Published under an exclusive license by AIP Publishing

scaling also emphasizes the density ﬂuctuations of the present MD

data, ten independent MD simulation runs were averaged for this rep-

resentation. Both MD and CFD show that early on, the density of the

bridge and its vicinity is signiﬁcantly below the saturation value. This

means that the initially rapidly growing bridge releases so much interfa-

cial energy within a very short time that the adjacent ﬂuid layers are

subjected to strong tensile forces that generate metastable states with a

negative pressure. These volume sections have the ambition to attain

physically stable states and hence demand for a supply of molecules

from the interior of the droplets to gain density. As a result, ﬂuid layers

further away from the bridge accelerate, cf. Fig. 6 (Multimedia view).

From the phase-ﬁeld model point of view, the sequence of cause

and effect can be traced by means of the stress tensor. At the begin-

ning, when the droplets rest, terms of Eq. (4) directly depending on

the density gradient or higher derivatives of density dominate the

behavior due to the strong curvature of the interface, where the drop-

lets are in contact. This situation gradually transforms to one where

the full static pressure contribution determines the evolution of the

system. The PeTS EoS is designed such that the Helmholtz energy can

also be safely evaluated in the metastable region, yielding the correct

negative pressure under such conditions. The pressure gradient is

relieved as adjacent liquid layers start to move, generating viscous

stress in the neighboring liquid layers at rest, and is ﬁnally relaxed by

the global motion of the droplets. The viscous part of the stress tensor

determines how fast the local motion can propagate to ﬁnally lead to

the mutual approach of the droplets. Therein, the shear viscosity

describes the resistance to deformation of a ﬂuid element. At low vis-

cosity, as in the present case with g0:144 m Pa s, the droplets

deform strongly before they start moving, which is opposed by inertial

forces. Consequently, the regions further away from the bridge remain

at their initial position for a certain period of time, as can be seen from

the “back of drop” positions x

. Therefore, the droplets must constrict

to meet the stress from the bridge region, cf. Fig. 7.

The velocity ﬁelds in Fig. 6 show that the global motion of the

droplets, including their back of drop positions x

, is not initiated

before the bridge has grown to approximately the initial radius R.

These observations are consistent with the experiments of Paulsen

et al.

with varying silicon oils. Droplets with a high viscosity almost

rigidly approached each other, whereas the ﬂuid motion was initially

limited to the vicinity of the growing bridge in the case of low

viscosity.

A quantitative comparison of the MD and CFD simulation

results over the entire coalescence process is shown in Fig. 8 by plot-

ting the bridge radius r

over time after contact s. A closer look at the

very early stages (s<100 ps) reveals that the ILV regime, in which

the bridge grows with the rate vbc=g, is preceded by another

FIG. 5. Density ﬁeld of droplets with

R¼25 nm at three instances of time

s¼200 ps (top), 400 ps (center), and

800 ps (bottom) sampled with MD (left)

and CFD (right). The color map was

scaled to a narrow range around the satu-

rated liquid density. Bridge regions are

highlighted by magniﬁed views.

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-7

Published under an exclusive license by AIP Publishing

regime that exhibits a substantially higher growth rate. That thermal

regime was identiﬁed by Perumanath et al.

who investigated coalesc-

ing water droplets with MD simulation. They proposed an order of

magnitude estimate for the thermal length lTlc,wherel

is the

width of the distribution function describing the probability of

occurrence of the site where the ﬁrst molecular bridge between the

droplets emerges, given by

lckBT



1=4

R1=2;(24)

with the Boltzmann constant k

. The probability of occurrence is rep-

resented by a normal distribution around the droplet’s line of

approach where the maximum probability meets the point of smallest

distance between the droplets at r¼0. This estimate agrees well with

the length l

that we determined at the point where the growth rate

transitions to viscous scaling, cf. Fig. 8. For the present argon droplets,

however, we found lTlcinstead of lT2lcas Perumanath et al.

did for water droplets. Interestingly, a higher growth rate was also

identiﬁed in the present CFD simulation results.

From the assumption lTlcand Eq. (24), it follows that

lT=R1=ﬃﬃﬃ

p; so, the thermal regime becomes negligible for macro-

scopic droplets. In the case of the present nanoscopic droplets, how-

ever, the thermal regime plays a dominant role in the early stages of

the coalescence process.

Hence, comparing the data shown in Fig. 6 with the asymptotic

behavior predicted by theory derived from investigations of macro-

scopic droplets, where the thermal regime is not considered, may lead

to confusion at ﬁrst sight. For the present droplets, characterized by

Oh 101,cf.Table III, the coalescence phase diagram proposed by

FIG. 6. Density and velocity ﬁelds of coa-

lescing droplets with an initial radius

R¼50 nm at four instances of time

s¼0:6;2:4;4:8;and 13:8 ns sampled

with CFD. Multimedia view: https://doi.org/

10.1063/5.0086131.2

FIG. 7. Superimposed contours of droplets with R¼100 nm in their initial state

(gray) and coalescing droplets with rb=R1 (yellow) sampled with MD. At the

later stage, global motion of the coalescing droplets just initiated such that the back

of drop positions x

at x=R62 started to move toward the center at x=R¼0.

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-8

Published under an exclusive license by AIP Publishing

Paulsen et al.

predicts a crossover from the ILV to the inertial regime

when the reduced bridge radius attains rb=R101. Based on the

condition Re 1, it can be estimated to occur when the bridge height

reaches the viscous length r2

b=Rlvwith

lv¼g2

qc ;(25)

which is equivalent to rbﬃﬃﬃﬃﬃﬃ

lvR

p.Comparingl

to ﬃﬃﬃﬃﬃﬃ

lvR

pin Table III

and Fig. 8 indicates that the thermal regime overrides almost half of

the ILV regime. This does not change the bridge radius where the iner-

tial regime is entered, but signiﬁcantly shortens the time to reach this

stage. Hence, before collapsing the present data by rescaling the bridge

radius and time axis with respect to the crossover radius rb!rb=rc

and crossover time s!s=scas proposed in Ref. 19,weshiftedthe

data by a certain amount of time. Thereby, it was pretended that the

thermal regime does not exist and the growth rate corresponds to that

of the viscous scaling right from the beginning, cf. Fig. 8 (bottom). The

collapse of the shifted data are shown in Fig. 9, exhibiting a rather

sharp crossover from the ILV to the inertial regime instead of a

smooth transition represented by the interpolation

r=rc¼21

s=scþ1

ﬃﬃﬃﬃﬃﬃﬃﬃﬃ

s=sc

1

;(26)

proposed in Ref. 19, which we assign to the nanoscopic size of the pre-

sent droplets.

FIG. 8. Temporal evolution of the bridge

radius r

sampled with MD (symbols) and

CFD (lines). The statistical uncertainties of

the MD data are typically 2% and thereby

within symbol size. From left to right, col-

umns show results of simulations with an

initial droplet radius R¼25, 50, and

100 nm. Shaded areas depict the thermal

regime l

(gray) and the ILV regime ﬃﬃﬃﬃﬃﬃ

lvR

(blue). Top: global view. CFD results

shown for R¼100 nm were obtained

from a simulation employing the interface

widening approach with n¼4 (short

dash), whereas that of R¼25 and 50 nm

were obtained from simulations without a

widened interface (solid). Horizontal lines

(dashed) mark the equilibrium state of the

merged droplets. Center: view at the early

stages. CFD results with and without wid-

ened interface are shown, corresponding

to n¼1 (solid), n¼2 (dashed), and

n¼4 (short dash). The dashed dotted

line depicts the growth rate vbc=gof

the bridge in the ILV regime, calculated

from the surface tension cand shear vis-

cosity ggiven in Table III. Bottom: data

shifted such that viscous scaling c=g

meets the origin.

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-9

Published under an exclusive license by AIP Publishing

Finally, the growth rate of the bridge vb¼drb=dsand the veloc-

ityofthebackofdroppositionsvbd ¼dxbd=dsevaluated from MD

and CFD simulation data were directly compared. Therefore, splines

were ﬁtted to the temporal evolution of the bridge radius r

and the

back of drop positions’ translation Dxbd ¼2Rjxbdjto obtain

smooth velocity proﬁles, cf. Fig. 10. For a better comparability of the

results for the three droplet sizes, units of length and time were

reduced by the initial radius Rand the inertial time scale s

.Thisillus-

tration emphasizes the agreement between the results of the two simu-

lation methods. Moreover, the similar dynamic behavior of all three

droplet sizes points to the applicability of the scaling laws. Evolved in

the thermal regime, the initially high bridge growth rate v

reduces

immediately to lower rates after passing the ILV and entering the iner-

tial regime at si1, were the back of drop positions start moving

>0). Once the maximum back of drop velocity v

has passed, the

bridge radius r

and the back of drop translation Dxbd reach their

maximum at the same time and subsequently the merged droplet

starts to oscillate around its equilibrium state, i.e., a resting droplet

with a radius of 21=3R. The amplitude and period of this weakly

damped oscillation increase for larger droplets due to higher inertial

forces.

Because of the rapid dynamics at the beginning, but also for a

generally better assignment of the states, a plot over the reduced bridge

radius is advantageous, cf. Fig. 11. It reveals that the deviations

between the MD and CFD data are mainly limited to the growth rate

at the early stages. This can be most clearly seen from the results of

R¼100 nm, which were obtained from a simulation with a widened

interface (n¼4), showing signiﬁcant deviations for v

up to

rb=R0:4. As can be concluded from Fig. 8, this is related to the fact

that with a more widened interface, the decrease in the high initial

growth rate caused by the strong curvature of the bridge and its menis-

cus becomes increasingly delayed. This effect seems to be related to the

different expression of the contact zone, which expands increasingly in

the radial direction with a more widened interface, cf. Fig. 3 (top). On

the time axis, however, the part that exhibits these deviations accounts

for only a short period due to the high velocity v

and therefore has

only a minor impact on the overall process, cf. Fig. 10. Compared to

MD, for CFD simulations the maximum bridge radius r

was reached

slightly more than hundred picoseconds earlier for R¼25 nm and

about two hundred picoseconds earlier for R¼50 nm, corresponding

to a relative deviation of about 2%.

Moreover, Fig. 11 shows that the back of drop positions of the

largest system start moving when the bridge radius reached the initial

radius rb=R¼1, whereas they do so increasingly earlier with smaller

droplet size. The oscillation of the velocities v

and v

after uniﬁcation

of the droplets manifests itself in this diagram as a spiral, leading to

the equilibrium point. The increasing size of the spiral with rising

droplet size, best seen for the velocity v

, again indicates the increas-

ingly strong expression of inertia-induced oscillation.

V. CONCLUSION

Coalescing argon droplets with a radius of 25, 50, and 100 nm

were studied with computational methods. Relatively large MD simu-

lations were carried out, containing up to 2 108molecules, to gener-

ate reference data. Moreover, a phase-ﬁeld model resting on a

Helmholtz energy EoS was devised that was at the core of CFD simula-

tions. Exactly the same scenarios in terms of geometry, ﬂuid, and state

were considered with these approaches. MD and CFD simulation

results showed an excellent agreement over the entire coalescence pro-

cess. Theoretical results on the asymptotic behavior of coalescence

process regimes were conﬁrmed. All considered scenarios crossed

from the inertially limited viscous regime over to the inertial regime

FIG. 9. Bridge radius r

over time after contact ssampled with MD (symbols) and

CFD (lines). The statistical uncertainties of the MD data are typically 2% and

thereby within symbol size. CFD results for R¼100 nm were obtained from a sim-

ulation employing the interface widening approach with n¼4, whereas those for

R¼25 and 50 nm were obtained from simulations without interface widening. The

data were shifted, cf. Fig. 8 and related text, and then collapsed by rescaling the

horizontal and vertical axes by the crossover time s

and the crossover radius r

each data set. The limiting behavior is proportional to s=sc(dotted line) at early

times and ðs=scÞ1=2at late times (dashed line).

A transition between these limit-

ing behaviors is described by Eq. (26) (solid black line).

TABLE III. Parameters characterizing the coalescing droplets, i.e., initial radius R, saturated liquid and vapor densities q0and q00, shear viscosity g, Ohnesorge number Oh,

thermal length l

, and an estimate for the bridge radius ﬃﬃﬃﬃﬃﬃ

lvR

pwhere the transition from ILV to inertial regime is expected.

R(nm) q0(kg m3)q00 (kg m3)g(mPa s) Oh (–) l

(nm) ﬃﬃﬃﬃﬃﬃ

lvR

p(nm)

25 1247.2 34.16 0.1451 0.3158 3.44 7.89

50 1245.6 33.56 0.1442 0.2222 4.87 11.11

100 1244.8 33.26 0.1438 0.1567 6.88 15.67

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-10

Published under an exclusive license by AIP Publishing

because of the low shear viscosity of argon. The present phase-ﬁeld

model captured the dynamics even at the initial stages of the coales-

cence process during the thermal regime that was only recently identi-

ﬁed by MD simulation. Moreover, a closer look at the liquid density

revealed that metastable states associated with negative pressure were

attained in the emerging liquid bridge between the coalescing droplets.

Thereby, it was demonstrated that the phase-ﬁeld model is even capa-

ble of adequately handling the onset of coalescence.

The vapor–liquid interface, where the surface tension acts, is the

crucial element for droplet coalescence and resides on the nanoscale in

general. Consequently, a nanoscopic discretization had to be speciﬁed

for the present CFD simulations resting on the phase-ﬁeld model.

With an interface widening approach, the phase-ﬁeld model was

transferred to coarser grids, while retaining the thermodynamic

properties including the surface tension. This allowed for a reduction

of the number of grid cells by a factor of 4 (n¼2) or even 16 (n¼4),

leading to almost the same simulation results with a considerably

reduced computational effort. Only in the case of extreme curvatures,

as they occur in the early stages of coalescence, increasing deviations

from the MD reference data were observed with rising interface wid-

ening. However, since strong curvatures are short-lived, they have a

minor effect on overall coalescence process.

Another advantage of the phase-ﬁeld model is its transferability

to other ﬂuids, including mixtures. The only requirements are an

appropriate Helmholtz energy EoS and expressions for the shear vis-

cosity and surface tension, where for the latter practically only a single

state point is needed to adjust the inﬂuence parameter. For future

work, it seems worthwhile to combine the present phase-ﬁeld model

FIG. 10. Temporal evolution of the bridge

radius r

, its growth rate v

, back of drop

positions’ translation Dxbd , and its velocity

. Proﬁles of r

and Dxbd were obtained

by ﬁtting splines to the simulation data

from MD (dashed line) and CFD (solid

line). Proﬁles of v

and v

were then

extracted from the ﬁrst derivatives. CFD

results for R¼100 nm were obtained

from a simulation employing the interface

widening approach with n¼4, whereas

those for R¼25 and 50 nm were

obtained from simulations without inter-

face widening. From left to right, columns

show results for an initial droplet radius

R¼25, 50, and 100 nm. Horizontal lines

mark the equilibrium state of the coa-

lesced droplet.

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-11

Published under an exclusive license by AIP Publishing

with an adaptive grid to further reduce the computational effort and to

address macroscopic droplets.

ACKNOWLEDGMENTS

This work was supported by the German Research Foundation

(DFG) through the Project SFB-TRR 75, Project No. 84292822

“Droplet Dynamics under Extreme Ambient Conditions” and SFB

926, Project No. 172116086 “Bauteiloberﬂ€

achen: Morphologie auf

der Mikroskala,” as well as the German Federal Ministry of

Education and Research (BMBF) under Grant No. 01IH16008

“TaLPas: Task-basierte Lastverteilung und Auto-Tuning in der

Partikelsimulation.” The MD simulations were performed on the

national supercomputer HPE Apollo (Hawk) at the High

Performance Computing Center Stuttgart (HLRS) and on the

cluster Cray CS500 (Noctua) at the Paderborn Center for Parallel

Computing (PC

).TheCFDsimulationswereperformedon

the cluster Elwetritsch at Regionales Hochschulrechenzentrum

Kaiserslautern (RHRK).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conﬂicts to disclose.

DATA AVAILABILITY

The data that support the ﬁndings of this study are openly avail-

able in the Data Repository of the University of Stuttgart (DaRUS) at

https://doi.org/10.18419/darus-2434,Ref.66 and https://doi.org/

10.18419/darus-2456,Ref.67.

REFERENCES

F. B. Wadsworth, J. Vasseur, E. W. Llewellin, J. Schauroth, K. J. Dobson, B.

Scheu, and D. B. Dingwell, “Sintering of viscous droplets under surface

tension,” Proc. R. Soc. A 472, 20150780 (2016).

J. Breitenbach, I. V. Roisman, and C. Tropea, “From drop impact physics to

spray cooling models: A critical review,” Exp. Fluids 59, 55 (2018).

J. Jin, C. Ooi, D. Dao, and N.-T. Nguyen, “Coalescence processes of droplets

and liquid marbles,” Micromachines 8, 336 (2017).

R. W. Hopper, “Plane Stokes ﬂow driven by capillarity on a free surface,”

J. Fluid Mech. 213, 349–375 (1990).

J. Eggers, J. R. Lister, and H. A. Stone, “Coalescence of liquid drops,” J. Fluid

Mech. 401, 293–310 (1999).

Y. D. Shikhmurzaev, “Coalescence and capillary breakup of liquid volumes,”

Phys. Fluids 12, 2386–2396 (2000).

A. B. Thompson and J. Billingham, “Inviscid coalescence in the presence of a

surrounding ﬂuid,” IMA J. Appl. Math. 77, 678–696 (2012).

H. J. Maris, “Analysis of the early stage of coalescence of helium drops,” Phys.

Rev. E 67, 066309 (2003).

L. Duchemin, J. Eggers, and C. Josserand, “Inviscid coalescence of drops,”

J. Fluid Mech. 487, 167–178 (2003).

FIG. 11. Growth rate v

, back of drop

positions’ translation Dxbd and its velocity

plotted over the bridge radius r

Proﬁles of r

(not shown) and Dxbd were

obtained by ﬁtting splines to the simulation

data from MD (dashed line) and CFD

(solid line). Proﬁles of v

and v

were

then extracted from the ﬁrst derivatives.

CFD results for R¼100 nm were

obtained from a simulation employing the

interface widening approach with n¼4,

whereas those for R¼25 and 50 nm

were obtained from simulations without

interface widening. From left to right, col-

umns show results for an initial droplet

radius R¼25, 50, and 100 nm. Horizontal

lines mark the equilibrium state of the coa-

lesced droplet.

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-12

Published under an exclusive license by AIP Publishing

T. G. Owe Berg, G. C. Fernish, and T. A. Gaukler, “The mechanism of coales-

cence of liquid drops,” J. Atmos. Sci. 20, 153–158 (1963).

T. M. Dreher, J. Glass, A. J. O’Connor, and G. W. Stevens, “Effect of rheology

on coalescence rates and emulsion stability,” AIChE J. 45, 1182–1190 (1999).

R. Ishiguro, F. Graner, E. Rolley, and S. Balibar, “Coalescence of crystalline

drops,” Phys. Rev. Lett. 93, 235301 (2004).

M. Wu, T. Cubaud, and C.-M. Ho, “Scaling law in liquid drop coalescence

driven by surface tension,” Phys. Fluids 16, L51–L54 (2004).

D. G. A. L. Aarts, H. N. W. Lekkerkerker, H. Guo, G. H. Wegdam, and D.

Bonn, “Hydrodynamics of droplet coalescence,” Phys. Rev. Lett. 95, 164503

(2005).

S. T. Thoroddsen, K. Takehara, and T. G. Etoh, “The coalescence speed of a

pendent and a sessile drop,” J. Fluid Mech. 527, 85–114 (2005).

K. Fezzaa and Y. Wang, “Ultrafast x-ray phase-contrast imaging of the initial

coalescence phase of two water droplets,” Phys. Rev. Lett. 100, 104501 (2008).

S. C. Case, “Coalescence of low-viscosity ﬂuids in air,” Phys. Rev. E 79, 026307

(2009).

S. T. Thoroddsen, B. Qian, T. G. Etoh, and K. Takehara, “The initial coales-

cence of miscible drops,” Phys. Fluids 19, 072110 (2007).

J. D. Paulsen, J. C. Burton, and S. R. Nagel, “Viscous to inertial crossover in liq-

uid drop coalescence,” Phys. Rev. Lett. 106, 114501 (2011).

M. M. Rahman, W. Lee, A. Iyer, and S. J. Williams, “Viscous resistance in drop

coalescence,” Phys. Fluids 31, 012104 (2019).

A. Menchaca-Rocha, A. Mart

ınez-D

avalos, R. N

nez, S. Popinet, and S. Zaleski,

“Coalescence of liquid drops by surface tension,” Phys. Rev. E 63, 046309

(2001).

J. D. Paulsen, J. C. Burton, S. R. Nagel, S. Appathurai, M. T. Harris, and O. A.

Basaran, “The inexorable resistance of inertia determines the initial regime of

drop coalescence,” Proc. Natl. Acad. Sci. U. S. A. 109, 6857–6861 (2012).

J. D. Paulsen, “Approach and coalescence of liquid drops in air,” Phys. Rev. E

88, 063010 (2013).

Y. Chen, C. Shen, and G. P. Peterson, “Hydrodynamics and morphologies of

droplet coalescence,” Ind. Eng. Chem. Res. 54, 9257–9262 (2015).

C. Shen, Y. Chen, C. Yu, and X. Liu, “Numerical study on the liquid-liquid

interface evolution during droplet coalescence,” Microgravity Sci. Technol. 32,

737–748 (2020).

P. Zimmermann, A. Mawbey, and T. Zeiner, “Calculation of droplet coales-

cence in binary liquid–liquid systems: An incompressible Cahn–Hilliard/

Navier–Stokes approach using the non-random two-liquid model,” J. Chem.

Eng. Data 65, 1083–1094 (2020).

L. Baroudi, M. Kawaji, and T. Lee, “Effects of initial conditions on the simula-

tion of inertial coalescence of two drops,” Comput. Math. Appl. 67, 282–289

(2014).

S. J. Lim, M. C. Choi, B. M. Weon, and B. Gim, “Lattice Boltzmann simulations

for water coalescence,” Appl. Phys. Lett. 111, 101602 (2017).

M. Gross, I. Steinbach, D. Raabe, and F. Varnik, “Viscous coalescence of drop-

lets: A lattice Boltzmann study,” Phys. Fluids 25, 052101 (2013).

J. E. Sprittles and Y. D. Shikhmurzaev, “Coalescence of liquid drops: Different

models versus experiment,” Phys. Fluids 24, 122105 (2012).

R. T. Eiswirth, H.-J. Bart, A. A. Ganguli, and E. Y. Kenig, “Experimental and

numerical investigation of binary coalescence: Liquid bridge building and inter-

nal ﬂow ﬁelds,” Phys. Fluids 24, 062108 (2012).

J.-C. Pothier and L. J. Lewis, “Molecular-dynamics study of the viscous to iner-

tial crossover in nanodroplet coalescence,” Phys. Rev. B 85, 115447 (2012).

H. Deka, G. Biswas, S. Chakraborty, and A. Dalal, “Coalescence dynamics of

unequal sized drops,” Phys. Fluids 31, 012105 (2019).

R. Stierle and J. Gross, “Hydrodynamic density functional theory for mixtures

from a variational principle and its application to droplet coalescence,”

J. Chem. Phys. 155, 134101 (2021).

A. H. Rajkotwala, H. Mirsandi, E. A. J. F. Peters, M. W. Baltussen, C. W. M.

van der Geld, J. G. M. Kuerten, and J. A. M. Kuipers, “Extension of local front

reconstruction method with controlled coalescence model,” Phys. Fluids 30,

022102 (2018).

J.-Y. Chen, P. Gao, Y.-T. Xia, E.-Q. Li, H.-R. Liu, and H. Ding, “Early stage of

delayed coalescence of soluble paired droplets: A numerical study,” Phys.

Fluids 33, 092005 (2021).

V. Cristini and Y.-C. Tan, “Theory and numerical simulation of droplet

dynamics in complex ﬂows—A review,” Lab Chip 4, 257–264 (2004).

J. Kamp, J. Villwock, and M. Kraume, “Drop coalescence in technical liquid/

liquid applications: A review on experimental techniques and modeling

approaches,” Rev. Chem. Eng. 33, 1–47 (2017).

S. Perumanath, M. K. Borg, M. V. Chubynsky, J. E. Sprittles, and J. M. Reese,

“Droplet coalescence is initiated by thermal motion,” Phys. Rev. Lett. 122,

104501 (2019).

J. Koplik, S. Pal, and J. R. Banavar, “Dynamics of nanoscale droplets,” Phys.

Rev. E 65, 021504 (2002).

L. Zhao and P. Choi, “Molecular dynamics simulation of the coalescence of

nanometer-sized water droplets in n-heptane,” J. Chem. Phys. 120, 1935–1942

(2004).

B.-B. Wang, X.-D. Wang, W.-M. Yan, and T.-H. Wang, “Molecular dynamics

simulations on coalescence and non-coalescence of conducting droplets,”

Langmuir 31, 7457–7462 (2015).

T. Li, Y. Xia, L. Zhang, X. Zhang, C. Fu, Y. Jiang, and H. Li, “Effects of stripy

surfaces with intervals on the coalescence dynamics of nano-droplets: Insights

from molecular dynamics simulations,” Appl. Surf. Sci. 481, 951–959 (2019).

S. Grottel, M. Krone, C. Muller, G. Reina, and T. Ertl, “MegaMol—A prototyp-

ing framework for particle-based visualization,” IEEE Trans. Visualization

Comput. Graph. 21, 201–214 (2015).

P. J. Hoogerbrugge and J. M. V. A. Koelman, “Simulating microscopic hydro-

dynamic phenomena with dissipative particle dynamics,” Europhys. Lett. 19,

155–160 (1992).

V. Akella and H. Gidituri, “Universal scaling laws in droplet coalescence: A dis-

sipative particle dynamics study,” Chem. Phys. Lett. 758, 137917 (2020).

N. Moelans, B. Blanpain, and P. Wollants, “An introduction to phase-ﬁeld

modeling of microstructure evolution,” CALPHAD 32, 268–294 (2008).

D. M. Anderson, G. B. McFadden, and A. A. Wheeler, “Diffuse-interface meth-

ods in ﬂuid mechanics,” Annu. Rev. Fluid Mech. 30, 139–165 (1998).

V. E. Badalassi, H. D. Ceniceros, and S. Banerjee, “Computation of multiphase

systems with phase ﬁeld models,” J. Comput. Phys. 190, 371–397 (2003).

D. Jacqmin, “Calculation of two-phase Navier–Stokes ﬂows using phase-ﬁeld

modeling,” J. Comput. Phys. 155, 96–127 (1999).

J. W. Cahn and J. E. Hilliard, “Free energy of a nonuniform system. I.

Interfacial free energy,” J. Chem. Phys. 28, 258–267 (1958).

J. Gross and G. Sadowski, “Perturbed-chain SAFT: An equation of state based

on a perturbation theory for chain molecules,” Ind. Eng. Chem. Res. 40,

1244–1260 (2001).

M. Heier, S. Stephan, J. Liu, W. G. Chapman, H. Hasse, and K. Langenbach,

“Equation of state for the Lennard-Jones truncated and shifted ﬂuid with a cut-

off radius of 2.5rbased on perturbation theory and its applications to interfa-

cial thermodynamics,” Mol. Phys. 116, 2083–2094 (2018).

K. Langenbach, “Co-oriented ﬂuid functional equation for electrostatic interac-

tions (COFFEE),” Chem. Eng. Sci. 174, 40–55 (2017).

H. Kahl and S. Enders, “Calculation of surface properties of pure ﬂuids using den-

sity gradient theory and SAFT-EOS,” Fluid Phase Equilib. 172, 27–42 (2000).

C. Miqueu, B. Mendiboure, A. Graciaa, and J. Lachaise, “Modelling of the sur-

face tension of pure components with the gradient theory of ﬂuid interfaces: A

simple and accurate expression for the inﬂuence parameters,” Fluid Phase

Equilib. 207, 225–246 (2003).

P. Rehner and J. Gross, “Predictive density gradient theory based on nonlocal

density functional theory,” Phys. Rev. E 98, 063312 (2018).

F. Diewald, M. Heier, M. Horsch, C. Kuhn, K. Langenbach, H. Hasse, and R.

M€

uller, “Three-dimensional phase ﬁeld modeling of inhomogeneous gas-liquid

systems using the PeTS equation of state,” J. Chem. Phys. 149, 064701 (2018).

S. Stephan, K. Langenbach, and H. Hasse, “Interfacial properties of binary

lennard-jones mixtures by molecular simulation and density gradient theory,”

J. Chem. Phys. 150, 174704 (2019).

F. Diewald, M. P. Lautenschlaeger, S. Stephan, K. Langenbach, C. Kuhn, S.

Seckler, H.-J. Bungartz, H. Hasse, and R. M€

uller, “Molecular dynamics and

phase ﬁeld simulations of droplets on surfaces with wettability gradient,”

Comput. Methods Appl. Mech. Eng. 361, 112773 (2020).

M. P. Lautenschlaeger and H. Hasse, “Transport properties of the Lennard-

Jones truncated and shifted ﬂuid from non-equilibrium molecular dynamics

simulations,” Fluid Phase Equilib. 482, 38–47 (2019).

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-13

Published under an exclusive license by AIP Publishing

J. Vrabec, G. K. Kedia, G. Fuchs, and H. Hasse, “Comprehensive study of the

vapour–liquid coexistence of the truncated and shifted Lennard–Jones ﬂuid

including planar and spherical interface properties,” Mol. Phys. 104,

1509–1527 (2006).

C. Niethammer, S. Becker, M. Bernreuther, M. Buchholz, W. Eckhardt, A.

Heinecke, S. Werth, H.-J. Bungartz, C. W. Glass, H. Hasse, J. Vrabec, and M.

Horsch, “ls1 mardyn: The massively parallel molecular dynamics code for large

systems,” J. Chem. Theory Comput. 10, 4455–4464 (2014).

N. Tchipev, S. Seckler, M. Heinen, J. Vrabec, F. Gratl, M. Horsch, M.

Bernreuther, C. W. Glass, C. Niethammer, N. Hammer, B. Krischok, M. Resch,

D. Kranzlm€

uller, H. Hasse, H.-J. Bungartz, and P. Neumann, “TweTriS:

Twenty trillion-atom simulation,” Int. J. High Perform. Comput. Appl. 33,

838–854 (2019).

F. Diewald, “Phase ﬁeld modeling of static and dynamic wetting,”

Forschungsbericht/Technische Universit€

at Kaiserslautern (Lehrstuhl Technische

Mechanik, 2020), Vol. 19.

M. Heinen and J. Vrabec, Coalescence of Two Nanoscopic Argon Droplets by

Molecular Dynamics Simulation (University of Stuttgart, 2022).

M. Hoffmann, F. Diewald, and K Langenbach, Coalescence of Two Nanoscopic

Argon Droplets by Phase-Field Modeling (University of Stuttgart, 2022).

Physics of Fluids ARTICLE scitation.org/journal/phf

Phys. Fluids 34, 042006 (2022); doi: 10.1063/5.0086131 34, 042006-14

Published under an exclusive license by AIP Publishing