Thermally activated vapor bubble nucleation: the Landau-Lifshitz/Van der Waals approach

Vapor bubbles are formed in liquids by two mechanisms: evaporation (temperature above the boiling threshold) and cavitation (pressure below the vapor pressure). The liquid resists in these metastable (overheating and tensile, respectively) states for a long time since bubble nucleation is an activated process that needs to surmount the free energy barrier separating the liquid and the vapor states. The bubble nucleation rate is difficult to assess and, typically, only for extremely small systems treated at atomistic level of detail. In this work a powerful approach, based on a continuum diffuse interface modeling of the two-phase fluid embedded with thermal fluctuations (Fluctuating Hydrodynamics) is exploited to study the nucleation process in homogeneous conditions, evaluating the bubble nucleation rates and following the long term dynamics of the metastable system, up to the bubble coalescence and expansion stages. In comparison with more classical approaches, this methodology allows on the one hand to deal with much larger systems observed for a much longer times than possible with even the most advanced atomistic models. On the other it extends contin- uum formulations to thermally activated processes, impossible to deal with in a purely determinist setting.

Thermal fluctuations play a dominant role in the dynamics of fluid systems below the micrometer scale. Their effects are significant in, e.g., the smallest microfluidic devices [1,2] or in biological systems such as lipid membranes [3], for Brownian engines and in artificial molecular motors [4]. They are crucial for thermally activated processes such as nucleation, the precursor of the phase change in metastable systems. Nucleation is directly connected to the phenomenon of bubble cavitation [5] and of freezing rain [6], to cite a few. There, thermal fluctuations allow to overcome the energy barriers for phase transitions [7][8][9]. Depending on the thermodynamic conditions, the nucleation time may be exceedingly long, the so-called "rare-event" issue. Classical nucleation theory (CNT) [10] provides the basic understanding of the phenomenon which is nowadays addressed through more sophisticated models like density functional theory (DFT) [11,12] or by means of molecular dynamics (MD) simulations [13]. These approaches need to be coupled to specialized techniques for rare events, like the string method [14], the forward flux sampling [15] and the transition path sampling [16], to reliably evaluate the nucleation barrier and determine the transition path [17]. For many real systems they are often computationally too expensive and therefore limited to very small domains.
Here we adopt a mesoscopic continuum approach, embedding stochastic fluctuations, for the numerical simulation of thermally activated bubble nucleation. Since the pioneering work of Landau andLifshitz (1958, 1959) [18] several works contributed to the growing field of "Fluctuating Hydrodynamics" (FH) [19]. More recently the theoretical effort has been followed by a flourishing of highly specialized numerical methods for the treatment of the stochastic contributions [20][21][22][23]. The present model is based on a diffuse interface [24] description of the two-phase vapor-liquid system [25] similar to the one recently exploited by Chaudhri et al. [26] to address the spinodal decomposition. The thermodynamic range of applicability of this approach is subjected to some restrictions: i) at the very first stage of nucleation the vapor nucleii, smaller than the critical size, need to be numerically resolved; analogously, ii) the thin liquid-vapor interface needs to be captured for the correct evaluation of the capillary stresses; iii) fluctuating hydrodynamics predicts that the fluctuation intensity grows with the inverse cell volume, ∆V , leading to intense fluctuations, contrary to the assumption of weak noise needed to derive the model ( δf 2 / f 1). Notwithstanding these restrictions, where it can be applied, this mesoscale approach offers a good level of accuracy (as will be shown when discussing the results) at a very cheap computational cost compared to other techniques. The typical size of the system we simulate on a small computational cluster (200 × 200 × 200 nm 3 , corresponding to a system of order 10 8 atomistic particles) is comparable with one of the largest MD simulations [27] on a tier-0 machine. Moreover the simulated time is here T max ∼ µs to be compared with the MD T max ∼ ns. The enormous difference between the two time extensions allows us to address the simultaneous nucleation of several vapor bubbles, their expansion, coalescence and, at variance with most of the available methods dealing with quasi-static conditions, the resulting excitation of the macroscopic velocity field.
The diffuse interface modeling adopted here has a strict relationship with more fundamental atomistic approaches, since it is based on a suitable approximation of the free energy functional derived in DFT [12]. It dates back to the famous Van der Waals square gradient approximation of the Helmholtz free energy functional F [ρ, θ] = V dV (f 0 (ρ, θ) + 1/2λ∇ρ · ∇ρ), where arXiv:1801.06817v1 [cond-mat.soft] 21 Jan 2018 f 0 is the classical bulk free energy, expressed as a function of density ρ and temperature θ. λ is the capillarity coefficient that controls the (equilibrium) surface tension γ and interface thickness, e.g.
with ω 0 (ρ, θ) the bulk grand potential per unit volume and the superscript sat denoting the (temperature dependent) saturation conditions, see Supplemental Material [28]. The mesoscopic fields obey mass, momentum and energy conservation, with the addition of stochastic contributions (Lifshitz-Landau-Navier-Stokes equations with capillarity): where u is the fluid velocity, p = ρ 2 ∂(f 0 /ρ)/∂ρ is the pressure, E is the total energy density, E = U + 1/2ρ|u| 2 + 1/2λ|∇ρ| 2 , with U the internal energy density. In the momentum and energy equations, Σ and q are the classical deterministic stress tensor and energy flux, respectively, and the terms with the prefix δ are the stochastic parts, required to satisfy a suitable fluctuationdissipation balance (FDB). The deterministic stress tensor Σ and the energy flux q follow by standard nonequilibrium thermodynamic considerations [29] as: Σ = (λ/2|∇ρ| 2 +ρ∇·(λ∇ρ))I −λ∇ρ⊗∇ρ+µ[(∇u+∇u T )− 2/3∇ · uI], q = λρ∇ρ∇ · u − k∇θ, with µ and k the dynamic viscosity and the thermal conductivity, respectively. Enforcing the FDB, the covariance of the stochastic fluxes follows as, δΣ( where Q Σ αβνη = 2k B θµ (δ αν δ βη + δ αη δ βν − 2/3δ αβ δ νη ) and Q q αβ = 2k B θ 2 kδ αβ , with k B the Boltzmann constant (see the Supplemental Material [28] for details). Thanks to the Curie-Prigogine principle [30], the crosscorrelation between different tensor rank fluxes vanishes, Bubble nucleation is investigated in a metastable liquid enclosed in a cubic box with periodic boundary conditions, with fixed volume, total mass and energy (NVE). The equation of state (EoS) we use, which can be chosen freely among available models, e.g. van der Waals or IAPWS [31] EoS for water, corresponds to a Lennard-Jones (LJ) fluid [32] to allow direct comparison with MD simulations. Quantities are made dimensionless according to ρ * = ρ/ρ r , θ * = θ/θ r , u * = u/U r , by introducing as reference quantities the parameters of the LJ potential, σ = 3.4 × 10 −10 m as length, = 1.65 × 10 −21 J as energy, m = 6.63 × 10 −26 kg as mass, U r = ( /m) 1/2 as velocity, T r = σ/U r as time, θ r = /k B as temperature, µ r = √ m /σ 2 as shear viscosity, c vr = mk B as specific heat at constant volume and k r = µ r c vr as thermal conductivity. The only dimensionless control group is the capillary number C = λρ r /σ 2 U 2 r , fixed here to C = 5.244 in order to reproduce the LJ surface tension [32]. The symbol * will be omitted in the following to simplify notation.
The system volume V = (600) 3 has been discretized on an equi-spaced grid with 50 cells per direction. Critical nucleus size, interface thickness and transition time are crucial to prepare a well resolved simulation. The preliminary set-up was based on the expected transition time derived from the free energy barriers and on the size R c of the critical nucleus at the different thermodynamic conditions, as estimated from the string method [14] applied to the present system. Technical details are available in the SM [28]. The comparison with the predictions of CNT are summarized in Tab. I. It turns out that, for the present relatively large system, there is no quantitative difference between the barriers estimated from CNT in the different ensembles, and we may use the grand canonical expression ∆Ω CN T = 4/3πγR 2 c . After a convergence analysis we found that a grid size ∆x = 12 is sufficient for a reliable simulation in these thermodynamic conditions, see SM [28]. Thanks to the extension of the simulated domain, ten runs for each condition provide a reasonably well converged statistics.
Among the different conditions we have investigated, we mainly focused on the initial temperature θ 0 = 1.25 at changing bulk density to explore the corresponding metastable range ρ L ∈ [ρ spin , ρ sat ] = [0.44, 0.51], where ρ sat and ρ spin are the saturation and spinodal densities, respectively. A few snapshots of the evolution for two different initial conditions are shown in the top panels of Fig.s 1 and 2. Starting from a homogeneous metastable liquid phase, the fluctuations lead the system to spontaneously nucleate vapor bubbles. The nucleii start out with a complex, far from spherical, shape, see, e.g., [13]. Roughly, when they happen to reach a size larger than critical they typically expand. Eventually, after a long and complex dynamics where bubbles expand and co-  alesce, stable equilibrium conditions are reached. This new configuration is characterized by several stable vapor bubbles in equilibrium with the surrounding liquid. The case at ρ L = 0.46, the closest one to the spinodal we considered here, is the more populated, Fig. 1 in comparison with Fig. 2. This system has a lower barrier, hence it nucleates faster. The initial (metastable) thermodynamic condition also influences the number and typical dimension of the bubbles in the final stage, bottom panels of Fig.s 1 and 2 providing the bubble number N b as a function of time. A tracking procedure has been put forward to follow the evolution of the distinct bubbles. By monitoring volume, mass, center of mass and its velocity, the tracking algorithm allows to detect coalescence events (see [28] for the procedure details). The actual number of collisionsÑ coll evaluated at each time step is characterized by a highly discontinuous fingerprint, hence we smoothed it out with a Gaussian kernel with standard deviation on the order of 50 time units to gain more robust indications. The smoothed number of collisions N coll , plotted with the blue line in the bottom panel of the Fig.s 1 and 2, shows a strong correlation with the number of bubbles throughout the initial nucleating stagewhen N b grows linearly with time (i.e. at a constant nucleation rate) -and the first part of the expansion phase -characterized by a rapidly decreasing number of bubbles mainly due to collapse -that we will call collapsing stage. These two stages are characterized by a competing-growth mechanism [33] due to the constraint of constant mass, explaining the high number of supercritical bubble collapses. The coalescence events start being less and less probable during the slowly-expanding regime that characterizes the long-time dynamics of the multi-bubble system. The inset of the Fig. 1 zooms into this regime showing that isolated collision events are still occurring, contributing to important acceleration toward the final equilibrium condition. The volume history of the distinct bubbles (in particular those survived up to the last time investigated) have been plotted in Fig. 3. Among the different bubble evolutions, we highlighted in red the volume histories of those bubbles that experienced intense coalescence events, characterized by a sudden increase in volume. It is apparent that the larger bubbles gained substantial part of their volume by coalescence. To substantiate this impression, for each bubble in the last configuration, the sum of the volumes acquired by coalescence throughout the whole evolution, V coll , was calculated, inset of Fig. 3.
The present mesoscale approach allows to access the statistics of bubble dimensions. The probability distribution function of bubble volumes f (V ) is plotted in Fig. 4. During both the nucleating and collapsing stages the pdf is sharply peaked at small volumes, of the order of 2-4 V c . The successive bubble expansion phase is substantially slower and calls for a much longer observation time to detect a significant growth (green dash-dotted curve at t = 34000). The intense coalescence events explain  the presence of the second peak in the pdf at very large volume (black curve in the inset of Fig. 4 at t = 163760).
The initial nucleating stage, where the bubble number increases linearly, gives access to the nucleation rate J in terms of bubbles formed per unit time and volume. It is here calculated as the slope of the linear fit to the curves of Fig.s 1 and 2 near the origin. The CNT expression for the (dimensional) nucleation rate is J CN T = n L 2γ/mπ exp(− ∆Ω/k B θ), where n L is the liquid number density. It is compared with the measured one in Fig. 5 which also provides the comparison with some MD simulations [13,34]. Our results are in reasonable agreement with molecular dynamic simulations in the NVE ensemble [13], as shown in the inset of Fig. 5, and consistent with the order of magnitude predicted by CNT. The apparent discrepancy that, contrary to expectation, the nucleation rate is smaller where the barriers are smaller than CNT, may be explained by the compression induced by bubble nucleation in the NVE ensemble.
In conclusion, the FH approach together with a diffuse interface modeling of the multiphase system have been exploited to study homogeneous nucleation of vapor bubbles in metastable liquids. We evaluated the nucleation rate and compared it favorably with state of the art simulations and theories. The present technique has revealed extremely cheaper with respect to MD simulations, allowing the analysis of the very long bubble expansion stage where bubble-bubble interaction/coalescence events turn out to determine the eventual bubble size distribution. The accurate results and the efficiency of the modeling encourage the exploitation to more complex conditions, like e.g. heterogeneous nucleation and multi-species systems, and could pave the way for the development of innovative continuum formulation to address thermally activated processes.
The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/ 2007-2013)/ERC Grant agreement no. [339446].