Modelling of free bubble growth with Interface Capturing Computational Fluid Dynamics

This paper presents simulations of the growth of stationary and rising vapour bubbles in an extend pool of liquid using an Interface Capturing Computational Fluid Dynamics (CFD) methodology coupled with a method for simulating interfacial mass transfer at the vapour-liquid interface. The model enables mechanistic prediction of the local rate of phase change at the vapour-liquid interface and is applicable to realistic cases involving two-phase mixtures with large density ratios. The simulation methodology is based on the Volume of Fluid (VOF) representation of the flow, whereby an interfacial region in which mass transfer occurs is implicitly identified by a phase indicator, in this case the volume fraction of liquid, which varies from the value pertaining to the “bulk” liquid to the value of the bulk vapour. The novel methodology proposed here has been implemented using the Finite Volume framework and solution methods typical of “industrial” CFD practice embedded in the OpenFOAM CFD toolbox. Simulations are validated via comparison against experimental observations of spherical bubble growth in zero gravity and of the growth of a rising bubble in normal gravity. The validation cases represent a severe test for Interface Capturing methodologies due to large density ratios, the presence of strong interfacial evaporation and upward bubble rise motion. Agreement of simulation results with measurements available in the literature demonstrates that the methodology detailed herein is applicable to modelling bubble growth driven by phase-change in real fluids.


Introduction
Advances in the simulation of two-phase flows with interface-resolving techniques, enabled by increasing computer power, have created the possibility to model two-phase processes such as boiling at the scale of the single bubbles directly via predicting from first-principles the basic interfacial phenomena that give rise to mass transfer. One crucial aspect is the requirement to develop mechanistic modelling of interfacial mass transfer due to phase change that is general enough to be coupled with popular Interface Capturing methodologies such as the Volume of Fluid (VOF) or Level Set methods (Tryggvason et al., 2011) and incorporated into broadly applicable Computational Fluid Dynamics (CFD) "toolboxes" such as the OpenFOAM code.
Whereas the vast majority of simulation methodologies for flows with mass transfer have been implemented in application-specific computational codes of limited applicability, development of such modelling capabilities in general calculation procedures (e.g., established general purpose CFD codes such as OpenFOAM) is still a real challenge. One possible modelling approach was demonstrated in Giustini and Issa (2021), where the authors use the standard Finite Volume discretization framework to demonstrate the applicability of interface capturing simulation with mass transfer to modelling boiling flows on arbitrary mesh configurations, via implementation of a general flow solver into the OpenFOAM code version 20.06. In this paper, Giustini and Issa's (2021) methodology is used as a vehicle to further demonstrate the capabilities of interface capturing CFD and it is applied to model bubble growth due to phase change in conditions of a set of Vol. 5, No. 4, 2023, 357-364 Experimental andComputational Multiphase Flow https://doi.org/10.1007/s42757-022-0139-5 G. Giustini,R. I. Issa 358 laboratory experiments (Florschuetz et al., 1969). Interfaceresolving CFD, and other approaches such the Lattice Boltzmann method (Li et al., 2020), are promising tools for simulating two-phase flows of real fluids with large density ratios and interfacial evaporation. The perspectives of the work are applications to modelling of boiling phenomena, which are integral part of power generation, e.g., water-cooled nuclear reactors (Giustini, 2020), and thermal management technologies (Karayiannis and Mahmoud, 2017). Welch and Wilson (2000) appear to be the first to introduce mechanistic phase change in Interface Capturing (VOF) two-phase flow simulation. The method assumes thermodynamic equilibrium at the vapour-liquid interface, which is taken at the saturation temperature at the prevailing pressure, and models phase change via explicit reconstruction of temperature gradients at the moving phase boundary, a method which has been adopted by other workers for similar applications (Gibou et al., 2007;Sato and Ničeno, 2013). Such methods often require unwieldy procedures to reconstruct the exact location of the phase boundary and compute the local temperature gradients, which can be challenging to implement into general calculation procedures on "unstructured" mesh configurations (Sahut et al., 2021) (e.g., tetrahedral, polyhedral, or hybrid meshes). A different approach relies on a so-called "continuum representation" of the interface, whereby interfacial phenomena such as phase change are accounted for as volumetric source terms located at an implicitly defined phase boundary e.g., a narrow transition region between the bulk phases. The template of such methods, described in Hardt and Wondra (2008), was found to require elaborate measures in order to tackle simulation of flows with large density differences and strong interfacial phase change (Kunkelmann and Stephan, 2009), a shortcoming that was eliminated by Giustini and Issa (2021).
The simulation framework established in Giustini and Issa (2021) is here adopted to conduct simulations of the growth of vapour bubbles in isolation, i.e., away from solid surfaces, in both zero-gravity and normal (Earth) gravity conditions.
Simulations are validated via comparison against experimental observations (Florschuetz et al., 1969). The validation cases represent a severe test for Interface Capturing methodologies due to large density ratios, the presence of strong interfacial evaporation, and bubble rise motion.

