A weighted particle scheme for Enskog-Vlasov equation to simulate spherical nano-droplets/bubbles

The Enskog-Vlasov equation has been proven successful in capturing the complex behaviour of ﬂuids undergoing a phase change. However, its numerical solution is computationally demanding, and this has restricted the studies to one-and two-dimensional planar ﬂows. In this work, a weighted particle scheme is developed for the numerical solution of the Enskog–Vlasov equation in spherically symmetric geometry. It is shown how weighted schemes designed for the Boltzmann equation can be extended to cope with the non-local structure of the Enskog collision integral, and a compact expression of the mean force ﬁeld is determined using the shell theorem. As an application, the growth rates of nano-droplets/bubbles in the bulk of a homogeneous metastable vapour/liquid are evaluated in a wide range of supersaturation ratios that was not possible until now due to the high computational cost required by alternate approaches. The proposed scheme signiﬁcantly broadens the range of problems that can be investigated via the Enskog-Vlasov equation, and it is a stepping stone towards the simulation of general three-dimensional liquid-vapour ﬂows.


Introduction
The formation of liquid droplets from vapour condensation and the appearance of vapour bubbles within liquids in cavitation/boiling are ubiquitous in natural phenomena [1][2][3] and have numerous applications, ranging from cloud chambers for visualizing the passage of ionizing radiation [4] to accelerated healing of bone fractures [5] and targeted drug delivery [6].
The Classical Nucleation Theory (CNT) provides the simplest theoretical framework to rationalize these nucleation processes [1].The basic idea is that the onset of a nucleus inside a homogeneous bulk phase requires the formation of an interface with an associated energy penalty which the system must overcome for the nucleus to grow into a critical size and trigger the phase transition.In CNT the interface is supposed to be sharp, but more recent approaches based on density functional theory [7] remove this assumption by treating the interface as a continuously changing profile from the value at the centre of the nucleus to the bulk density far from it [8].However, notwithstanding the considerable theoretical ef-forts, the width of the energy barrier and the waiting time of nucleation for specific systems can only be predicted at the qualitative level, especially for multi-component substances [1].
Experimental apparatuses and techniques for studying nucleation processes have been considerably improved over the years, but accurate measurements are still challenging.Indeed, nucleation is a rare event that is associated with a very long waiting time and, once it starts to happen, the process itself unfolds quickly.Furthermore, experiments cannot provide insights into the detailed nucleation mechanism on the molecular level.
These theoretical and experimental challenges have spurred an increased interest in "numerical experiments" performed by means of molecular dynamics (MD) simulations.As the time scale required for the critical nuclei to form exceeds the MD time scale which is limited to nanoseconds, a seeding technique is typically employed, namely nuclei of different sizes are inserted beforehand in the metastable fluid at a fixed temperature, and the critical nucleus is identified as the one that neither grows nor collapses.Using this seeding technique, the homogeneous nucleation and subsequent growth rates have been investigated for droplets in metastable vapour [9][10][11][12].Bubbles growth in superheated liquid has also been studied both with reference to homogeneous nucleation [13][14][15][16] and plasmonic resonance which induces cavitation dynamics around intensely heated solid [17][18][19].Although MD simulations have been shown to provide invaluable insights into the homogeneous nucleation processes, they are also plagued by serious limitations, mostly the high computational burden that prevents systematic investigations.
In recent years, the Enskog-Vlasov (EV) equation has been proven a viable alternative to MD to study fluid flows undergoing a phase transition [20][21][22].This kinetic equation provides a statistical description of fluids via the distribution function in space and time of atoms interacting through the Sutherland potential resulting from the superposition of the repulsive hard sphere core and an attractive smooth tail [23].The short-range repulsive forces are treated as in the standard or revised Enskog equation theory [24], whereas the long-range attractive forces enter the equation in the form of a mean-field Vlasov term.The EV equation has been extensively used in mathematical studies of the liquid-vapour phase transition [25][26][27][28] and to elucidate the physics behind evaporation processes of simple fluids composed of monatomic [22,[29][30][31] and polyatomic [32,33] molecules, as well as mixtures [34,35].
The numerical solution of the EV equation is typically achieved by a particle scheme that extends the classical DSMC to account for the non-local structure of the Enskog collision integral [36,37].The computational effort required by this scheme is two orders of magnitude less than that of MD for simulating one-and two-dimensional flows [38], and substantial computational savings are envisioned for three-dimensional flows as well, provided suitable approximations of the mean force field are used [39].
As the growth of a small liquid/vapour nucleus in the bulk of a homogeneous metastable vapour/liquid can be treated as a one-dimensional spherically symmetric flow, the EV equation appears to be ideally suited to tackle this process.However, the particle numerical scheme that is currently employed for the solution of the EV equation must be properly redesigned to address a peculiar challenge posed by the study of flows with spherical symmetry.Indeed, unlike problems in planar geometry, there is no freedom to choose the cross-section area of the simulation box, and, therefore, the number of the simulation particles must match the number of the real atoms.This leads to a large computational burden as, in spherical geometry, the number of simulation particles increases with the cube of the radial coordinate.In order to circumvent this issue, weighted particles can be adopted where one simulation particle represents many real atoms.Weighted particle schemes have been developed for the Boltzmann equation [40], but their extension to the Enskog collision integral is not straightforward due to its non-local structure.
The present paper aims at developing a weighted particle scheme for solving spherically symmetric liquid-vapour flows described by the EV equation.This scheme is then used to investigate the growth of a small liquid/vapour nucleus in the bulk of a homogeneous metastable vapour/liquid.
The key novelty of the paper is that, for the first time, a numerical scheme is developed to solve the EV equation (a) in spherical symmetry and, (b) by employing weighted particles.The systematic study of the droplet/bubble growth rates is an additional novelty worth being stressed.Note that the proposed numerical scheme has a broad interest not limited to studies on nucleation, as the Enskog theory is increasingly used for many applications, including shale gas exploitation [41,42] and fluidized beds [43] to cite a few, and paves the way to the simulation of three-dimensional dense and/or liquid-vapour flows.
The rest of the paper is organized as follows.Section 2 presents the mathematical formulation of the EV equation for spherically symmetric fluid flows.Section 3 discusses the weighted particle numerical scheme.Section 4 outlines the simulation setup.Section 5 provides the validation of the code.Section 6 presents the results of droplet growth in supercooled and superheated vapour, and bubble growth in superheated liquid.Section 7 summarizes and comments on the main findings of the paper.

Mathematical formulation
Let us consider a fluid composed of identical spherical atoms of mass m and diameter a interacting through the Sutherland potential [23] defined by the superposition of the hard sphere potential and an algebraic attractive soft potential tail: where ρ = ||r 1 − r 2 || is the distance between the interacting atoms at position r 1 and r 2 , and the two positive constants φ a and γ define the depth of the potential well and the range of the soft interaction, respectively.Note that an exponential attractive tail can also be considered [30,37].The fluid can be described statistically via the distribution function which gives the number of atoms per unit volume in the single-particle phase space.A closed-form evolution equation for the distribution function can be derived by making two simplifying assumptions, namely long-range particle correlations are neglected and short-range particle correlations are approximated by following the Enskog theory of dense gases [20][21][22]25].By using spherical coordinates and assuming spherical symmetry, the one-particle distribution function simplifies to f (r, v r , v t , t), being r the radial space component, v r and v t represent the radial and tangential velocity components, respectively.The evolution equation of the particle distribution thus reads: where F r [n] is the radial component of the self-consistent force field generated by the soft attractive potential tail, which depends on the number density field n(r, t), and C E [ f ] is the hard-sphere collision integral (square brackets are used to highlight that these quantities are functionals of their arguments).
The radial component of the self-consistent force field can be computed by using the shell theorem.For ease of reference, we report here the expression obtained for γ = 6 as this value will be used in the following simulations: In Eq. (2b), the first term on the right-hand side refers to the inner and outer shells' contributions, whereas the second term refers to the contribution of spherical shells cut out by a sphere of radius a around r (see Appendix A for more details).The hard-sphere collision integral is given by: where (•) + indicates that the surface integral is restricted to the half-sphere for which v rel • k > 0, r is the unit vector in the radial direction, and χ [n] is the contact value of the pair correlation function in a hard-sphere fluid in equilibrium with a given number density [24] (see also Appendix B).At variance with standard Enskog theory (SET), the pair correlation function at contact is computed by following the Fisher and Methfessel prescription [44] as a functional of n: where In Eq. (3a), χ SET is evaluated by using the approximated equation of state of the hard sphere fluid proposed by Carnahan and Starling [45]: Note that, although SET suffers from some theoretical deficiencies (e.g., the proof of the H-theorem is still lacking), numerical solutions show a very good agreement with exact molecular dynamics simulations (see Refs. [26,46] for a more detailed discussion).

