1. Introduction
Droplet ripening, commonly known as Ostwald ripening, is a collisionless coarsening mechanism in the phase-separation process, characterised by the growth of larger droplets at the expense of smaller ones (Akira Reference Akira1982; Voorhees Reference Voorhees1985). Understanding its underlying mass transfer mechanisms is crucial for diverse applications, ranging from biological crystallisation (Joseph et al. Reference Joseph, Bernard, Jean, Anne, Daniel and Richard1996) and emulsion stability (Ariyaprakai & Dungan Reference Ariyaprakai and Dungan2010) to nanomaterial synthesis (Yec & Zeng Reference Yec and Zeng2014).
The first quantitative framework for this process was established by the classical Lifshitz–Slyozov–Wagner (LSW) theory (Lifshitz & Slyozov Reference Lifshitz and Slyozov1961; Wagner Reference Wagner1961). By assuming that droplets interact with a statistically uniform background, it distinguishes two limiting regimes: bulk diffusion-controlled ripening, where the growth rate scales inversely with the droplet radius, and interface-controlled ripening, where the growth rate is independent of the radius (Marqusee & Ross Reference Marqusee and Ross1984; Hardy & Voorhees Reference Hardy and Voorhees1988). Building upon this foundation, significant efforts have been made to investigate the process under broader governing mechanics (Burlakov Reference Burlakov2006; Tateno & Tanaka Reference Tateno and Tanaka2021) and bubble dynamics (Lemlich Reference Lemlich1978; Chieco & Durian Reference Chieco and Durian2023). To capture detailed kinetics beyond analytical limits, numerical approaches have also been extensively employed. Macroscopic methods like the lattice Boltzmann method efficiently reproduce large-scale evolution and droplet interactions (Shan & Chen Reference Shan and Chen1993; Chibbaro et al. Reference Chibbaro, Falcucci, Chiatti, Chen, Shan and Succi2008; Chen et al. Reference Chen, Kang, Mu, He and Tao2014, Reference Chen, Feng, Hu, Du and Wang2022; Qin et al. Reference Qin, Xie, Shan, Chen and Hu2025), while molecular dynamics provides crucial atomistic insights into droplet behaviours often masked in population-averaged predictions (Kraska Reference Kraska2008; Roy & Das Reference Roy and Das2012; Midya & Das Reference Midya and Das2017).
Despite these insights, current research faces several limitations. Physically, classical analytical theories including LSW and standard macroscopic simulations such as the lattice Boltzmann method typically rely on continuum description and often neglect the intrinsic microscopic structures of the liquid–vapour system, specifically the finite diffuse interface and the Knudsen layer. While acceptable for macroscopic fluids, the simplifications render continuum descriptions inadequate at the micro- and nanoscales where non-equilibrium interfacial effects govern mass transfer. Computationally, while molecular dynamics simulation resolves these features, its high cost prevents statistical analysis of late-stage ripening. To rigorously resolve these structures, molecular kinetic theory offers a powerful solution. Among available models, the Enskog–Vlasov equation is particularly effective, which naturally incorporates intermolecular attractions, enabling the capture of both the diffuse interface and the non-equilibrium Knudsen layer. This capability has been extensively validated in phase-change dynamics and nanodroplet growth (Frezzotti, Gibelli & Lorenzani Reference Frezzotti, Gibelli and Lorenzani2005; Zhang et al. Reference Zhang, Xu, Qiu, Wei and Wei2020; Takata, Matsumoto & Hattori Reference Takata, Matsumoto and Hattori2021). As the Enskog–Vlasov equation is still computationally expensive to solve, we utilise our recently developed simplified Enskog–Vlasov model to study droplet ripening in this work (Li et al. Reference Li, Su, Shan, Li, Gibelli and Zhang2024, Reference Li, Gibelli and Zhang2026). Validated for non-equilibrium evaporation flows, our model accurately resolves mesoscopic interfacial details with significantly enhanced computational efficiency, making it ideal for investigating the complete evolution of micro- and nanoscale ripening.
In this paper, we investigate droplet ripening in a two-dimensional single-component system. First, we focus on two-droplet systems to elucidate the fundamental ripening mechanism. Then, we derive a theoretical growth rate based on molecular flux imbalance and validate its linear scaling with the thermodynamic driving force. Crucially, we identify a late-stage deviation driven by the overlap of diffuse interfaces and establish a robust power-law scaling for the kinetic coefficient with temperature. Finally, simulations of three-droplet systems confirm the universality of the proposed theory, demonstrating that the ripening kinetics remain valid regardless of droplet number or spatial configuration.
The remainder of this paper is organised as follows. Section 2 introduces the molecular kinetic model for van der Waals fluids. Section 3 presents numerical results, derivation of the theoretical growth law and analysis of the nonlinear ripening regime and temperature dependence. Finally, § 4 summarises the main findings and discusses their implications.
2. Kinetic model of van der Waals fluids
Based on the mean-field theory, we consider a fluid composed of identical and spherical molecules of mass
$m$
and diameter
$\sigma$
, interacting by the Sutherland potential
$\phi (r)$
, resulting from the superposition of a repulsive hard sphere core and an attractive smooth tail, i.e.
\begin{align} \phi (r) = \left \{ \begin{aligned} &+\infty , \quad &r&\lt \sigma \\ & -\phi _{\sigma } \left (\frac {r}{\sigma }\right )^{-\gamma _{\sigma }}, &\quad r &\geqslant \sigma ,\end{aligned} \right . \end{align}
where
$\phi _{\sigma }$
and
$\gamma _{\sigma }$
are two positive constants related to the depth of the potential well and the extent of the soft interaction, respectively, and
$r=||\boldsymbol{x}_{1}-\boldsymbol{x}||$
is the distance between the interacting atoms at
$\boldsymbol{x}_{1}$
and
$\boldsymbol{x}$
. In this work, the fluid is described statistically in terms of molecular distribution functions. As droplet ripening is being investigated in a two-dimensional system here, two reduced distribution functions,
$G_{1}(\boldsymbol{x}, \boldsymbol{\xi }, t)$
and
$G_{2}(\boldsymbol{x}, \boldsymbol{\xi }, t)$
, are introduced to save computational cost. These functions represent the number of atoms at time
$t$
in the elementary volume of the single-particle phase space around the position
$\boldsymbol{x} = (x, y)$
and the velocity
$\boldsymbol{\xi } = (\xi _{x}, \xi _{y})$
, i.e.
\begin{equation} G_{1} = \int f {\rm d} \xi _{z}, \quad G_{2} = \int \xi ^{2}_{z} f {\rm d} \xi _{z}, \end{equation}
where
$f$
is the full molecular distribution function in the three-dimensional velocity space and
$G_{1}$
satisfies
\begin{equation} \frac {\partial G_{1}}{\partial t} + \boldsymbol{\xi } \boldsymbol{\cdot }\frac {\partial G_{1}}{\partial \boldsymbol{x}} + \frac {\boldsymbol{F}}{m} \boldsymbol{\cdot }\frac {\partial G_{1}}{\partial \boldsymbol{\xi }} = S_{a}^{G_{1}}+S_{b}^{G_{1}}. \end{equation}
The
$S_{a}^{G_{1}}$
contains the Shahkov-like terms:
\begin{equation} S_{a}^{G_{1}} = \frac {G_{1}^{eq}-G_{1}}{\tau }+\frac {G_{1}^{eq}}{\tau }\frac {2m(1-Pr)\boldsymbol{q^{K}} \boldsymbol{\cdot }\boldsymbol{C}}{5n(k_{B}T)^{2}}\left (\frac {mC^{2}}{2k_{B}T}-2\right )\!, \end{equation}
while
$S_{b}^{G_{1}}$
includes all the derivative terms of the macroscopic quantities originating from the Enskog collision term:
\begin{equation} \begin{aligned} S_{b}^{G_{1}} = -\rho b\chi G_{1}^{eq}&\left \{\boldsymbol{C} \boldsymbol{\cdot }\left [\frac {2}{n}\frac {\partial n}{\partial \boldsymbol{x}}+\frac {1}{T}\frac {\partial T}{\partial \boldsymbol{x}}\left (\frac {3mC^{2}}{10k_{B}T}-\frac {1}{5}\right )\right ]+ \frac {2m}{5k_{B}T}\boldsymbol{CC}:\frac {\partial }{\partial \boldsymbol{x}}\boldsymbol{u} \right . \\[3pt] & \left . +\frac {2}{5}\left (\frac {mC^{2}}{2k_{B}T}-2\right )\left (\frac {\partial }{\partial \boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{u}\right ) \right \}-\rho b \boldsymbol{C} G_{1}^{eq}\boldsymbol{\cdot }\frac {\partial \chi }{\partial \boldsymbol{x}}\\[3pt] &+\frac {\partial }{\partial \boldsymbol{x}} \boldsymbol{\cdot }\left [G_{1}^{eq} \frac {\varpi }{nk_{B}T}\left (\frac {\partial }{\partial \boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{u}\right ) \boldsymbol{C} \left (\frac {mC^{2}}{2k_{B}T}-1\right )\right ]\!. \end{aligned} \end{equation}
Here,
$G_{1}^{eq}$
is the local Maxwellian distribution function:
\begin{equation} G_{1}^{eq} = \frac {nm}{2\pi k_{B}T} \textrm {exp} \left (-\frac {m C^{2}}{2k_{B}T}\right )\!. \end{equation}
Similarly, the reduced distribution function
$G_{2}$
satisfies
\begin{equation} \frac {\partial G_{2}}{\partial t} + \boldsymbol{\xi } \boldsymbol{\cdot }\frac {\partial G_{2}}{\partial \boldsymbol{x}} + \frac {\boldsymbol{F}}{m} \boldsymbol{\cdot }\frac {\partial G_{2}}{\partial \boldsymbol{\xi }} = S_{a}^{G_{2}}+S_{b}^{G_{2}}. \end{equation}
The terms on the right-hand side of (2.4a ) can be obtained by
\begin{equation} S_{a}^{G_{2}} = \frac {G_{2}^{eq}-G_{2}}{\tau }+\frac {G_{2}^{eq}}{\tau }\frac {2m(1-Pr)\boldsymbol{q^{K}} \boldsymbol{\cdot }\boldsymbol{C}}{5n(k_{B}T)^{2}}\left (\frac {mC^{2}}{2k_{B}T}-1\right )\!, \end{equation}
\begin{equation} \begin{aligned} S_{b}^{G_{2}} = -\rho b\chi G_{2}^{eq}&\left \{\boldsymbol{C} \boldsymbol{\cdot }\left [\frac {2}{n}\frac {\partial n}{\partial \boldsymbol{x}}+\frac {1}{T}\frac {\partial T}{\partial \boldsymbol{x}}\left (\frac {3mC^{2}}{10k_{B}T}+\frac {2}{5}\right )\right ]+ \frac {2m}{5k_{B}T}\boldsymbol{CC}:\frac {\partial }{\partial \boldsymbol{x}}\boldsymbol{u} \right . \\[5pt] & \left . +\frac {2}{5}\left (\frac {mC^{2}}{2k_{B}T}-1\right )\left (\frac {\partial }{\partial \boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{u}\right ) \right \}-\rho b \boldsymbol{C} G_{2}^{eq}\boldsymbol{\cdot }\frac {\partial \chi }{\partial \boldsymbol{x}} \\[5pt] &\frac {\partial }{\partial \boldsymbol{x}} \boldsymbol{\cdot }\left [G_{2}^{eq} \frac {\varpi }{nk_{B}T}\left (\frac {\partial }{\partial \boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{u}\right ) \frac {m\boldsymbol{C}C^{2}}{2k_{B}T}\right ]\!, \end{aligned} \end{equation}
\begin{equation} G_{2}^{eq} = \frac {k_{B}T}{m}G_{1}^{eq}. \end{equation}
In the above equations,
$\chi$
is the pair correlation function of van der Waals fluid at equilibrium evaluated at the contact point, which accounts for correlations between colliding particles. The pair correlation function is related to the local density
$\rho$
and can be calculated by (Chapman & Cowling Reference Chapman and Cowling1990)
\begin{equation} \chi (\rho ) = \frac {1-0.125b\rho }{(1-0.25b\rho )^3}, \end{equation}
where
$b = 2\pi \sigma ^{3}/3m$
, which is related to the reduced density
$\eta = b\rho /4$
,
$n$
is the number density,
$\boldsymbol{u}$
is the bulk velocity,
$T$
is the temperature and
$\boldsymbol{q}^{K}$
is the kinetic heat flux. These macroscopic quantities, as well as the kinetic stress tensor
$\boldsymbol{p}^{K}$
, are determined by the velocity moments:
\begin{align} [n, \boldsymbol{u}, T, \boldsymbol{q}^{K}, \boldsymbol{p}^{K}] = \int \! \left [ G_{1}, \frac {\boldsymbol{\xi }G_{1}}{n}, \frac {m}{2nc_{v}}(C^{2}G_{1}+G_{2}), \frac {m\boldsymbol{C}}{2}(C^{2}G_{1}+G_{2}), m\boldsymbol{C}\boldsymbol{C}G_{1}\right ] \textrm {d}\boldsymbol{\xi }, \end{align}
where
$c_{v} = 3k_{B}/2m$
is the specific heat capacity at constant volume,
$k_B$
is the Boltzmann constant and
$\boldsymbol{C} = \boldsymbol{\xi }-\boldsymbol{u}$
is the peculiar velocity. The relaxation time
$\tau$
and Prandtl number
$Pr$
can be obtained from
\begin{equation} \tau = \frac {\mu }{nk_{B}T}\frac {1}{1+\frac {2}{5}\rho b\chi }, \end{equation}
\begin{equation} Pr = \frac {5}{2}\frac {k_{B}}{m}\frac {\mu }{\kappa }\frac {1+\frac {3}{5}\rho b \chi }{1+\frac {2}{5}\rho b \chi }. \end{equation}
The macroscopic transport equations and the associated transport coefficients can be derived using the Chapman–Enskog expansion, details being given in our previous works (Shan et al. Reference Shan, Su, Gibelli and Zhang2023; Su et al. Reference Su, Gibelli, Li, Borg and Zhang2023). To account for real-gas effects, the transport coefficients for dense gases, shear viscosity
$\mu$
, thermal conductivity
$\kappa$
and bulk viscosity
$\varpi$
, are evaluated using established empirical formulations. Specifically,
$\mu$
and
$\kappa$
incorporate collision integrals with density-dependent corrections following Chung, Lee & Starling (Reference Chung, Lee and Starling1984) and Chapman & Cowling (Reference Chapman and Cowling1990), while
$\varpi$
is evaluated using the approach for van der Waals fluids proposed by Hoover et al. (Reference Hoover, Ladd, Hickman and Holian1980). Explicit expressions can be found in Shan et al. (Reference Shan, Su, Gibelli and Zhang2023).
The final term to be defined is the mean-field force
$\boldsymbol{F}$
in two-dimensional form. From our previous work (Barbante, Frezzotti & Gibelli Reference Barbante, Frezzotti and Gibelli2015), the force term can be calculated by
\begin{equation} \begin{aligned} \boldsymbol{F}(x,y,t) = \int &[(x_{1}-x)\boldsymbol{i} + (y_{1}-y)\boldsymbol{j}]n(x_{1},y_{1},t) \\ & K(x_{1}-x, y_{1}-y)\,{\rm d}x_{1}\,{\rm d}y_{1}, \end{aligned} \end{equation}
where
$\boldsymbol{i}$
and
$\boldsymbol{j}$
are the unit vectors in the
$x$
and
$y$
directions. The kernel
$K$
depends on the attractive tail of the chosen potential, which can be calculated as
\begin{equation} K(x_{1}-x, y_{1}-y) = \int ^{+\infty }_{-\infty } \frac {1}{r}\frac {\partial \phi (r)}{\partial r}{\rm d}z, \end{equation}
where
$r = \sqrt {(x_{1}-x)^2+(y_{1}-y)^2+z^2}$
.
3. Results and discussion
3.1. Simulation set-up
The simulations described in the following sections use a rectangular computational domain with periodic boundary conditions containing two or three circular droplets surrounded by their vapour. We investigate two configurations: a two-droplet system to elucidate fundamental dynamics and validate the linear growth law, followed by a three-droplet system to confirm the generality of the proposed scaling and ensure that ripening kinetics remain independent of droplet number and spatial arrangement.
All the simulations use dimensionless quantities, with the molecular diameter as the reference length
$\ell _0$
, the ambient temperature as the reference temperature
$T_0$
and the most probable molecular velocity at
$T_0$
as the reference velocity
$u_0 = \sqrt {2k_B T_0 / m}$
. For argon at ambient temperature of
$273\,\rm K$
,
$\ell _0 \approx 0.3\,\rm nm$
,
$u_0 \approx 350\rm\,m\,s^{-1}$
. Accordingly, time
$t$
is normalised by
$\ell _0/u_0$
, number density
$n$
by
$n_0=1/\ell _0^3$
, the mean force field by
$u_0^2/\ell _0$
and the distribution function by
$n_0/u_0^3$
. The dimensionless quantities used in the following text are distinguished by a tilde. The molecular kinetic model is solved using the discrete velocity method (Wang et al. Reference Wang, Wu, Ho, Li, Li and Zhang2020), with a uniform mesh for the spatial discretisation, and Gaussian–Hermite quadrature nodes for the velocity discretisation. The derivatives are approximated by second-order finite differences. An explicit first-order Euler scheme is used for time integration. Based on the grid-independence test, all the simulations use a spatial grid size of
$\ell _0 / 9$
and
$8$
velocity nodes in each direction. To reduce computational cost while maintaining numerical accuracy, the cutoff radius for evaluating the mean-field force was set to five molecular diameters through numerical analysis.
3.2. Case 1: two-droplet system
In this section we describe the evolution of a two-droplet system from initial equilibration to ripening. The system consists of two droplets of different radii, positioned symmetrically to ensure they maintain their initial positions under periodic boundary conditions. Figure 1 shows the evolution at
$\tilde {T} = 0.6$
, which is below the dimensionless critical temperature
$T_{c}/T_{0} = 0.7546$
for the chosen intermolecular potential, ensuring phase separation. The initial diameter of the larger droplet is 20, while that of the smaller droplet is 16. The droplets are arranged with a centre-to-centre separation distance of 34. During the initial stage, both droplets relax towards quasi-steady states. Non-equilibrium effects in the vapour near the interface manifest as deviations from the local Maxwellian distribution function, which is not shown here for brevity’s sake. At this stage, the thickness of the diffuse liquid–vapour interface is measured to be approximately four molecular diameters. Subsequently, the larger droplet grows at the expense of the smaller one without forming a liquid bridge. This process accelerates as the smaller droplet shrinks, eventually leaving a single stationary droplet. In our simulations, the magnitude of maximum spurious velocities is ensured to be negligibly smaller than the ripening velocity, evaluated at the midpoint between the two droplets along
$y=H/2$
. Therefore, the influence of discretisation-induced spurious velocities on the reported ripening dynamics can be neglected.
(a–f) Time-lapse images depicting the evolution of morphology and velocity during the droplet ripening process at
$\tilde {T} = 0.6$
. The time of each image is displayed in the top left corner, and the velocity vectors are represented by white arrows. The black arrow at
$\tilde {t} = 50$
indicates the reference velocity vector, corresponding to a dimensionless magnitude of 0.015.
Figure 1. Long description
The image sequence illustrates the evolution of droplet morphology and velocity during the ripening process. Each frame displays a different time point, with the specific time indicated in the top left corner. White arrows represent velocity vectors, showing the direction and magnitude of movement. A black reference arrow indicates a dimensionless velocity magnitude of 0.015. The droplets change in size over time, demonstrating the growth of larger droplets at the expense of smaller ones.
Distinct from classical multicomponent ripening driven by concentration gradients, this single-component process is governed by pressure differences described by the Gibbs–Thomson relation (also referred to as the Kelvin equation), which links the local equilibrium vapour pressure to the droplet radius,
$p_{v1} = p_{eq} \exp { ( {B}/{R} )}$
, where
$p_{v1}$
is the saturation vapour pressure at the interface of a droplet with radius
$R$
and
$p_{eq}$
is the equilibrium pressure at a planar interface. The constant
$B$
represents the competition between surface tension and thermal agitation. This equation, also called the Ostwald–Freundlich equation, clearly shows that smaller droplets have higher vapour pressures at equilibrium. Consequently, in a system with two droplets of different sizes, a pressure gradient forms between them, with the smaller droplet being surrounded by a higher local vapour pressure.
Our simulation results provide direct quantitative evidence for the ripening process driven by the curvature-induced pressure difference. Figure 2 depicts the number density profile along the centreline of the computational domain at
$\tilde {t} = 1800$
. The profile reveals a distinct density gradient between the two droplets, establishing the direction of mass transfer from the smaller to the larger droplet. Crucially, the ambient vapour density is observed to be very low. This state confirms that mass transfer occurs through a collisionless vapour phase rather than via liquid bridge formation or droplet coalescence. Not only do these results substantiate our theoretical understanding of the pressure-driven ripening mechanism in single-component systems but also provide a clear visualisation of the transport dynamics.
(a) Number density profile along the x direction at
$y = H/2$
. Enlarged views of the number density (b) and pressure (c) in the vapour region, as indicated by the red dashed box in (a), are also shown.
Figure 2. Long description
The image contains three line graphs. The first graph (a) shows the number density profile along the x direction, highlighting regions of large and small droplets and a vapour region. The second graph (b) provides an enlarged view of the number density in the vapour region, as indicated by a red dashed box in the first graph. The third graph (c) shows an enlarged view of the pressure in the same vapour region. The x-axis represents the position along the x direction, while the y-axes represent number density and pressure respectively. The graphs illustrate the changes in number density and pressure within the vapour region during the droplet ripening process.
We now proceed to a detailed quantitative analysis of the ripening kinetics, focusing on the evolution of individual droplet radii and the governing mechanism driven by interfacial exchanges. Consider a general system containing
$N$
droplets. In the regime where interfacial resistance dominates bulk transport, the mass transfer is fundamentally governed by the molecular exchange kinetics at the liquid–vapour interface. The net mass flux
$J_i$
at the interface of droplet
$i$
arises from the intrinsic imbalance between the evaporation of molecules from the liquid surface and the condensation of vapour molecules from the surroundings. According to gas kinetic theory, this net flux is linearly proportional to the difference between the local saturation number density and the ambient vapour number density (Ytrehus & Østmo Reference Ytrehus and Østmo1996):
\begin{equation} J_{i}=K_{\textit{kin}}(n_{i}-n_{vap}), \end{equation}
where
$K_{\textit{kin}}$
is the kinetic coefficient. Here,
$n_{i}$
denotes the saturation vapour number density at the surface of droplet
$i$
and
$n_{vap}$
represents the ambient vapour number density determined by the collective behaviour of the droplets. Considering mass conservation for the closed system, the total rate of mass change vanishes, i.e.
\begin{equation} \sum _{i=1}^{N} \frac {\textrm {d}M_i}{\textrm {d}t} = 0. \end{equation}
For a two-dimensional system, the mass of a droplet is proportional to its area (
$M_i \propto R_i^{2}$
), and the mass change rate is related to the flux by
$\mathrm{d}M_i/\mathrm{d}t \propto -J_i \boldsymbol{\cdot }R_i$
. Combining (3.1) and (3.2), we obtain the expression for
$n_{vap}$
:
\begin{equation} n_{vap}=\frac {\sum R_{i}n_{i}}{\sum R_{i}}. \end{equation}
Then, we define the critical radius
$R_c$
such that the ambient density
$n_{vap}$
corresponds to the saturation density of a droplet with radius
$R_{c}$
. By applying the Gibbs–Thomson equation and performing a first-order Taylor expansion, we derive the rate of radius change for a single droplet:
\begin{equation} \frac {\textrm {d}R_{i}}{\textrm {d}t}=K\!\left (\frac {1}{R_{c}}-\frac {1}{R_{i}}\right )\!, \end{equation}
where
$R_{c}$
is found to be the arithmetic mean of the droplet radii:
\begin{equation} R_{c} = \frac {\sum R_{i}}{N}, \end{equation}
and
$K \propto K_{\textit{kin}} \gamma n_{\textit{sat}}/{n_l^2}$
, where
$\gamma$
is surface tension,
$n_{\textit{sat}}$
is equilibrium vapour density and
$n_l$
is liquid density. From (3.4), it is evident that droplets smaller than the critical radius
$R_c$
(i.e. the average radius) will shrink, while larger droplets will grow.
(a) Relationship between the radius change rate and the thermodynamic driving force (
$1/\tilde {R}_c - 1/\tilde {R}$
). The red circles represent simulation data, while the dashed line indicates the best linear fit crossing the origin point. The vertical line at
$\tilde {t} = 4300$
marks the transition from the linear kinetic regime to the nonlinear deviation regime. (b) Cross-sectional number density profile of the small droplet along the centreline. (c) Relationship between the change of radius and
$(1/\tilde {R}_{c}^{\textit{diff}} - 1/\tilde {R})/\tilde {R}$
under the diffusion-controlled assumption.
Figure 3. Long description
The image contains three graphs. The first graph shows the relationship between the radius change rate and the thermodynamic driving force. Red circles represent simulation data, and a dashed line indicates the best linear fit crossing the origin point. A vertical line marks the transition from the linear kinetic regime to the nonlinear deviation regime. The second graph displays the cross-sectional number density profile of a small droplet along the centerline, with different colored lines representing various time points. The third graph illustrates the relationship between the change of radius and the diffusion-controlled assumption, with red circles representing simulation data and dashed lines indicating linear fits.
(a) Dependence of the droplet radius change rate on the thermodynamic driving force for temperatures
$\tilde {T} = 0.55,\ 0.6$
and
$0.65$
. The dashed lines indicate the best linear fits to the simulation data. (b) Variation of the kinetic coefficient
$K$
with temperature
$\tilde {T}$
. The dashed line represents the best power-law fit, given by
$K = 0.525\tilde {T}^{6.52}$
.
Figure 4. Long description
The image contains two graphs. The first graph on the left shows the dependence of the droplet radius change rate on the thermodynamic driving force for different temperatures. The x-axis represents the inverse of the normalized droplet radius, and the y-axis represents the rate of change of the droplet radius. Three datasets are plotted with different symbols and colors: red squares for a temperature of 0.55, blue circles for 0.6, and yellow triangles for 0.65. Dashed lines indicate the best linear fits to the simulation data. The second graph on the right shows the variation of the kinetic coefficient with temperature. The x-axis represents the temperature, and the y-axis represents the kinetic coefficient. Red circles represent simulation results, and a dashed line represents the best power-law fit, given by the equation K equals 0.525 times T to the power of 6.52. The graphs illustrate how the kinetic coefficient increases with temperature and how the droplet radius change rate varies with the thermodynamic driving force at different temperatures.
To validate the theoretical model, we evaluated the relationship between the radius change rate and the thermodynamic driving force, defined as
$1/R_{c}-1/R$
, in a two-droplet system. For a system governed by the kinetics of molecular exchange, this relationship is predicted to be linear with a constant slope determined by the fluid properties. As shown in figure 3(a), the simulation results for
$\tilde {T} = 0.6$
exhibit a remarkable linear dependency up to
$\tilde {t} = 4300$
that passes through the origin, confirming that the ripening dynamics is indeed driven by the imbalance of molecular fluxes at the liquid–vapour boundary. Beyond this point, a distinct deviation from linearity is observed, as shown in figure 3(b). By monitoring the cross-sectional density profiles, we identified that this transition coincides with disappearance of bulk liquid of the shrinking droplet, where the diffuse interface is no longer distinguished from the bulk phase, rendering the theoretical assumption of a constant liquid number density invalid and indicating that the droplet is no longer a distinct thermodynamic phase. To quantify this breakdown, we introduce the dimensionless number
$\zeta = R/\delta$
, where
$R$
is the droplet radius and
$\delta$
is the diffuse interface thickness. The onset of interface overlap defines the critical threshold
$\zeta = 0.5$
. It should be noted that the above analysis differs from the classical LSW theory, in which ripening is assumed to be diffusion-controlled. Within that framework, the mass flux is governed by Fick’s law, leading to the following expressions for the droplet growth rate and the critical radius:
\begin{equation} \frac {{\rm d}R_i}{{\rm d}t} \propto \frac {1}{R_i} \left ( \frac {1}{R_c^{\textit{diff}}} - \frac {1}{R_i} \right )\!, \quad R_c^{\textit{diff}} = \frac {N}{\sum _{i=1}^{N} \frac {1}{R_i}}. \end{equation}
To examine whether this mechanism applies to the present system, we further analysed the simulation data by plotting
${\rm d}R/{\rm d}t$
and
$({1}/{R}) (1/R_c^{\textit{diff}} - 1/R)$
. As shown in figure 3(c), the resulting data do not pass through the origin, indicating a clear deviation from the behaviour predicted by the classical diffusion-controlled LSW theory. This observation suggests that the ripening process at the present scale is not governed by diffusion.
To complement our understanding of the ripening dynamics, we investigated the temperature dependence of the kinetic coefficient
$K$
, i.e. the slope of the linear growth regime. Figure 4(a) presents the simulation results for
$\tilde {T} = 0.55,\ 0.6$
and
$0.65$
. The linear relationship remains robust across these temperatures, with
$K$
increasing significantly as temperature rises. Theoretically,
$K$
is a composite parameter determined by multiple temperature-dependent fluid properties. To capture this complex dependence in our simulations, its value is obtained via empirical data fitting. As shown in figure 4(b), the variation of
$K$
is well described by a power-law correlation,
$K = 0.525\tilde {T}^{6.52}$
, over the investigated temperature range. This scaling law provides a practical means to estimate the ripening rate for the temperature range; the same functional form can be applied to other systems, with the exponent and prefactor adjusted to the specific fluid.
3.3. Case 2: three-droplet system
In this section, we extend our investigation to a three-droplet system to examine the proposed theory, with a specific focus on the determination of the critical radius
$R_{c}$
. The computational domain has dimensions of
$60 \times 54$
and contains three droplets with radii
$\tilde {R}_{1} = 9$
,
$\tilde {R}_{2} = 8$
and
$\tilde {R}_{3} = 7$
. The droplets are initially distributed arbitrarily within the domain at
$\tilde {T} = 0.6$
. The initial centre coordinates for the large, intermediate and small droplets are (16, 14), (46, 14) and (37, 42), respectively. According to our theoretical derivation for a two-dimensional system,
$R_c$
is determined by the arithmetic mean of the droplet radii. Consequently, for this initial configuration, the critical radius is exactly
$\tilde {R}_{c} = 8$
. Figure 5(a) shows the evolution of the droplet radii over time. Initially, the critical radius
$R_{c}$
of the system decreases gradually, and the slightly larger middle droplet grows slowly. At
$\tilde {t} \approx 4500$
, the small droplet fully dissolves, causing
$R_{c}$
to rise sharply above the middle droplet’s radius. This reverses the growth driving force for the middle droplet, causing its rate of change (
${\rm d}R/{\rm d}t$
) to become negative. Consequently, the middle droplet shrinks and feeds the remaining large droplet. The quantitative validation in the linear region is presented in figure 5(b), where the radius change rates are plotted against the thermodynamic driving force, calculated using the arithmetic mean radius. The data collapse into a well-defined linear curve. This result serves as a robust verification of our theoretical model, confirming that the radius change rate follows the predicted linear scaling law and that the critical radius is correctly defined by the arithmetic mean, independent of the droplet number or spatial configuration.
(a) Time evolution of the droplet radius change rates. (b) Dependence of the radius change rate on the thermodynamic driving force. The red dashed line represents the best linear fit to the combined simulation data of the three droplets.
Figure 5. Long description
The image contains two line graphs. The first graph on the left shows the time evolution of the droplet radius change rates for large, middle, and small droplets. The x-axis represents time, while the y-axis represents the rate of change of the droplet radius. The red solid line represents the large droplet, the blue dashed line represents the middle droplet, and the yellow dotted line represents the small droplet. The second graph on the right shows the dependence of the radius change rate on the thermodynamic driving force. The x-axis represents the inverse of the critical radius minus the inverse of the droplet radius, while the y-axis represents the rate of change of the droplet radius. The red circles represent simulation results, and the red dashed line represents the best linear fit to the combined simulation data of the three droplets. The graph is divided into regions for large, middle, and small droplets.
4. Conclusions
In this study, we employed the Enskog–Vlasov kinetic model to investigate droplet ripening in a two-dimensional single-component liquid–vapour system. This kinetic approach effectively captures non-equilibrium interfacial phenomena, which are inherently overlooked by macroscopic models, while offering greater computational efficiency than molecular dynamics. Through theoretical analysis and numerical simulations, we demonstrated that mass transport is fundamentally governed by the intrinsic kinetics of molecular exchange at the liquid–vapour interface when interfacial effects play a dominant role. Specifically, the net mass flux arises from the imbalance between evaporation and condensation rates. Consequently, the droplet-radius change rate exhibits a linear dependence on the thermodynamic driving force. We observed that this linear regime persists until the droplet radius becomes comparable to the interface width, where the overlap of diffuse interfaces leads to a reduction in liquid density and subsequent thermodynamic destabilisation. Furthermore, the kinetic coefficient of the linear relationship is found to follow an empirical power-law scaling with temperature. Importantly, the validity of the proposed theoretical model, including the determination of the critical radius, is confirmed to be robust across systems with different numbers of droplets and spatial configurations. Looking ahead, future research efforts will focus on extending the current kinetic framework to three-dimensional systems to accurately reproduce real-world experimental scenarios.
Declaration of interests
The authors report no conflict of interest.
T~=0.6
t~=50
y=H/2
1/R~c−1/R~
t~=4300
(1/R~cdiff−1/R~)/R~
T~=0.55, 0.6
0.65
K
T~
K=0.525T~6.52