Methodology
Simulations presented in this paper build on Giustini and Issa's (2021) methodology, which comprises a two-phase fluid flow solver (Section 2.1.1) coupled to a method for tracking the evolution of the liquid-vapour interface (Section 2.1.2) and to a heat transport model (Section 2.1.3) augmented by a calculation of the rate of phase change (Section 2.1.2.1) at the liquid-vapour interface.

Model formulation
The current methodology is based on the single-fluid formulation (Tryggvason et al., 2011) of two-phase fluid dynamics whereby one single set of momentum balance equations with variable fluid properties is used to describe the instantaneous behaviour of a two-phase flow. An indicator function, in this case the volume fraction of the liquid phase (equal to 1 in the liquid phase and to 0 in the vapour phase), is used to track the phase distribution and to assign the relevant fluid properties to each phase. Interfacial phenomena, such as surface tension forces and mass transfer, are captured in an interfacial region, where the indicator function varies sharply but continuously from the value pertaining to the bulk liquid to the value of the bulk vapour, over a thickness of a few computational cells. Surface tension (Brackbill et al., 1992) and mass transfer (Kunkelmann and Stephan, 2009) are calculated in the interfacial region using suitable implementations of established models and are incorporated into the solution of the numerical model as volumetric source terms. As such, they do not require special treatment, which would otherwise be needed with other approaches (Gibou et al., 2007;Rajkotwala et al., 2019) such as "immersed boundary" techniques. As outlined in Section 2.2, straightforward representation of interfacial terms as volumetric sources enables implementing the current calculation procedure using standard solution methods (Ferziger et al., 2002), typical of single-phase CFD, and the Finite Volume discretization approach, embodied in the OpenFOAM CFD code (Greenshields, 2017).

Fluid flow model
The incompressible Navier-Stokes equations read ( μ being the fluid viscosity), p is the pressure, and the body force g and the surface tension force σ f (the latter to be discussed in Section 2.2.1).
The velocity divergence ⋅ u is linked to the volumetric rate of interfacial phase change 3 (kg / (m s)) m ⋅  by a continuity constraint (Tryggvason et al., 2011): where the subscripts l and v indicate, respectively, the liquid and vapour phases, both being considered to be incompressible.
With the current Interface Capturing methodology (Giustini and Issa, 2021), broadly based on the VOF method (Welch and Wilson, 2000), the phases are distinguished with an indicator function here taken as the volume fraction of liquid α , which is used to compute "mixture" fluid properties, denoted by ρ , as (fluid viscosity μ , specific heat capacity c, and thermal conductivity k are treated in the same way):

Interface capturing method
The Interface Capturing method of Giustini and Issa (2021), uses an advection equation for the indicator function α originally introduced in Olsson and Kreiss (2005) and augmented by terms arising from interfacial phase change: n is a compressive term required to maintain the sharp interfaces, ( ) ε α ⋅  is a numerical diffusion term to control the interface thickness and n is the unit vector normal to the vapour-liquid interface. The parameter ε , which can be varied to adjust the thickness of the interfacial region, is a user specified parameter linked to the local grid size Δx . In this work ε is set equal to Δ / 2 x , which is the recommended value (Olsson and Kreiss, 2005) for vapour-liquid flow simulations. The terms , which represent mass transfer due to phase change (Badillo, 2012), are functions of the volumetric rate of phase change m  , which is computed as follows.

Modelling of phase change
Carey's (2020) formula, originally introduced in Schrage (1953), is used to compute the volumetric rate of interfacial phase change: where T denotes the temperature and SAT T is the value of the saturation temperature at the prevailing system pressure.
Here the value of the "accommodation" coefficient σ is set equal to one for all simulations presented, as recommended in Carey (2020). In the above expression for the evaporation coefficient φ , lv h is the latent heat of vaporization, M is the molar mass, and g R is the universal gas constant. The rate of phase change, which is zero in the bulk phases and non-zero in the interface cells where 0 α  ¹ , is positive for evaporation ( SAT 0 T T -> ) and negative for condensation