Numerical scheme
In this work, the spherically symmetric EV equation, Eqs.(2), is numerically solved by using a particle scheme.
The simulation domain is a one-dimensional interval [0, R w ] taken along a fixed radial direction that, without loss of generality, can be identified with the z-axis, and boundary conditions are imposed at r = R w (please refer to Sec. 4 for further details).The domain is divided into N c cells of equal length which represent spherical shells of thickness R w /N c .
The distribution function is approximated by N p simulation particles represented by Dirac's δ functions: where r i , v i , and w i are the position, the velocity, and the weight of the ith particle at time t.Particles' weights give the number of real atoms represented by each simulation particle and are assumed to be equal to the volume of the spherical shell for the ease of the algorithm design.
As stressed in the introduction, using weighted particles is crucial for simulating spherically symmetric fluid flows.Indeed, unlike the planar one-dimensional case, in spherical geometry one cannot freely choose the number of simulation particles, and match the number of the real atoms by a proper choice of the cross-section area of the computational domain.Instead, the number of simulation particles must be equal to the number of atoms which scales as O(R 3 w ).Since the computational cost increases linearly with the number of simulation particles, it is apparent that one simulation particle must represent many real atoms to keep the computational burden at a reasonable level when simulating flows in large domains.
To fully leverage the one-dimensionality of the spherically symmetric flows, the simulation particles are initially distributed along the same radial direction, and the algorithm is designed to keep them constantly there.The distribution function given by Eq. ( 4) is updated by a fractional-step method based on the time-splitting of the evolution operator in two sub-steps, namely free streaming and collision.Note that the design of weighted particle schemes entails some additional complexities over standard unweighted schemes.Indeed, in the streaming step, particles crossing the boundary of a cell must be removed or duplicated depending on whether they move into a cell with a larger or smaller weight, respectively.Furthermore, as discussed below, in the collision step, if particles with different weights collide, the velocity of the particle with a larger weight must be changed with a probability defined as the ratio of the particles' weights.
Numerical schemes that fully address these difficulties were designed for the Boltzmann equation [47,48], but none copes with particles' interactions modelled via the nonlocal Enskog collision operator.The extension of the streaming and collision sub-steps to this case is described in the following subsections.

