{'fraction_non_alphanumeric': 0.027316018580735634, 'fraction_numerical': 0.009991527651991001, 'mean_word_length': 4.988453463960812, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A method for simulating interfacial mass transfer on arbitrary meshes

This paper presents a method for modelling interfacial mass transfer in Interface Capturing simulations of two-phase flow with phase change. The model enables mechanistic prediction of the local rate of phase change at the vapour-liquid interface on arbitrary computational meshes 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. The proposed methodology for capturing mass transfer is applicable to arbitrary meshes without the need to introduce elaborate but artificial smearing of the mass transfer term as is often done in other techniques. The method has been validated via comparison with analytical solutions for planar interface evaporation and bubble growth test cases, and against experimental observations of steam bubble growth.


INTRODUCTION
Computation of two-phase flows with phase change requires combined solution of the hydrodynamics and heat transfer processes in the liquid and vapour phases, plus some form of tracking of the interface motion and of the local rate of interfacial mass transfer. A mathematical model that is general enough to be applied to any phase change process does not exist yet due to the multi-physics and multi-scale nature of the problem. Nevertheless, the development of Interface Capturing modelling techniques [1], such as the Volume Of Fluid (VOF) and Level Set methods, at the continuum scale enabled significant progress to be made in microscopic modelling of phase change phenomena in two-phase flow, at least in laboratory conditions. The vast majority of methods [2][3][4] for modelling mass transfer have been implemented in research codes applicable only to structured meshes and elementary geometries, and the only method [5] so far conceived to be applicable to arbitrary meshes and complex geometries typical of 'industrial' CFD is affected by serious shortcomings.
Here we present an Interface Capturing framework for simulation of multiphase flows with phase change on arbitrary meshes. The method's cornerstone is a straightforward representation of interfacial mass transfer that is applicable to any vapour-liquid phase change phenomenon and arbitrary fluids and is amenable to simulation on arbitrary meshes. For demonstration, the methodology was implemented using the open-source CFD toolbox OpenFOAM and was verified with a set of test cases that entails simulation of real fluids with large density difference and strong interfacial evaporation. The aim of this article is to demonstrate that the proposed method retains the capabilities of state-of-the-art methods for capturing mass transfer while extending their range of applicability to arbitrary meshes and eliminating the elaborate measures introduced by the Ref [5] methodology.
Ref [6] is considered to be the first to apply the VOF method to simulate two-phase flow with interfacial mass transfer using mechanistic evaluation of the rate of phase change at the vapourliquid interface. Similar VOF-based methods were later adopted by a number of workers and applied to model nucleate boiling [7][8][9] [10], film boiling [5] and boiling in microchannels [11,12]. The most popular VOF implementation, based on the work of Ref [5], has been implemented in the Computational Fluid Dynamics (CFD) codes ANSYS [12] and OpenFOAM [8], and is therefore in principle applicable to arbitrary geometries and computational meshes typical of industrial CFD simulation. However, that methodology presents serious shortcomings due to the introduction of numerical smearing of the interfacial rate of phase change, which is some distance removed from the true interfacial region, being artificially relocated in order to achieve numerical stability. Despite shortcomings introduced by such a spurious numerical 'cropping' of the physical rate of phase-change, which the current method eliminates, the work of Ref [5] is so far the only Interface Capturing method for phase change that has been implemented in industrial CFD codes. The Level-Set [3,13] [14] and Front Tracking [2,15] methods were also applied to various phase change problems, however these methods have not yet been used in industrial CFD practice and their application has been limited to idealised test cases; they also suffer from problems in conserving volume and mass.
A different approach, based on the Constrained Interpolation Profile (CIP) method [16], has been implemented in a general-purpose two-phase flow solver [4] [17], which has been extensively used for modelling boiling phenomena [18]. The accuracy of the methodology was demonstrated both via comparison with analytical solutions and experimental data [19] and that method is therefore considered as a benchmark for numerical simulation of phase change in two-phase flow. Despite these superior attributes, the method is still not applicable on arbitrary meshes. Likewise, in Ref [20] a capable Phase Field model for quantitative prediction of mass transfer phenomena, albeit only applicable to elementary geometries, was developed.
In all the above methods, a common approach to modelling phase-change at the vapour-liquid interface is based on explicit evaluation of the heat flux jump ('Stefan condition' [21] [1]) at the phase boundary, which requires explicit reconstruction of the interface and of the interfacial temperature gradients, typically resulting in unwieldy procedures [3,4,7] in order to extrapolate the temperature gradient from/to each side of the interface. An alternative approach has been derived from first principles within the framework of near-equilibrium thermodynamics and applied to Phase Field simulations of boiling phenomena [20]. The treatment of mass transfer in Ref [20] is based on the principles of equilibrium thermodynamics and was used to determine consistent formulations for the rate of phase change (as a function of the local deviations of interface temperature from its saturation value) and for the interfacial heat balance equation in near equilibrium conditions. An appropriate version of the general approach, suitable for incorporation into a VOF-based methodology, is that of Ref [22], in which the rate of phase change depends explicitly on the interface thickness, which cannot be justified on physical grounds. Different values of the interface thickness are used in the two implementations of the method reported so far [22] [20]; however, neither author provided justification for any particular imposed value of the interface thickness. Use of other methods, such as Lattice Boltzmann simulation [23] [24], has so far been limited to idealised fluids and test cases; also narrow is the range of applicability of conventional unstructured-mesh CFD approaches, which have typically been developed as application-specific tools [25,26].
The majority of the Interface Capturing methods discussed above are embodied in research codes of limited industrial applicability; their implementation in the generally applicable framework of 'industrial' CFD for simulation on arbitrary geometries and computational meshes is still in the early stages of development. It is therefore highly desirable to generalise the phase-change modelling capabilities of current CFD approaches and enable their application to the simulation of two-phase flow on arbitrary meshes. The perspectives of such an enabling technology are applications to modelling boiling phenomena in industrial applications, requiring an ability to model mass transfer in geometrically complex flow configurations [27,28] and to solve accurately the boiling model equations for real fluids on irregular grids, which are outside the current capabilities of CFD.
In this paper, we present the development of a methodology to capture mass transfer in Interface Capturing simulations of two-phase flows on arbitrary grids, using the CFD toolbox OpenFOAM [29]. The proposed simulation methodology is a novel approach that builds on the strengths of current general modelling techniques and extends their applicability to arbitrary meshes through implementation in a flow solver based on the VOF Interface Capturing method.
The methodology is suitable for modelling boiling and bubble growth in industrially relevant fluids (e.g. water boiling at atmospheric pressure), which represent a real challenge for most Interface Capturing numerical methodologies, due to the large density difference between the liquid and vapour phases and the presence of strong evaporation. Therefore, a set of planar surface evaporation and bubble growth test cases has been chosen to test the proposed methodology. The objective of this paper is to demonstrate the accuracy of the current simulation methodology, via comparison with analytical solutions and laboratory observations of steam bubble growth in realistic boiling conditions.

MATHEMATICAL MODEL
The proposed simulation framework comprises a fluid flow model (Section 2.1), a method for interface capturing (Section 2.2) and a thermal model (Section 2.4). These are augmented in order to simulate, at the vapour-liquid interface, the effect of mass transfer due to phase change.
The rate of phase change is computed mechanistically at the interface based on the predicted local temperature distribution (Section 2.3).

Fluid flow model
The incompressible Navier-Stokes equations can be written as where is the density, the fluid velocity, = (2).
In equation (2) and in the following, the subscript l indicates the liquid phase and v the vapour phase.
The current methodology is based on the Volume Of Fluid (VOF) method [1], whereby the phases are identified with the volume fraction of liquid , which is used to compute a generic (or mixture) fluid property, say the density , as .
The phase indicator varies from 1 (in the liquid) to 0 (in the vapour) across the fluid interface, whereby the phase indicator (in this case the liquid volume fraction) is used to model a transition [5] between the bulk phases, as opposed to other methods ( [3,4] [30]) that are based on reconstructing explicitly the geometry of the phase boundary.

Interface capturing model
It is assumed that the vapour and liquid phases are separated by a sharp interface. The interface is advanced through time using the model of Ref [20], which is applicable to conditions where thermally driven phase change occurs. The model equation for a generic conserved phase indicator derived in Ref [20] is used here to advance in time the liquid volume fraction distribution: The last term in the equation is operative only in the interface region which is inevitably diffused numerically over a few mesh cells and it arises from the consideration that the liquid and vapour phases have different velocities within that region, as derived in Ref [20].
Using the continuity relation (2), the ∇ • term can be eliminated from the right-hand side of equation (4), which can be expressed in terms of the rate of phase change ̇ and rewritten as .
Ordinary interpolation schemes fail to maintain a constant interface thickness when used to discretize an advection equation such as (5) and typically cause unphysical diffusion of the interfacial region. The problem of interface diffusion is well known in computational science and a popular method to alleviate it has been proposed by Ref [31], who introduced a modified advection equation, augmented by additional compressive terms for maintaining the sharp interfaces. Equation (5) is therefore augmented to include interface 'sharpening' terms |∇/| is the interface unit normal vector and is the desired interface thickness, which needs to be specified as an input to the model. Proof of the effect of the additional terms in equation (6) was illustrated from first principles in the original work [31] in a Level-Set framework (in the absence of phase change); the Ref [31] sharpening technique was also successfully applied to flows with phase change in Ref [4], using the Constrained Interpolation Profile [16] method for advancing the interface. In this work, as in the previous studies that used the Ref [31] sharpening method, the interface thickness is related to the grid size. The interface transport model here proposed has been solved with the standard MULES solver available in the OpenFOAM code [32].

Calculation of the rate of phase change
The rate of phase change is computed at the interface from the local temperature distribution.
Various approaches can be used to this end.
A straightforward method [5] is based on considerations at the molecular level [33] and uses a basic kinetic theory description of the interfacial phase change process to link local deviations of the interface temperature from its saturation value to the rate of phase change pertaining to the continuum fluid model. With the kinetic approach, the rate of phase change can be computed in the interfacial region as where 34' is taken as the saturation temperature at the externally imposed pressure and the coefficient is computed as where ℎ -. is the latent heat of vaporization, M the molar mass and ( the universal gas constant. The assumption implicit in equation (7) is that the portion of surface area of vapour-liquid interface contained in one cell can be expressed as <=' = |∇ | >?--. ̇≠ 0 only in the interfacial region where |∇ | ≠ 0 and ̇= 0 elsewhere; the rate of phase change is positive for the case of evaporation ( > 34' ) and negative for the case of condensation ( < 34' ). A theory for predicting the 'accommodation' [34] coefficient M does not yet exist, and therefore the appropriate value for M needs to be extracted from experiment. Various experimental sources indicate a value close to 1, as summarised in a review [34] that covers experimental works on evaporation and condensation [35] spanning several decades. Here a value of 1 is retained for the coefficient M. It was found in [5,11] that flow solvers using equation (7) for modelling mass transfer were subject to numerical instabilities in the case of strong evaporation and large density ratios -/ . , as is the case in most flows of practical interest. In order to alleviate the issue, a 'cropping' procedure was introduced, whereby the source terms due to mass transfer are artificially removed from the interfacial region and located some distance away. The current solver is immune to such issues and it has therefore been possible to eliminate the unphysical 'cropping' step.
For all simulations presented in this paper, the volumetric rate of phase change ̇ is computed with equation (7).

Thermal model
In most flows of practical interest, interfacial mass transfer is driven by diffusion of heat towards or away from the vapour-liquid interface; in typical boiling conditions, for example, evaporation is driven by the flow of heat from the liquid towards the interface. Therefore, it is necessary to solve an appropriate enthalpy balance equation in order to obtain the temperature field required for evaluation of the rate of phase change. To this end, a single-fluid enthalpy is where c is the specific heat capacity and a possible choice of the reference temperature is the saturation temperature at the prevailing system pressure, that is, Enthalpy transport is modelled with the following balance equation is the mixture thermal diffusivity.
The volumetric heat sink ; is due to phase change and is modelled after Refs [20,22] as where, as noted in the derivation presented in Ref [20], is the same quantity appearing also in the interface advection equation (6). Equation (10) can be solved with any standard iterative method. In the current simulations, use of a pre-conditioned bi-conjugate gradient solver, available in OpenFOAM [29], was made. Note that owing to the current formulation of mass transfer it is possible to model the effect of phase change on the energy balance with a volumetric sink term, which enables one to dispense with explicit reconstruction of the interfacial heat flux jump.

Advecting the interfacial region
The interface is advanced through time using the standard MULES solver available in the OpenFOAM code. In semi-discrete form, equation (6) can be written for a control volume O as Which is solved explicitly using the MULES solver In the above equations, K is the control volume, f indicates the surfaces separating the control where (∇ N ) M is the component of the volume fraction gradient normal to surface f.

Pressure-velocity coupling
Coupling between pressure and velocity is achieved with a PISO-based [36] [37] segregated algorithm.
With O as the control volume under consideration and i indicating neighbouring control volumes, the advection-diffusion matrix is unpacked into its diagonal and off-diagonal components: , and the details of the particular interpolation schemes used to generate the coefficients P are specified in a later section. The central coefficient K is augmented by a contribution K Q = I -J# due to the temporal derivative. The form of the central coefficient A so computed, which includes K Q , depends on the chosen temporal discretization scheme, also to be specified in a later section.
As is common in modelling buoyant flows, the hydrostatic pressure is lumped together with the static pressure to give the piezometric pressure as a main working variable R = − • with the body force being computed as ( = −( • )∇ , x being the position vector.
In the adopted segregated solution procedure, the nodal velocity K is advanced in time and corrected as follows: First, the component of the velocity normal to face f (or flux) is estimated as M = ℒ 6 ( ) with the face density ̅ M taken as the harmonic average [38] of densities in the two adjacent cells: where K and P interpolation factors based on the distance of each cell centre i or O from the face centre f. For a uniform grid, the face density would be Here f indicates the surface separating the control volume O from each one of its neighbours i, ℒ indicates linear interpolation from cell centres to cell faces, M is the face area vector (normal to face f). The surface tension force ) is computed from the interface curvature using the wellestablished Continuum Surface Force (CSF) method of Ref [39] as where is the surface tension coefficient and the curvature is computed as The pressure field is computed by solving a modified equation that includes the source term ̇6 $ % " − $ % # 7 due to mass transfer, that is: where ∆ M R is the piezometric pressure difference across face f. The flux is then corrected using the newly calculated pressure field, and used to update the nodal velocity as n being an integer used to count corrector steps.
The corrected velocity can be used to recalculate ( ) and repeat the corrector step. The reconstruction operator [40] [41] is defined as ℛl M n = B∑ . ⨂ .
3 . M  For unsteady simulations, two corrector steps are sufficient for convergence within a time step [36,42]. The sequence of steps for enforcing coupling of the pressure and velocity fields highlighted in the preceding text is in the same spirit as that which may be employed by any standard segregated flow solver; the modifications introduced here are for extending the solution procedure to two-phase flows (equation 16) and for including the effect of mass transfer (equation 21).

Solution loop
The solution procedure entails solving the Navier-Stokes equations, including the volumetric CSF approximation of the surface tension force and the body force, the energy and interface advection equations, and evaluating the mass transfer model.
For every time step, the segregated solution sequence is as follows: 1. Evaluate equation (13) to advance the interface using a known rate of phase change from the previous time step.
3. Compute the surface tension force with the CSF method.
4. Compute the rate of phase change with equation (7).  Table 1. Unless specified otherwise in Table 1, linear interpolation has been used for divergence terms and the Gauss theorem for gradient terms. Further details on OpenFOAM discretization schemes are available online [29] and in the Ref [32] thesis.
Term Scheme 'Interface compression' [32] ∇ • ( ℎ) Van Leer Table 1 Discretization schemes employed. If not specified in the table, linear interpolation and the Gauss theorem have been used for discretization of, respectively, divergence and gradient terms.

VALIDATION -PLANAR INTERFACE TEST CASES
In order to validate the current method, simple one-dimensional planar interface test cases are considered first.
In the Stefan problem test (Section 4.1), the evaporation-driven motion of a flat interface separating two fluids with equal density is considered; for this test, owing to the equal density assumption, there is no fluid dilatation due to evaporation, the domain boundaries are impermeable and the interface moves solely due to the local effect of mass transfer.
In the sucking interface problem (Section 4.2), another one-dimensional flat-interface benchmark, saturated water at atmospheric pressure is the test fluid, with a density ratio

Stefan problem
The Stefan problem concerns the prediction of the evaporation-driven motion of a planar interface due to transport of heat from superheated vapour. .
The fluid properties listed in Table 2 have been used to evaluate the above expressions. For uniform grids considered here and in the following, =   Table 2 Physical properties of the fictitious fluid considered for the Stefan problem.
As shown in Figure 2, the time history of the interface position matches the analytical solution, confirming the ability of the current methodology to capture interface motion due to interfacial mass transfer for the case of fluids with equal density.

Sucking interface problem
The physical problem considered here is that of predicting the motion of a planar interface, and The corresponding temperature distribution in the liquid is For evaluation of the analytical solutions, the thermophysical properties of saturated water at 1 bar, as indicated in Table 3, have been used. i is set equal to 383.15K.

Property
Vapour Liquid  Table 3 Physical properties of saturated water at 1 bar.
In the simulation, the initial temperature distribution and interface position correspond to the analytical solution after 0.1 seconds; the simulation is continued for 1.1 seconds.
Mesh-convergent behaviour is observed in Figure 4

VALIDATION -ISOLATED BUBBLE GROWTH IN A POOL OF LIQUID
A set of isolated bubble growth problems is used to validate predictions of mass transfer at a curved interface.
The physical problem is that of predicting the growth of a spherical bubble in a uniformly superheated pool of liquid, away from solid surfaces and in the absence of gravity, as originally nevertheless, sodium is an important substance widely used in the energy conversion industry.
The main difference between sodium and fluids, such as water or refrigerants, more commonly studied using CFD is that liquid sodium has a much higher thermal conductivity than water or any refrigerant.
Using a Scriven test for the case of water, the current method is validated on irregular polyhedral grids in Section 5.1.3. Such a qualitative test, using low-resolution grids to limit the computational effort, is used to demonstrate applicability of the method to arbitrary meshes.
Finally, the current method is validated in Section 5.2 and Section 6 against laboratory observations [44] [45] of steam bubble growth at a distance from walls in a uniformly heated pool of liquid.
In the test cases that follow, axisymmetric simulation has been used to compare results with spherical analytical solutions and experimental observations.

Analytical solution
The growth of an isolated bubble in a uniformly superheated pool of liquid is one of the few boiling test cases for which a solution for the bubble radius as a function of time can be derived analytically [43]. The theory 1 assumes that a bubble grows due to the diffusion of heat from the remote liquid towards the curved surface of the bubble [46]. Under the assumption that the vapour pressure inside the bubble is the same as the remote liquid pressure (i.e. the externally imposed system pressure), it is assumed that the body of vapour in the bubble is at the saturation temperature corresponding to the externally imposed pressure, and that saturated vapour is generated at the bubble surface. Neglecting the effect of gravity and surface tension, the growth of the bubble is limited purely by the rate of diffusion of heat towards the bubble surface at 34' , from the remote liquid assumed at some higher temperature i . Scriven 1 The theory is predicated on the assumption that bubble growth is driven purely by thermal effects. Ref [37] showed that assuming thermally-driven bubble growth is appropriate at times after bubble inception larger than a characteristic time scale computed as = ( ) = 2 t -(bubble radius) (31).
The growth constant is evaluated via solution of the following transcendental equation (a discussion of the applicability of the Scriven approach to computing the growth constant of bubbles is presented in [47]): Once the growth constant and the temporal variation of bubble radius are known, the radial temperature distribution in the liquid (as noted, the vapour is assumed at the saturation temperature) can be computed at any time t via evaluation of the following expression: (temperature distribution in the liquid) (33). , and the latent heat transfer due to evaporation can be expressed as and -w# vv must be equal and therefore ' ≅ The above set of expression can be used to compute bubble radius, temperature distribution and thermal boundary layer thickness at any time and is used in the following to compute the initial conditions of the simulations presented.

Modelled bubble growth vs analytical solution -water
Fluid properties are those of saturated water at atmospheric pressure indicated in Table 3.
Initial conditions for the simulation are generated with the Scriven model equations (described in the preceding text) at Q = 0.17 into bubble growth (sufficiently large compared to ≈ 7 × 10 6h ). This is an arbitrary time chosen in order to obtain initial radius Q and initial  This test case is more demanding, from the point of view of numerical resolution, than the planar interface case (Section 4.2), as the grid size needs to be fine enough to resolve the very thin initial thermal boundary layer on the liquid side of the interface, which is initially about 5 microns thick. In practice it was found that only a few mesh points are required to resolve the initial thermal boundary layer, as indicated in Figure 6, where good agreement with analytical solution is observed for the two finest meshes that correspond, respectively, to 3 and 6 grid points in the initial thermal boundary layer. The spherical bubble tests were affected to some extent by the build-up of spurious currents, a well-known numerical artefact discussed, for example, in Refs [4] [38] [48]. The currents are due to irreducible errors inherent (according to Ref [39]) in the surface tension model, however they were of small enough magnitude for the numerical solutions not to be affected significantly.

Figure 6
Comparison with analytical solution of modelled time history of bubble radius for the case of bubble growth in isolation in a uniformly superheated pool of liquid water.

Simulation using an irregular polyhedral mesh
The above test case was repeated using an irregular two-dimensional translationally symmetric grid, consisting of polyhedral cells of average size equal to approximately ∆ ----≈ 4.5 , shown in Figure 7, where the initial volume fraction distribution is also indicated. For this test, = ∆ ----/2 and the results are compared with the case using a regular two-dimensional translationally symmetric grid with ∆ = 4.8 .   Figure 9. A qualitative comparison of time histories of the bubble radius modelled using the polyhedral and uniform quadrilateral grids is shown in Figure 10 and it demonstrates that the current method retains on irregular grids approximately the same accuracy attained using uniform grids.

Figure 8
Volume fraction distribution at a sequence of times for bubble growth simulation using a polyhedral mesh.

Figure 9
Magnified view of the volume fraction distribution in the interfacial region for bubble growth simulation on a polyhedral mesh.

Figure 10
Comparison between time histories of the bubble radius modelled using a uniform translationally symmetric grid ('2D', dashed line) and a non-uniform polyhedral translationally symmetric grid ('2D', solid line).

Modelled bubble growth vs analytical solution -sodium
Comparison with analytical solution was repeated for the case of a sodium vapour bubble.
Liquid sodium is used alongside water as coolant in nuclear reactors [49] and therefore, as for the case of water, it is important to study its behaviour in conditions when nuclear reactions generate enough heat to cause the fluid to boil; the fluid is also used in single-phase liquid and phase-change receivers of solar power plants [50]. The properties of saturated sodium [51,52] in conditions typical of its industrial applications are listed in ), corresponding to an initial radius Q of 2 mm and initial liquid thermal boundary layer thickness of 0.386 mm. As for the previous case of water, agreement with the analytical solution is observed as the grid is increasingly refined, as shown in Figure   11.