Thermal model
Heat transfer is modelled via a transport equation for the fluid specific enthalpy , c being the "mixture" specific heat capacity), here defined as which takes the value of the liquid specific enthalpy Here the reference temperature REF T is taken as the saturation temperature at the prevailing system pressure SAT .
T The thermal transport model equation (Giustini and Issa, 2021) here adopted reads The right-hand side of Eq. (8) represents a sink term due to phase change and is modelled as where n is the unit vector normal to the vapour-liquid interface. The first term in square brackets in Eq. (9) represents latent heat transfer, and the second term accounts for sensible heat transfer to due phase change, e.g., the formation of vapour and removal of liquid in evaporation. The ε parameter takes the same value of Δ / 2 x adopted in the interface advection equation (4).

Model solution
The current flow solver incorporates a standard segregated (Ferziger et al., 2002) solution procedure based on the PISO algorithm (Issa, 1986).
For each control volume (cell) O, the semi-discrete form of the momentum balance equation (1) can be written in compact notation as and o O u denotes the velocity at the previous time step. The Euler scheme is adopted for temporal discretization (as is the case for all simulations presented in this paper), δt being the time step. Here the hydrostatic pressure is lumped together with the static pressure to give the piezometric pressure where ( )  denotes linear interpolation from cell centroid O to face centroid f and f S is the surface area vector normal to face f.
The pressure field is calculated via solution of an equation for z p : where Δ f z p is the piezometric pressure difference across face f and the second term on the right-hand side is a source term due to mass transfer.
The flux is then corrected using the pressure field calculated from iterative solution of Eq. (12) and the velocity at the cell centroid is updated using the corrected flux (Perot, 2000). The updated velocity field is then used to recalculate ( ) H u and repeat the pressure-velocity correction sequence for a desired number of times.

Calculation of the surface tension force
The Continuum Surface Force (CSF) method (Brackbill et al., 1992) is here adopted to compute the surface tension force at the vapour-liquid interface. The CSF method expresses the surface tension force σ f in terms of the local curvature of the interface κ = -⋅ n: where σ is the surface tension coefficient, here taken as a constant.
In standard implementations of the CSF method, normal n and curvature κ = -⋅ n are computed from the indicator function gradient α  via the relationship (Tryggvason et al., 2011): With the approach, it is often the case that numerical errors, implicit in evaluating the interface unit normal from a steeply varying indicator function distribution, introduce considerable distortion of the calculated curvature field. Such a numerical imperfection is the likely root cause of inaccurate estimates of the surface tension force commonly reported in the literature (Popinet, 2018) which usually manifest themselves as "spurious velocities" (sometimes called "parasitic currents"), i.e., unphysical velocity oscillations near the vapour-liquid interface. In order to alleviate the issue, a smoothed gradient  α  of the indicator function is here used in the surface tension calculation. The smoothed gradient  α  is calculated recursively from the "raw" gradient α  computed from the indicator function distribution resulting from solution of Eq. (4). The recursive formula can be expressed as denotes linear interpolation from cell centroid O to face centroid f.
With gradients of the indicator function so computed, no significant velocity oscillations were observed at the vapour-liquid interface, as demonstrated by simulation results presented in Section 3.2.

Sequence of solution steps and OpenFOAM numerical solution settings
For every time step, the following sequence is used to solve the model equations presented in the preceding text: (1) Advance the interface via solving Eq. (4) using a known rate of phase change from the previous time step.
(3) Compute the surface tension force with the CSF method.
(6) Update the rate of phase change.
(7) Compute the velocity and pressure fields via the desired number of PISO steps (as outlined in Section 2.2).
The interface advection equation (4) is solved with the standard MULES solver (Damiàn, 2013) embedded in the OpenFOAM code.
For all simulations presented in this work, the time step was dynamically updated based on the CFL condition, where Δx is the smallest cell characteristic dimension, u is the maximum value of the velocity magnitude in the domain, and CFL S is a safety factor here set equal to 0.1. Typical observed flow velocities are of the order 0.1 m/s, which correspond, for cell sizes as small as 0.78 μm here employed, to time step values around 10 -6 s.
Linear interpolation was used for the divergence terms and for the diffusive terms in the model equations, via selection of the relevant discretization schemes in the OpenFOAM input specifications. Further details on OpenFOAM discretization schemes are available online (Damiàn, 2013;Greenshields, 2017).

Validation: Simulations of free bubble growth
Predictions of the current method are validated via comparison with experimental observations of bubble growth in ethanol from Florschuetz et al.'s (1969) dataset. Firstly, spherical bubble growth in an extended pool of liquid ethanol in zero-gravity conditions is simulated. More challenging normal (Earth) gravity conditions, whereby a bubble grows while rising through stationary liquid, are considered next. The new simulations presented in this paper extend the set of test cases beyond conditions considered in Giustini and Issa (2021). For all the simulations presented in this work, uniform meshes in axisymmetric configuration have been used.