Free-streaming sub-step
In the free-streaming sub-step, the collisions between particles are neglected and the distribution function is advanced from t to t + t based on the equation: Equation ( 5) is solved by the velocity Verlet integration [32,49] which is a second order time marching scheme comprising two stages.In the first stage, the half-step velocity is evaluated using the value of the mean force field at time t: Note that the half-step tangential velocities remain unchanged as the mean force field only acts along the radial direction.
The radial position of each atom, r i (t + t), is then updated based on the radial r i and tangential r ⊥,i displacements computed using the velocities given by Eqs.(6): where The radial and tangential velocity components given by Eqs. ( 6) must then be transformed because the simulation particles move off the initial radial axis whereas it is required that stay there for the numerical scheme to be truly one-dimensional.The velocity is thus decomposed along the new local radial and tangential components, and these latter are rotated to bring the radial velocity component back to the initial radial axis: Algorithm 1: Streaming sub-step.
Variables: N p : simulation particles w i : weight of particle i w j : weight of cell j for i = 1 to N p do Evaluate the mean force field, Eq. (2b); Compute the half-step velocity components, Eqs. ( 6); Compute the new radial position, Eq. ( 7); Bring the velocity components back to the radial axis, Eqs. ( 9); Compute the full-step velocity components, Eqs.(10); Identify the cell j occupied by particle i; if (w i > w j ) then Change the weight of particle i to w j ; Create Int[w i / w j ] − 1 copies of weight w j ; Create a copy of weight w j with probability: w i / w j − Int[w i / w j ]; else Delete particle i; Create 1 particle of weight w j with probability w i / w j ; end end vr,i = vr,i cos θ + vt,i sin θ, where θ = arccos ((r i + r i )/r i ).
Finally, the half-step velocity is updated to the full step velocity: where, as before, the tangential velocity components stay the same since the tangential component of the mean force field is zero.
As mentioned, weights are essentially attributed to cells, and particles inherit the weight of the cell where they are.
Accordingly, when a particle with weight w i moves to a cell with weight w j , two cases can be distinguished: if w i > w j , then the particle is decomposed in Int[w i / w j ] particles of weight w j , and another particle of weight w j is generated with probability w i / w j − Int[w i / w j ] (here Int[x] denotes the largest integer that does not exceed the magnitude of x); otherwise, if w i < w j , the particle is changed to a particle of weight w j with probability w i / w j [47].
The streaming sub-step is summarised by Algorithm 1.