Figure 11
Comparison with analytical solution of modelled time history of bubble radius for the case of bubble growth in isolation in a uniformly superheated pool of liquid sodium.

Comparison with experiment
The experiments of Ref [44] consist of optical observations of steam bubbles growing in isolation in a uniformly superheated pool of water. These are very unusual measurements because normally steam bubbles are observed growing at a heated wall, which introduces the Comparison between numerical and experimental bubble growth is shown in Figure 12, showing all three experimental runs from the Ref [44] data set, for the case of 3.1K of liquid superheat, corresponding to a bubble growth constant = 9.4. Numerical values of the bubble radius follow closely the experimental data are always well within the bounds of the experimental uncertainty of ±10% estimated in Ref [44].

Figure 12
Comparison between numerical (red solid line) and experimental (black/white symbols) bubble growth.

RISING STEAM BUBBLE UNDER GRAVITY
The final validation test considers the simultaneous growth and rise of a steam bubble in an initially stationary pool of superheated water under gravity. Experimental data for comparison are taken from Ref [45], reporting three series of bubble growth observations at similar superheat levels of 3.2, 3.5 and 3.6 K. The axisymmetric simulation domain is a rectangle of 6.0 due to mass transfer and rises in the pool of liquid due to buoyancy, as indicated in Figure 13, showing close-up views of temperature and velocity distributions, for a uniform discretization with ∆ = 5.0 . The fields are displayed at different times until 0.048 s into the simulation i.e., the latest time at which an experimental value for the bubble radius is available. Drag due to relative motion between the rising bubble and the surrounding stationary liquid causes the bubble shape to depart from a sphere, while currents generated by bubble rise under the effect of gravity distort the temperature distribution around the bubble, which loses its spherical symmetry, and generate a wake of cold liquid at positions swept by the bubble's vertical motion. Comparison between modelled and measured bubble growth is indicated in Figure 14.
Values of the bubble radius have been divided by the group 2 tfor comparison with the Ref [45] experimental data, which are presented in such a form. The comparison increases confidence in applying the current method to situations where interface motion is not solely due to phase change. External forces such as gravity and hence non-trivial velocity fields causing convective heat transport are captured adequately, at least in conditions of interest and to the extent where available experimental data on bubble growth provide a meaningful comparison.

Figure 14
Comparison between modelled and measured bubble radius for the Ref [45] rising steam bubble test case.

CONCLUSIONS
This paper presented development of a numerical method to compute interfacial mass transfer for Interface Capturing simulation of two-phase flow. Implementation of the current methodology using the standard interpolation and solution techniques embodied in the OpenFOAM code enabled its application to arbitrary computational meshes, such as polyhedral meshes here employed that are typical of industrial CFD practice. The accuracy of the method, applicable to arbitrary fluids, was tested via application to evaporation and bubble growth problems representative of the boiling of real fluids, in this case water and sodium. The chosen set of problems represents a severe test for Interface Capturing methodologies due to large density ratios and the presence of strong interfacial evaporation; the proposed framework is however not limited to boiling and may be extended to other phase change phenomena.
Comparison with laboratory observations of steam bubbles growing in an extended pool of uniformly superheated water demonstrated that the proposed methodology is applicable to modelling phase-change phenomena in realistic conditions.

This work was supported by the United Kingdom Engineering and Physical Sciences Research
Council (EPSRC) through grant EP/T027061/1.