Simulation of spherical bubble growth in zero-gravity conditions
The test case considers bubble growth in a pool of stationary liquid at a uniform temperature T ¥ . The initial condition of the simulation has been calculated with the theory of Scriven (1959), which is applicable to modelling bubble growth due to heat transport from surrounding uniformly superheated liquid in zero gravity. The theory assumes that the vapour in the interior of a spherical growing bubble is at the saturation temperature SAT T at the prevailing system pressure and that heat is transported from the liquid towards the curved surface of the bubble through a thermal boundary layer where the temperature varies from the remote value T ¥ (at some distance into the liquid) to the saturation value SAT T (prevailing in the interior of the bubble). Under these assumptions, Scriven (1959) provides closed-form expressions for the bubble radius as a function of time, and for the temporal variation of the spatial distribution of the temperature around the bubble as a function of the radial distance r, reproduced below: The above expressions are used to initialize a small bubble in the axisymmetric simulation domain as in Fig. 1, where the initial bubble radius 0 R and the indicative temperature distribution are imposed as initial condition of the simulation calculated with Scriven's (1959) theory at a short time 0 t from bubble inception. The subsequent growth of the bubble is captured by the current method, and the time history of the bubble radius is extracted from the simulation and compared with experimental data.
Thermophysical properties of saturated ethanol in atmospheric pressure conditions here used are listed in Table 1 The simulated time history of the bubble radius is compared with experiment in Fig. 2, showing results for increasingly fine uniform grids. Very similar results obtained with the two finest grids corresponding to uniform cell sizes with Δ 1.56μm x = and Δ 0.78μm x = indicate mesh convergence of the current method if sufficiently fine grids are used. The simulation using a grid with Δ 1.56μm x = has been repeated on a larger domain to accommodate the growing bubble at later time into the simulation and enable comparison with experimental data up to approximately 180 ms into bubble growth, corresponding to the final point in the experimental time series, as shown in Fig. 3. Very good agreement with experimental data is observed.

Growth of a rising bubble in normal gravity
The test case considers the growth of an ethanol bubble rising in a uniformly superheated pool of liquid in conditions of Florschuetz et al.'s (1969) experiments.  Normally in CFD simulations of rising bubbles, the continuous liquid phase is kept stationary and the bubble rises through it. The approach may imply the necessity to use a very tall computational domain in the vertical direction to simulate the entire bubble travel time. In order to eliminate the issue, here the computational frame of reference is attached to the bubble, which is kept stationary while the surrounding liquid is moved instead. The corresponding simulation setup is illustrated in Fig. 4 ① . The velocity imposed as a uniform inlet boundary condition at the top end of the domain, directed downwards along the gravitational acceleration vector, is recalculated at every time step to keep the center of gravity of the bubble stationary, in the same position occupied by it at the beginning of the simulation.
The remote liquid superheat is equal to 3.1 K and the initial conditions of the simulation have been determined as outlined in the preceding text, using the analytical solution for spherical bubble growth, which has been used to initialize a small bubble and the surrounding temperature distribution at the center of the computational domain. In the region of the computational domain where the bubble develops, a uniform mesh with Δ 3.12μm x = is used. The initial bubble radius is equal to 208 μm. Figure 5 shows the temperature distribution around the growing bubble at different times into the simulation. Details of the temperature and velocity distribution at 85 ms into the simulation are shown in Fig. 6, showing a wake of cold liquid downstream of the bubble and a thin thermal boundary layer around the portion of the bubble-curved surface exposed to the incoming liquid. Figure 7 compares the equivalent bubble radius extracted     (Florschuetz et al., 1969) time histories of the rising bubble radius. from simulation with two experimental runs from Florschuetz et al.'s (1969) experimental dataset, indicating good agreement between modelling and measurement.

Conclusions
This paper presented a numerical methodology for simulating mass transfer at phasic interfaces and its application to modelling of bubble growth in conditions of an experimental dataset available in the literature. The method, which is applicable to arbitrary fluids, was implemented into the general purpose CFD code OpenFOAM and validated via simulations of vapour bubbles growing in a stationary pool of liquid away from solid boundaries. The chosen validation cases represent a challenge for CFD modelling of boiling flows due to large density ratios and strong interfacial evaporation. Comparison with experimental data indicated good agreement between modelled and measured bubble growth in zero-gravity and normal (Earth) gravity conditions.