Collision sub-step
In the collision sub-step, the positions of particles are frozen, and the distribution function is modified by evaluating the effect of short range hard-sphere interactions over the time interval t based on the equation: supplemented by the initial condition given by the solution of the streaming sub-step.The peculiar aspect of the collision process between weighted particles is that the velocity of the particle with a larger weight must be changed only in probability.Indeed, if the velocity components of the collision partners are updated regardless of their weights, every collision between simulation particles would correspond to a number of collisions between atoms given by the product of the particles' weights instead of the minimum between the two as it should be.The more natural probabilistic rule for updating the velocity of the colliding pair is to assume that the particle with lower weight always changes its velocity, whereas the particle with larger weight does it with a probability equal to the ratio of the two weights [47]: It is worth stressing that, under equilibrium conditions, Eqs. ( 12) conserve linear momentum and energy in a statistical sense, i.e. the conservation of these properties is only approximately satisfied over a sufficiently large number of collisions.In DSMC simulations for solving the Boltzmann equation using weighted particles, a split-merge algorithm has been developed to explicitly enforce momentum conservation at each collision, albeit with the associated cost of an energy loss proportional to w j /w i (w i − w j ).The energy can be conserved on average by adding the energy loss in one collision to the centre of mass of the subsequent collision pair, on a cell by cell basis [48].However, this solution cannot be straightforwardly implemented in the present case due to the non-local collision structure of the collision operator.Indeed, collisions involve particles located in two different cells and, therefore, one should define some criterion to redistribute the energy loss between these two cells.
In order to account for the probabilistic nature of collisions expressed by Eqs.(12), the collision frequency of simulation particles must be computed based on a collision cross section that is different from the real one.The key point is that two colliding weighted particles can be regarded as two molecular beams composed of w i and w j atoms.Therefore, the expected number of collisions between the atoms of the beams per unit of time is given by: N c = w i w j σ ij v rel , (13) where σ ij is the total cross section of atoms, and v rel is the relative speed of the two beams.As weighted particles are made to collide only in probability, the expected average number of collisions based on Eqs. ( 12) is: where σ ij is the total cross section of the weighted particles.By equating Eqs. ( 13) and ( 14), it follows that: Except for the extra-factor that appears in the total cross-section of the weighted particles, the collision sub-step closely follows the one developed in Ref. [36].The key quantity to evaluate is the expected number of collisions in the simulation box over the time step t, namely: where M i is the number of cells containing a portion of the sphere of radius a centred at the position of the particle i, N p,m is the number of particles in cell m, and where In Eq. ( 16), N c,i denotes the expected number of the collisions that particle i suffers in the time step t regardless of the position of the collision partner, whereas N c,i,(m) and N c,i,(m, j) represent the expected number of the collisions of particle i during the same interval of time with particles in the cell m, and with the particle j located in the cell m, respectively.Note that, as expected, if weights are constant, i.e. w i = w j , ∀i, j, then Eq. ( 16) simplifies to Eq. ( 9) of Ref. [36].The direct evaluation of the total number of collisions N c is computationally demanding, as it is proportional to N 2 p .Therefore, N c is estimated by a stochastic process based on a majorant frequency method.The following upper bound can be easily obtained: where A i and C i satisfy the following inequalities: The resulting Nc is obviously greater than the real number of collisions, but it is less computationally demanding to evaluate.The difference will be represented by false collisions that do not alter the velocities of the colliding particles.Since the term Nc,i,(m, j) does not depend on index j, an upper bound of the total number of collisions (real and false) in the i-cell is given by: Algorithm 2: Collision sub-step.
Variables: Nc : total (real+false) number of collisions Evaluate Nc , Eq. ( 22) ; for i = 1 to Nc do Select a particle in probability, Eq. ( 23); Generate a vector k uniformly over the unit sphere; Select the collision partner at random in the cell pointed by k; Hence, the total number of computational collisions (real + false) reads: The per-cell values of the constants A i and C i are shared by all particles in the same cell, and are updated at each time step to avoid an excessive number of false collisions.According to Eq. ( 22), the probability for a particle to be selected for a collision (real or false) is given by: Nc . ( Once a particle is selected, a vector k is drawn uniformly from the unit sphere, and the collision partner is chosen at random from the cell corresponding to the radial shell pointed by it.In general, this unit vector identifies a point off the radial axis (z-axis) and, therefore, the velocity components of the selected particle must be rotated to bring its z-component to be on the radial axis going through that point.The rotations are performed using the rotation matrix associated with the Euler angles.This process is depicted in Fig. 1, and more details about the rotation of the velocity components are given in Appendix B. The collision between particle i and particle j is real with a probability given by: where Since the inequality 0 ≤ ( k) ≤ 1 holds for ( k), Eq. ( 24) can be interpreted as the probability that a random number uniformly distributed over the unit interval is less than ( k).If the collision is accepted as real, the velocities of the collision pair are finally updated according to Eqs. (12).The collision sub-step is summarised by Algorithm 2.

Computational setup
In all the simulations discussed below, the interaction parameters were chosen to be φ a /(k where k B is the Boltzmann constant, and γ = 6 to match the same far field behaviour of the 12-6 Lennard-Jones potential [50].
The domain [0, R w ] was divided into a mesh of spherical shells of constant size r = a/10.This equidistant discretisation was adopted to maintain the same resolution throughout the system, which is useful to map out the growth curves of droplets/bubbles.In order to mimic equilibrium conditions in the far field, the particles that are within one molecular diameter from the outer boundary of the simulation domain, r i > R w − a, are not made to collide but rather their velocities are directly re-sampled from a Maxwellian distribution with the initial temperature of the system using the Box-Muller algorithm.
Although in principle the cell size and the cell weight can be selected independently, the weights were set proportional to the shell volume, i.e., w j = αV j .This ensures that when the density is constant, each shell contains the same number of simulation particles, resulting in a constant level of statistical noise throughout the domain.Unless otherwise stated, the value of α was set to 10 −4 to provide accurate results even with short time-averaging intervals.
The time-step was assumed to be t = 10 −2 a/(R g T 0 ) 1/2 , where T 0 is the reference temperature and R g = k B /m is the specific gas constant.The time-averaging intervals are chosen to be of length a/(R g T 0 ) 1/2 (one time unit) which permits one to average over 10 2 samples.The reference time is t 0 = a/(R g T 0 ) 1/2 , while the reference pressure is given by p 0 = n 0 a 3 R g T 0 .
The size of the simulation box value was chosen differently based on the problem such that results were not marred by finite size effects, e.g. from 35a (growth/collapse of liquid droplets in metastable vapour) to 125a (growth of vapour bubbles in metastable liquid).
The typical runtime for a simulation is about 2-3 h (N p = 4 − 6 × 10 5 simulation particles) for a droplet growth in metastable vapour and about 11-15 h (N p = 3 − 3.6 × 10 6 simulation particles) for bubble growth in metastable liquid on a single core of an Intel(R) Xeon(R) Gold 6230 CPU running at 2.10 GHz.

Code validation
In this section, three test cases are presented to provide the code validation.Firstly, a sample simulation of a droplet growing in its metastable vapour is considered to assess how much the conservation properties of the numerical scheme and the computational time are affected by the weights of the simulation particles (an in-depth analysis of the droplet growth fluid dynamics is presented in Section 6.1).Afterwards, it is shown that the code provides a very good prediction of the fluid surface tension using the Young-Laplace equation as well as of the critical size of droplet nuclei in their vapour as a function of the supersaturation ratio.

Case study I: Properties of the weighted particle scheme
A spherical droplet of radius R in = 5a was considered, centred at the origin of a system with size R w = 35a and surrounded by its vapour.The liquid density was set equal to the value predicted by the equilibrium coexistence (binodal) curve at T in /T c = 0.729, while the vapour density to a value n v /n 0 = 0.073 which is larger than the equilibrium one at that temperature (i.e., it corresponds to a supersaturation ratio of S = 2.21 -see Section 6.1).Fig. 2(a) shows the simulation results using simulation particles with constant and space-varying weights.The number of simulation particles was varied by two orders of magnitude to assess the sensitivity of the results on this parameter.A very good agreement is obtained, and this showcases that the weighted scheme can capture the system evolution even using a

Table 1
Validation case study I: Number of simulation particles, N p , error on the density field, n , and computational time, T r , of the runs reported in Fig. 2; n was computed as the 2 -norm error assuming the equally-weighted solution as baseline for comparison.number of simulation particles one hundred times smaller (i.e., increasing the particles' weights by a factor of 10 2 ).Note that using weighted particles leads to considerable savings in computational time while keeping a high level of accuracy as illustrated in Table 1.
It is worth stressing that the proposed numerical scheme conserves mass and kinetic energy only in a statistical sense.Indeed, during the streaming sub-step, the weights of the simulation particles crossing cells' boundaries are changed in probability and, therefore, the mass of the system varies.The kinetic energy of the system fluctuates as well because it is not explicitly conserved by both individual collisions and the subsequent rotation of the post-collisional velocities needed to preserve the spherical symmetry of the system.By contrast, the momentum of the system can be considered to be constantly equal to zero despite the fact that the momentum of the particles is not conserved in the streaming nor in the collision sub-steps because, due to the imposed spherical symmetry, one may think that each simulation particle is evenly distributed in the three-dimensional spherical shell to which it belongs.The time variation of mass and kinetic energy is shown in Fig. 2 for three values of the weighting factor α = {10 −2 , 10 −3 , 10 −4 } (introduced in Section 4).As expected, the fluctuations of these quantities with respect to the initial values are directly proportional to α (i.e., inversely proportional to the number of simulation particles), and they turn out to be within 1% when α = 10 −4 which is the value chosen to run the simulation campaign presented in Section 6.

Case study II: Fluid surface tension
According to the Young-Laplace equation, the pressure difference across the curved interface formed between a spherical liquid droplet and its vapour reads: where p and p v are the bulk pressures in the centre of the droplet and far outside, respectively, R is the radius of the spherical droplet, and γ s is the fluid surface tension.
In this subsection, γ s is estimated based on Eq. ( 26) by computing p for liquid droplets of different R in equilibrium with their vapour.Simulations were carried out for droplets centred at the origin surrounded by their vapour.Initially, the droplets have a sharp interface located at R/a = {6, 7, 8, 9, 10, 12, 14} and the liquid and vapour densities were set equal to the equilibrium values predicted by the coexistence curve at the selected working temperatures, i.e.T /T c = {0.729,0.795, 0.861}.numerical results based on the EV equation; Solid lines: predictions of Eq. ( 27).Solid lines are drawn up to the maximum supersaturation ratio, S max = p(T , n s )/p eq (T ), where n s is the density value on the spinodal line corresponding to each temperature.
The simulation domain was taken sufficiently large to mitigate possible boundary effects, i.e., r ∈ [0, R + 30a].The system evolution is computed until the liquid-vapour interface forms and the system equilibrates.The Andersen thermostat [51] is applied throughout the system during the simulation to keep the temperature of the system constant during the transient evaporation phase.This thermostat was chosen for its simplicity and computational efficiency.Fig. 3 shows the pressure difference as a function of the droplet radius for different temperatures (solid symbols), and the best linear fit of the results (solid lines) that, according to Eq. ( 26), provides the estimated surface tension.The values of the surface tension are listed in Table 2 and compared against the predictions reported in Ref. [30] that were obtained directly by computing the surface tension of a planar interface based on its mechanical definition and indirectly via the Young-Laplace equations for cylindrical droplets.An excellent agreement is found for all the considered temperatures.

Case study III: Droplet critical radius
The critical radius of a droplet is defined as the minimum radius that the droplet must have to be in a stable equilibrium with its vapour.According to CNT, the critical radius, R * , is given by: where n is the equilibrium liquid number density at the temperature T , and S = p(T , n v )/p eq (T ) is the ratio of the vapour pressure to the equilibrium vapour pressure, which is referred to as supersaturation ratio.The CNT prediction of the critical radius is in good agreement with the numerical results of MD [52] and Monte-Carlo [53] simulations of Lennard-Jones fluids for temperatures sufficiently lower than the critical temperature.
In this section, the critical radius is estimated numerically and compared against the theoretical prediction provided by Eq. (27).
Simulations were carried out for droplets centred at the origin surrounded by their vapour, and R * is computed by reducing the droplet radius up to the value for which the droplet starts to evaporate.The liquid density was set equal to the equilibrium values predicted by the coexistence curve at the selected working temperatures, i.e., T /T c ∈ {0.663, 0.729, 0.795}.The vapour density is based on supersaturated ratio S > 1 which is varied by increasing the density up to values close to spinodal curve at constant temperature (see also Fig. 5).The simulation domain was defined as sufficiently large to mitigate possible boundary effects.
As shown in Fig. 4, an excellent agreement is found for the smaller values of the temperatures, whereas the CNT prediction underestimates the numerical results at the highest temperature, which is not unexpected.Indeed, the width of the liquid-vapour interface increases as the temperature becomes closer to the critical temperature, and, accordingly, the sharp interface assumption of the CNT becomes increasingly inaccurate.

Droplets and bubbles growth rates
In this section, we present results of droplet growth (collapse) in supercooled (superheated) vapour (Subsection 6.1), and the growth of bubbles in superheated liquid (Subsection 6.2) for a representative set of temperatures of the system, i.e., T in /T c = {0.663,0.729, 0.795, 0.861}.
The initial conditions explored in the simulation campaign are represented by solid symbols in Fig. 5 alongside the binodal and spinodal curves.The binodal curve was obtained by numerical solving the two coupled nonlinear algebraic equations that enforce the conditions of equal pressure and chemical potential at constant temperature across the liquid and vapour phases [26].The spinodal curve was obtained by requiring a vanishing pressure gradient with respect to density at constant temperature [22,26].

Liquid droplets in metastable vapour
In this subsection, the growth (collapse) of liquid droplets in supercooled (superheated) vapour is investigated for different temperatures of the system.Initially, the liquid droplet was centred at the origin with a sharp interface at R in = 5a, and the value of the density was set to the one predicted by the equilibrium coexistence (binodal) curve at the temperature of the system.The droplet was surrounded by a vapour whose density is varied based on the supersaturation ratio (see Fig. 5).The length of the simulation domain was set to R w = 35a.
Figs. 6(a)-(c) show the time evolution of the density, radial velocity, and temperature profiles, respectively, for the exemplary scenario of a growing droplet in supercooled vapour at T in /T c = 0.663 and S = 3.67.During an initial transient stage, t/t 0 50, the liquid-vapour interface forms, and the droplet density decreases due to the heating up driven by the strong condensation of the vapour, as shown in Fig. 6(b) where the large negative radial velocity of the vapour is apparent.Afterwards, the temperature profile levels off, and the droplet enters the linear growth regime as reported in Fig. 6(c) and in Fig. 6(d) which shows the temporal evolution of the bulk droplet temperature as a function of the radial distance for different supercooling ratios S. Note that a similar temperature behaviour has been found based on MD simulations [54].
The droplet growth rate Ṙ was estimated from the best linear fit of the growth curve, namely the displacement of the liquid-vapour interface as a function of time, being the position of the liquid-vapour interface defined by the point where the density takes the mean of the liquid and vapour densities.Note that the step-like behaviour of the growth curve shown in Fig. 7(b) is due to the space resolution dictated by the cell length.The results for Ṙ are plotted in Fig. 7(b) with respect to the supersaturation ratio S for different temperatures of the system.It is worth stressing that the growth curves do not go through the origin (i.e., Ṙ = 0 for S = 1), as one could have expected.Indeed, the supersaturation ratio is defined based on the equilibrium vapour density provided by the coexistence curve and does not account for the effect of the Laplace pressure on the actual density of the liquid droplet.As it can be observed, the growth rate increases with the supersaturation ratio and decreases with the temperature.Similar qualitative conclusions have also been drawn in previous works based on simplified kinetic models [55], and MD simulations [12,56].However, it is worth stressing that a much smaller portion of the problem parameter space has been explored in MD simulations due to the high computational burden.For completeness, in Fig. 7 results are also included for collapsing droplets, i.e., droplets with a negative growth rate.

Vapour bubbles in metastable liquid
In this subsection, the growth of vapour bubbles in superheated liquid is investigated for different temperatures of the system.
Initially, the vapour bubble was centred at the origin with a sharp interface at R in = 5a, and the density was set to the one predicted by the equilibrium coexistence (binodal) curve at the temperature of the system.The bubble was surrounded by a liquid whose density is slightly smaller than the equilibrium one (see Fig. 5).Compared to simulations of droplet growth, a larger simulation domain was considered, i.e., R w = 125a, to mitigate the spurious effects due to liquid compressibility.
Figs. 8(a)-(c) show the time evolution of the density, radial velocity, and temperature profiles, respectively, for the exemplary scenario of a growing bubble in superheated liquid at T in /T c = 0.729 and liquid pressure p /p 0 = −0.1098.Note that this latter is not a gauge pressure (i.e. it is negative relative to the initial pressure), but rather an absolute negative pressure due to the fact that the system is in a metastable superheated state.From the molecular standpoint, in these states the force field due to the particles' long-range attractive potential tail overcomes the net repulsive force due to hard-sphere collisions and leads to negative pressure values.
During the initial rapid expansion stage, the bubble bulk density takes a value smaller than the equilibrium one, as also observed in MD simulations [16].The temperature drops to about 20% of the equilibrium value within the bubble (see Fig. 8(d)) whereas overshoots this value at the interface by as much as 12 − 13% in the case of very fast growing bubbles.
This behaviour can likely be attributed to the local rapid liquid compression [16].As shown in Fig. 8(b), the radial velocity profile exhibits a maximum near the interface, that initially increases and then reaches a quasi-constant value when the bubble enters the linear growth regime.Note that the numerical results almost perfectly match the theoretical predictions provided by the model proposed in Refs.[57,58]: where x = R(t)/R w , being R(t) the radius of the bubble at time t.
Fig. 9(a) shows the temporal evolution of the bubble radius for different liquid pressures.It is apparent that, after the initial transient, the bubble enters the linear growth regime.As done in the droplets case, the growth rates are estimated by the slope of the linear best fit of the growth curves, and the results are reported in Fig. 9(b) with respect to the external liquid pressure p for different temperatures.The average growth rate decreases as the liquid pressure p increases and reduces with the temperature.As discussed in Appendix C, for the higher negative pressures, the bubble's growth rate can be estimated based on the Rayleigh-Plesset theory, and, as expected, the closer the liquid pressure is to the spinodal value, the better the agreement is between the theoretical prediction and the numerical result.Figs.10(a-c) show the temperature profiles near the interface for T in /T c = {0.663,0.729, 0.795} and various values of the surrounding liquid pressures.The maximum values of the temperature are attained close to the interface and decrease as the liquid pressure p becomes lower, whereas the vapour temperature remains relatively constant.Temperature profiles flatten as the simulation temperature T increases.An interesting behaviour is found when comparing relative temperature profiles obtained at different equilibrium temperatures but similar values of the growth rate.As it can be seen in Fig. 10 (d) all the results superimpose, suggesting that the relative temperature is well characterised by the bubble growth rate.

Conclusions
In this work, for the first time, the EV equation is solved in spherical geometry using a weighted particle scheme.Particle weights must be introduced to circumvent the large number of simulation particles that would be otherwise required as, in spherical geometry, the number of simulation particles scales with the volume of the simulation domain, which in turn is proportional to the cube of the radius.The collision algorithm for weighted particles is designed to cope with the nonlocality of the Enskog collision term, and a compact expression of the Vlasov mean force field is determined using the shell theorem.
For a sample simulation of the liquid droplet evolution in its metastable vapour, it was shown that the weighted particle scheme matches the simulation results obtained by using simulation particles with constant weights even using a number of simulation particles considerably lower.Furthermore, using weighted particles leads to a great reduction in computational time.
Two benchmark cases were further considered for code validation.First, the surface tension was estimated indirectly using the Young-Laplace equation and the numerical results were compared against values reported in the literature.Second, the critical radius of a droplet in its vapour was evaluated, and the numerical results were compared against the predictions of CNT.In both cases, an excellent agreement was found, thus providing a sound validation of the numerical scheme.
An extensive simulation campaign was then carried out to investigate the growth of droplets/bubbles in metastable vapour/liquid.
The droplet growth was found to be driven by vapour condensation.This causes an initial rapid increase of the temperature and a drop of the density in the bulk of the droplet.Afterwards, the liquid bulk temperature levels off, and this marks the beginning of the linear growth regime.The growth rate increases with the supersaturation ratio, and decreases with the temperature of the system.The results of droplet growth in metastable vapour were found in agreement with previous works that used simplified kinetic models [55] and with MD simulations [12,56].
The bubble growth in superheated liquid was found to be driven by the vapour expansion.This compresses the surrounding liquid and leads to an increase in the density and temperature close to the interface.As in the droplets case, after an initial transient, bubbles grow linearly in time.The growth rate decreases with the liquid pressure and with the temperature of the system.For states close to the spinodal curve, the numerical results match the predictions of the Rayleigh-Plesset equation in the limit of linear growth dominated by the liquid pressure term.Interestingly, the relative temperature profiles corresponding to different initial temperatures but similar growth rate, was found to superimpose each other, thus suggesting that the growth rate characterises the liquid behaviour during the bubble growth.
This work extends the range of liquid-vapour flows that can be studied based on the Enskog-Vlasov, and broadens the range of problems that can be tackled by the Enskog equation.Indeed, weighted particles boost the efficiency of simulations and permit one to tackle problems that would not be feasible otherwise, e.g., three-dimensional problems and/or problems with a large computational domain.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.shells.cutout sphere has its origin on a cell centre and a diameter of 2a.In this figure the value of the radial spacing is r = a, for better visual interpretation (the value adopted in the main text is r = a/10).

Appendix B. Collision dynamics
As discussed in Section 3.2, once a particle is selected for a collision, a random vector k is drawn uniformly on the unit sphere centred at the particle's position, and a collision partner is picked up with uniform probability from the cell corresponding to the spherical shell pointed by k.Afterwards, the velocities of the particles are updated according to the collision rule: where (v 1 , v 2 ) and (v 1 , v 2 ) are the pre-collisional and post-collisional velocities entering the kinetic equation, respectively, and v rel = v 1 − v 2 is the relative velocity.
It is worth stressing that, before applying Eqs.(B.1), the velocity of the collision partner must be rotated in such a way that the tail of the velocity vector coincides with the tip of the unit random vector.The rationale behind this procedure is that, by symmetry consideration, particles at the same radial distance from the origin have the same velocity components, namely the velocity vectors are the same up to a rotation that brings particles to superimpose each other.In the Cartesian reference frame, the velocity of the selected collision partner must thus be multiplied by the rotation matrix based on the Euler's angles corresponding to that rotation: where k x , k y , k z are the components of the unit random vector, and z 1 , x are the radial positions of the collision pair.
Finally, the velocity of the collision partner after the collision must be rotated back to bring radial component to the z-axis: Note that, due to symmetry one may use Eq.(B.4) to compute the radial component whereas the magnitude of the tangential one is given by: ṽ 2t = v 2  2 − ṽ 2 2r . (B.5)

Fig. 1 .
Fig.1.Selection of the collision pairs.The radial direction is identified with the z-axis.A particle at z i is selected for collision from the spherical shell (b).A random unit vector k is drawn to find the collision partner (the red unit sphere defines the set of possible directions).The collision partner is then selected at random in the spherical shell (d) pointed by k, with radial coordinate r.The velocity components of the collision partner ṽ2 , originating on the radial axis (z-axis), are rotated to the endpoint of k using the rotation matrix R resulting in v 2 , Eq. (B.2).Afterwards, the elastic collision is performed, Eq. (B.1), and the post-collisional velocities components are rotated back to the radial axis, Eq. (B.4).(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Validation case study I: Sample droplet growth simulation.(a) Comparison of runs with different number of simulation particles and constant/varying weights.(b) Time evolution of normalised total mass and total kinetic energy for different values of α; m 0 and E 0 represent the total initial mass and total kinetic energy, respectively.

Fig. 5 .
Fig. 5. Binodal and spinodal curves for the EV equation.Solid symbols are the initial conditions considered in the simulation campaign.Top-pointing triangles: superheated vapour; bottom-pointing triangles: supercooled vapour; empty rhombi: superheated liquid.The insets represent the initial density profile and the one in the linear growth regime in their respective phase diagram region.Empty symbols represent pairs of temperature and density used in: Fig. 2 (•); Fig. 6(a-c) ( ); Fig. 8(a-c) ( ).

Fig. 6 .
Fig. 6.(a) Density, (b) radial velocity, and (c) temperature profiles during the growth of the droplet at T in /T c = 0.663, and S = 3.67.(d) Time evolution of the droplet bulk temperature at T in /T c = 0.663 for different values of the supercooling ratio S.

Fig. 7 .
Fig. 7. (a) Typical droplet growth curves at T in /T c = 0.729 for different supercooling ratios S, along with the corresponding best linear fit lines.(b) Average droplet growth rate Ṙ with respect to the supercooling ratio S = p(T , n v )/p eq (T ) for different temperatures.

Fig. 8 .
Fig. 8. (a) Density, (b) radial velocity, and (c) temperature profiles during the growth of the bubble at T in /T c = 0.729 and p /p 0 = −0.110.(d) Time evolution of the bubble growth at T in /T c = 0.729 for different values of the liquid pressure p .Dashed lines in (b) are the theoretical predictions provided by the model proposed in Refs.[57,58]; Solid lines in (b)-(d) represent fitting curves of the numerical data obtained using cubic smoothing splines.

Fig. 9 .
Fig. 9. (a) Typical bubble growth curves at T in /T c = 0.663 for different liquid pressures p , along with the corresponding best linear fit lines.(b) Average bubble growth rate Ṙ with respect to the surrounding liquid pressure normalised to its absolute spinodal value p s ≡ |p(T in /T c , n s )|.

Fig. 10 .
Fig. 10.Temperature profiles during the bubble growth at (a) T in /T c = 0.663, (b) T in /T c = 0.729, and (c) T in /T c = 0.795 for different liquid pressures.(d) Temperature profiles during the bubble growth with Ṙt 0 /a ≈ 0.3 for different temperatures.In all cases, the radius of the bubble is R 0 = 25a.

Fig. A. 11 .
Fig. A.11. Evaluation of the radial component of the self-consistent force field (pictorial representation).The solid lines represent the computational cell edges, while the dashed lines represent the cell The spherical shells are partitioned into 3 groups: inner (A.2), outer (A.3), and cutout (A.4)

Table 2
Validation case study II: Comparison of the surface tension computed directly based on its mechanical definition for planar interfaces, and indirectly via the Young-Laplace equation for cylindrical and spherical liquid droplets.
Fig. 3. Validation case study II: Pressure difference at the liquid-vapour interface, p, versus the reciprocal of the droplet radius, a/R, for different temperatures.