Numerical simulation of bubble growth in a supersaturated solution

In this paper, a Volume of Fluid (VOF) based approach to simulate the growth of a pre-existing bubble in a supersaturated solution is developed and implemented in OpenFOAM R ©. The model incorporates the Compressive Continuous Species Transfer approach to describe the transport of dissolved gas and surface tension is treated using the Sharp Surface Force method. The driving force for bubble growth is defined using Fick’s 1 st law and a Sherwood number based correlation. The source terms for the governing equations are implemented by extending the work by Hardt and Wondra, J. Comp. Phys. 227 (2008) 5871–5895. The predictions of the proposed solver is compared against theoretical models for bubble growth in supersaturated solutions. The effect of spurious currents, which are generated while modelling surface tension, on bubble growth is also investigated. The proposed approach is used to model the growth of a rising bubble in the supersaturated solution. © 2020 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The growth of gas bubbles in supersaturated solutions occurs when the amount of dissolved gas ( C ) is higher than the saturation concentration ( C sat ) in the liquid given by Henry's law [1] . Some processes where this phenomena is relevant are effervescence in beer [2] , champagne [3] and carbonated beverages [4] , decompression sickness [5] , and electrochemical systems like electrolysis of water [6] and chloralkaline processes [7] .
Bubble evolution in supersaturated solutions consists of nucleation of the bubble, which is an atomistic event, whereas bubble dynamics and the associated convection are continuum scale events [8] . Different modes by which nucleation can occur in a supersaturated solution has been proposed by Jones et al. [9] , Vachaparambil and Einarsrud [10] and the critical radius ( R c ) required is where C / C sat indicates the supersaturation ratio ( S ), σ is the surface tension and P is the operating pressure [11] . The current paper delves into the continuum scale phenomena of bubble evolution in supersaturated systems, but not atomistic scale events. The works by [8,[11][12][13] provides a brief overview of various aspects of research in atomistic scale events relevant in bubble evolution.
The literature that has investigated bubble growth in supersaturated solutions can be divided into theoretical, experimental and numerical approaches. The two landmark theoretical papers are those by Epstein-Plesset [14] and Scriven [15] . Epstein and Plesset [14] studied the evolution of the bubble radius at various supersaturation levels when the effect of convection associated with the bubble growth is neglected. Scriven [15] considered the effect of this convection and obtained an asymptotic expression for the evolution of bubble radius. To account for the growth of bubble from a pre-existing bubble, Scriven's asymptotic solution is extended in [16] . On the other hand, experimental studies has reported the bubble radius to evolve as R = At x where A and x are growth and time coefficients, respectively [15,[17][18][19] . A value of x approximately equal to 0.5, derived by Scriven [15] , indicates that the rate determining step is the diffusion of dissolved gas across the interface [17] . When rate determining step for bubble growth is the electrochemical reactions, x has been reported to be equal to 1/3 [18,19] . The effect of bubble growth on rise velocities of bubbles is studied in works like [20] .
Computational Fluid Dynamics (CFD) is a widely used numerical approach to model continuum scale processes like bubble evolution. The main multiphase approaches used in CFD are Euler-Euler (EE), Euler-Lagrange (EL) and Volume of Fluid (VOF) methods. Although EE [21][22][23][24][25][26] and EL [8,25,[27][28][29][30][31] approaches have been used to model gas evolution in supersaturated systems, they rely on the so called closure laws to model the interaction between the phases [32] . Compared to EE and EL approaches, the VOF method advects the volume fractions of the phases without using any approximation for the shape of the interface [32] , which provides a versatile method to study the bubble evolution. In the work by Liu et al. [33] , evolution of a single hydrogen bubble from the electrode surface has been modelled using the VOF approach and validated by comparison against evolution of bubble radius obtained from experiments. Other studies have used the VOF approach to model large carbon dioxide bubbles but the smaller bubble (which are below the mesh resolution) are described using a population balance model [34][35][36] and EL method [37,38] which are then validated using experimental data on electrical potential and bubble thickness respectively. Interestingly, these previous works [33][34][35][36][37][38] are developed based on a commercial solvers like FLUENT R . The VOF based approaches typically use a Sherwood number ( Sh ) based correlation to describe the mass transfer across the interface [33][34][35][36] . Modelling mass transfer based on Sh can lead to erroneous results as these correlations are dependent on both the bubble shape and local velocity field [39,40] . A more universal approach is to use the Fick's 1 st law of diffusive mass transfer instead [5,41,42] .
Although considerable progress has been made in stimulating of bubble evolution driven by supersaturation, certain aspects of the phenomena such as bubble dynamics, mass transfer and momentum exchange between the phases are usually modelled rather than resolved to reduce the computational costs and are often suitable for very specific applications. To the best knowledge of the authors, there is a lack of generic solvers that can model bubble dynamics in supersaturated solutions and are validated using theoretical benchmarks. To further improve the computational approaches used, insights from the simulating boiling can be relevant. A review of the various computational approaches used for simulating boiling and condensation phenomena is available in [43] . One such approach, reported in the study by [44] , developed a VOF based model for evaporation by using a continuum-field representation of source terms in FLUENT R which was subsequently implemented in OpenFOAM R 1.5.x [45] as 'evapVOFHardt' [46,47] . Apart from evapVOFHardt, other approaches developed using OpenFOAM R to model phase change driven by temperature are [4 8 , 4 9] and 'phaseChangeHeatFoam' [50,51] .
In this paper, we propose a VOF based solver to model the bubble evolution in a supersaturated solution. The proposed solver is based on interFOAM (available in OpenFOAM R 6 [52] ) and the mass transfer approach, developed for evaporation by Hardt and Wondra [44] as well as Kunkelmann and co-workers [46,47] , is extended and implemented for the bubble growth driven by supersaturation of dissolved gases. The driving force for the bubble growth is modelled using both the Fick's 1 st law and a correlation using Sh . The transport of the dissolved gas is modelled using a Compressive Continuous Species Transfer (C-CST) model [53] and the surface tension force is modelled using the Sharp Surface Force approach [54] . Following the derivation of the governing equations and its implementation in OpenFOAM R , the mass transfer model in the proposed solver is validated against theoretical models like the approximate solution of Epstein-Plesset [14] , asymptotic solution of Scriven [15] and the extended Scriven model [16] . Spurious currents, produced due to modelling surface tension, and its effects on bubble growth are investigated. Finally, the proposed model is used to simulate the growth of a rising bubble in a supersaturated solution.

Governing equations and its implementation
The VOF method uses a scalar function to represent the interface and individual phases as In the proposed solver, α 1 represents the volume fraction of Phase 1 (or liquid) and the volume fraction of Phase 2 (or bubble) can be written as α 2 = 1 − α 1 .
Fluid properties like density ( ρ) and viscosity ( μ) are determined using the volume fraction as a weighting function: The continuity equation for a growing bubble is written as where ˙ m = ˙ m 1 α 1 + ˙ m 2 α 2 , with ˙ m 1 and ˙ m 2 representing source/sink terms for Phase 1 and Phase 2, are defined in Eq. (24) . The continuity equation ( Eq. (4) ) can be used to derive a corresponding equation to describe the volume fraction as shown below Based on Eqs. (4) and (5) , the evolution of α 1 can be computed as which is analogous to the volume fraction equation described in [46,47] . In order to render a sharp interface, the advection term in Eq. (6) is modified to include a compressive flux as where U r is compressive velocity defined in Eq. (8) .
where φ, is volume flux that is computed based on velocity and S f is cell face area vector, C α and δ are velocity flux, cell surface area, a user-defined compression factor and a small non-zero value which is calculated as 10 −8 / respectively. C α controls the interface smearing and is typically set to a value between zero and four [55] . In order to efficiently compute the pressure boundary condition and density jump at the interface, the solver uses a modified pressure ( p rgh ), defined as which is explained in [56] . Due to the use of modified pressure, −∇ p + ρ g can be written as −∇ p + ρ g = −∇ p rgh − g · x ∇ ρ. (10) The momentum equation, which is solved to obtain the velocity field in both the phases, is calculated as where last term represents the surface tension force which in this solver is computed based on the Sharp Surface Force (SSF) model [54] . Substituting Eq. (10) in Eq. (11) gives the momentum equation as implemented in OpenFOAM R , The SSF formulation used to describe the surface tension force can be defined as where κ final is curvature of the interface obtained using a three step procedure described in [54] and α sh is a sharpened volume fraction of liquid defined as where C sh is a sharpening coefficient. If C sh = 0, α sh is equivalent to α 1 whereas when C sh = 1 describes a very sharp interface that is numerically unstable [54] . Based on preliminary simulations, SSF model is preferred over the commonly used Continuum Surface Force model [57] because of its ability to reliably simulate a sub-millimeter bubble.
As the bubble growth is driven by supersaturation, the dissolved gas concentration ( C i ) can be represented as Eq. (15) represents a concentration of dissolved gas that causes bubble growth when C i > 0 as it represents concentration greater than C sat . Assuming that the interface of a bubble growing is saturated, the interface is represented by C i equal to zero, based on Eq. (15) . The governing equation for the transport of C i is based on the Compressive Continous Species Transfer (C-CST) model, proposed by Maes and Soulaine [53] , as where B is equal to (1 − He i ) / (α 1 + α 2 He i ) , S i is the sink term for the dissolved gas in the liquid at the interface (described in Eq. (25) ), U r is the compressive velocity (defined in Eq. (8) ) and ˆ D i is the harmonic interpolation of the diffusion coefficients [53] : where D i ,1 is the diffusion coefficient of dissolved gas in Phase 1 and D i ,2 is the self-diffusion coefficient of the gas [58] . In Eq. (16) , He i represents dimensionless Henry's constant that accounts for the jump in Ci across the interface that is expressed as C i, 2 = He i C i, 1 . In order to obtain the saturation condition of the interface when C i is defined based on Eq. (15) , He i should be equal to zero. The concentration of the dissolved gas ( C i ) represents the α The mass transfer across the interface is described using two approaches: Fick's 1 st law and Sherwood number based correlation. Fick's 1 st law, which is applicable for any arbitrary flow scenario, is implemented as where M i is the molar mass of the species i . The alternative method, which is typically implemented in practical applications, utilizes a Sherwood number correlation which can be written as where k is the mean mass transfer coefficient defined as k = D i, 1 Sh/ 2 R, where the mean Sherwood number ( Sh ) is derived from Fick's 1 st law for the mass transfer at interface for a spherical bubble in creeping flow [39] as In Eq. (20) , The local mass transfer rate ( ψ 0 ) can be written based on j (defined in Eqs. (18) and (19) ) as where N is a normalization factor is defined in Eq. (26) . The term α 1 is multiplied Eq. (22) to ensure that the mass flux that drives the bubble growth is based only on the supersaturation of the liquid.
In order to increase numerical stability, ψ 0 is smeared over a few computational cells by solving where ψ describes the smeared ψ 0 at every time step [44] . In Eq. (23) , D t is a product of a diffusion constant and parameter with dimension of time which controls the amount of smearing in Eq. (23) [44] . The source term of the continuity equation ( Eq. (4 )), ˙ m = ˙ m 1 α 1 + ˙ m 2 α 2 reduces to ˙ m = ˙ m 2 α 2 because liquid is not consumed when during bubble growth due to supersaturation. Based on the value of ψ, source term ˙ m (in Eq. (4 )) can be calculated as ˙ m = Aα 2 ψ, (24) where A is a normalization factor defined in Eq. (26) . The source term ˙ m , which is defined in Eq. (24) , is calculated in the region where α 1 < 0.001 by artificially moving them away from the interface based on volume fraction of the liquid (further described in [44,47] ). The sink term in dissolved gas transport ( S i ), in Eq. (16) , which represents the dissolved gas lost into the bubble due to mass transfer across the interface is calculated as where N is the normalization factor defined in Eq. (26) and α 1 in Eq. (25) ensures that S i is non-zero only at the liquid side of the interface. The normalization factors used in the solver can be defined as where is the domain for the flow computation. The overall algorithm used to solve the governing equations ( Eqs. (4) , (7), (12), (16) and (23) ) is summarized in Fig. 1 . The volume fraction equation ( Eq. (7) ) is solved the Multidimensional Universal Limiter with Explicit Solution (MULES) method [57] . The Pressure Implicit with Splitting of Operator (PISO) algorithm is used to solve the momentum ( Eq. (12) ) and continuity ( Eq. (4) ) equations [56,57,59] . The PISO algorithm recasts the continuity equation into a pressure Poisson equation (or 'pressure correction equation') which is solved and then used to update the predicted velocity fields [56,57] .

Computational setup
The computational domain is a square of dimensions 3 cm × 3 cm with a pre-existing bubble of radius equal to 250 μm located at the center of the geometry. The boundaries are set to zero gradient conditions for U , ψ, C i and α 1 while p rgh is equal to 101,325 Pa. The phenomena modelled in this paper corresponds to a pre-existing bubble of radius equal to 250 μm growing in a solution supersaturated by carbon dioxide. The saturation concentration of dissolved carbon dioxide in water is calculated using Henry's law at 25 • C and 101,325 Pa to be 33.44 mol/m 3 [60] . The initial bulk concentration of carbon dioxide in water is set to C i = 200 . 64 mol/m 3 , corresponding to a supersaturation ratio (defined as C / C sat ) equal to seven. To validate the mass transfer model implemented in the proposed solver, some of the simulations neglect the gravity and surface tension. When gravity and surface tension are treated in selected simulations, g is assigned as (0 −9.81 0) m/s 2 and σ equal to 0.0468 N/m [31] , respectively.
The settings used in the simulations are described in Table 1 and internal fields used for the simulations is tabulated in Table 2 . Values for D t and He i are set equal to 10 −6 m 2 and 10 −4 respectively and the maximum Courant number ( Co max ) is set equal to 0.05. In simulations where surface tension is neglected, the maximum time step is allowed to adjust  Table 2 Initial internal field of parameters used for the simulations.   automatically based on the Co max . When surface tension is considered, the time step constraint, which is required to prevent the temporal growth of spurious currents, is calculated as where τ μ and τ ρ are time scales which are defined as μ avg x / σ and ρ a v g ( x ) 3 /σ respectively, μ avg and ρ avg are defined as the average dynamic viscosity and density between the phases whereas x is the mesh resolution which is equal to 7.5 μm [57] . The values of C 1 and C 2 has been reported to be equal to 0.01 and 10, respectively [57] . For the M4 mesh used in the simulations, the constraint on the time step is calculated, based on Eq. (27) , to be equal to 7.2 ×10 −7 s. A parametric study to the investigate the effect of pre-existing bubble size, D t and He i on the solution is discussed later in the paper. The governing equations ( Eqs. (4) , (7), (12), (16) and (23) ) are discretised using schemes as mentioned in Table 3 with the relevant settings summarized in Tables 4 and 5 . The tolerence criteria while solving the pressure correction equation is set to 10 −20 to reduce the force imbalance (between surface tension and pressure forces) which is generated due to the iterative procedure used in the solution algorithm [57] . The maximum number of iterations while solving the governing equations are set such that the tolerance criteria set in Table 4 is met at every time step.

Convergence studies
The convergence of the simulations is studied using the Fick's 1 st law as the driving force (in Eq. (18) ) for bubble growth and neglecting gravity. The surface tension is also neglected in these simulations as the spurious currents generated are Table 5 Other parameters used in solving the discretised equations.
Used for interface compression in Eq. (8) .
momentumPredictor no Controls solving of the momentum predictor; typically set to 'no' for multiphase and low Reynolds number flows [70] . nOuterCorrectors 1 PISO algorithm is selected by setting this parameter equal to unity in PIMPLE algorithm [70] .
The number of times the PISO algorithm solves the pressure and momentum equation in each step; usually set to 2 or 3 [70] .
Specifies the under-relaxation factors; set equal to one for transient simulations [54] . C sh 0.3 Sharpening coefficient in Eq. (14) ; to model a reliable pre-existing bubble of radius 250 μm [74] .
dependent on the mesh resolution and a 'true' mesh convergence is difficult to achieve especially for capillary dominant flows [54,61,62] . The convergence is analyzed based on three criteria: independence of solution from the mesh, conservation of phases and dissolved gas, effect of domain size of the solution, and monitoring supersaturation/pressure data away from the bubble. The grid convergence is analyzed using six different hexahedral meshes, as shown in Table 6 , based on the radius normalized growth rates at t = 4 . 5 s. The growth rates and bubble radii for various meshes, shown in Fig. 2 , at t = 4 . 5 s show a relative deviation of nearly 2% between the M4 and M5 meshes. This discrepancy is due to the discontinuous initial condition of C i (in Table 2 ) which leads to a larger initial growth rate, offsetting the bubble radius for further time steps, see Fig. 2 (a). The radius normalized growth rate, equivalent to the mass flux in 2D, is not influenced by this offset. As the relative change of the radius normalized growth rate at t = 4 . 5 s between two finest meshes is lower than 0.25%, see Table 6 , M4 has been used for simulations in the paper.
The imbalance in the phases, illustrated in Fig. 3 , is lower than 0.1% of the volume of the liquid initially in the system. Calculated based on the amount of dissolved gas which is numerically in the bubble, the imbalance in the dissolved gas is lower than 0.1% of initial amount of dissolved gas in the system as shown in in Fig. 4 .
The region far stream from the bubble growth should have a constant concentration and pressure as it is not effect by the bubble growth. Monitoring a point away from the interface shows that C i and p rgh remains equal to 200.64 mol/m 3 and 101,325 Pa, respectively for the duration of the simulation.
Using smaller domain for the simulations has been observed to affect the radially symmetric nature of the velocity field as shown in Fig. 5 . In order reduce the effect of this error, the domain used in the simulations is nearly twelve times the bubble diameter at t = 10s i.e. 3 cm × 3 cm, which has been showed in Fig. 5 (c).

Validation
In this section, the proposed solver is validated by comparing against two theoretical cases i.e. the approximate solution of Epstein-Plesset [14] , the asymptotic solution of Scriven [15] and Extended Scriven [16] . It is worth pointing out that surface tension is neglected in the simulations discussed in this section unless explicitly specified. The simulations used for comparison can be divided based on the description of driving force that drives the mass transfer across the interface as • Constant driving force: The definition of j based on Sh correlation ( Eq. (19) ) is modified to account for driving force for bubble growth that is dependent on the bulk concentration of the supersaturated solution ( C ∞ ) which is expressed as • Local driving force: Eq. (19) provides a Sh correlation based definition of driving force that is dependent on the local concentration of the dissolved gas around the interface.

When convection of species concentration is neglected
The simulations consider a two dimensional growing bubble, a rate of bubble that is growing can be described as where R is the radius of the bubble and j is the constant driving force that causes the bubble growth (described in Eq. (28) ), i.e.
When the convection caused by bubble growth is neglected, i.e. Re = 0 , the above equation can be written as Integrating the above equation from R 0 at t = 0 s to R at t gives which is the 'approximate solution of Epstein-Plesset' [14] when the effect of surface tension and convection is neglected. Eq. (32) provides bubble growth when the bubble growth is driven by a the bulk concentration of the dissolved gas.
As the effect of convection on the bubble growth is neglected in the theoretical benchmark, the solver is modified by implementing the mass transfer coefficient as Using the mass transfer across the interface that is governed by Eqs. (33) and (28) in the proposed solver provides a prediction of bubble growth that is equivalent to the Epstein-Plesset solution as shown in Fig. 6 . The predictions from the proposed solver marginally under-predicts the final radius of the bubble by less than 0.01%. The concentration of the dissolved gas, interface and velocity distribution in the domain for bubble growth driven by a constant driving force (using mass transfer coefficient described in Eq. (33) ) is illustrated in Fig. 7 .

When convection of species concentration is considered
The asymptotic solution of bubble growth in a supersaturated solution, proposed by Scriven [15] , when convection is considered, is given by where β is the growth parameter, valid only for diffusion controlled growth of spherical bubble in an unbounded medium.
Based on the work by [16] , an 'Extended Scriven model' that modifies Eq. (34) to account for the growth from a pre-existing bubble of radius R 0 can be described as Previous works by Wang and others [4,9,19] , reported an analytical expression to determine β for a spherical bubble as which has been shown to agree with the predictions by Scriven [15] . The formulations proposed in [4,9,15,19] to calculate β is derived for a spherical bubbles which is different from the 2D bubbles simulated in the current paper. For direct comparison to the simulations, β 2 D , derived in Appendix A based on the work by Wang et al. [19] , can be calculated as where a is defined the same as Eq. (36) . For the phenomena modelled in this paper, β 3 D and β 2 D can be calculated, using Eqs. (36) and (37) , are equal to 5.3346 and 4.0509, respectively. The velocity, concentration of dissolved gas and interface position is compared for Fick's 1 st law ( Eq. (18) ) and local driving force ( Eqs. (19) and (21) ) in Fig. 8 . The temporal evolution of bubble size from simulations are compared to the Scriven's asymptotic ( Eq. (34) ) and Extended Scriven ( Eq. (35) ) in Fig. 9 . As expected, the bubble growth predicted by the Scriven's asymptotic expression and Extended Scriven theories using β 2 D provides a better representation than β 3 D .
The driving force based on the Sh correlation under predicts the bubble radius by nearly a factor of 2 at t = 10 s while Fick's 1 st law provides a better agreement to the theoretical predictions using β 2 D , with an error less than 2.5% at t = 10 s.
The discrepancy between the simulations can be explained using the growth rate in Fig. 10 , where the growth rate is nearly 6.5 times smaller for the Sh correlation. Although the growth rate predicted by the driving force described by Fick's 1 st law is initially different from the corresponding theory, it seems to asymptotically match the theoretical prediction. The initial discrepancy between the two is due to the discontinuous nature of concentration while initializing, see Table 2 .

Influence of surface tension
Due to small length scales associated with bubble growth, surface tension dominates the flow physics. Modelling surface tension has its own challenges, namely the generation of spurious currents around the interface [44,50,54,57] . Due to the time step constraint, described in Eq. (27) , the simulations are run only until 2 ms.
The spurious currents ( U sc ), calculated as max( | U | ), are observed on both sides of the interface as illustrated in Fig. 11 and its temporal variation plotted in Fig. 12 (a). The convection that is generated by spurious currents seems to remove the dissolved gas at the interface which results in the reduction of the growth rate (calculated as volume integral of ψ 0 in the computational domain) observed in Fig. 12 (c). Due to the generation of large magnitude of spurious currents initially, there is a substantially drop in the growth rate of the bubble in the first few time steps as seen in Fig. 12 (c). The average growth rate of the bubble during 2 ms reduces by approximately 32% compared to when surface tension is not treated. The reduction in the growth rate is also reflected in the temporal variation of bubble radius, as shown in Fig. 12 (b). At t = 0 . 002 s, the bubble radius differs by nearly 0.15% from the bubble radius when surface tension is not treated.
where p 0 is the operating pressure is equal to 101,325 Pa. For the two-dimensional bubble, the Laplace pressure in the bubble is calculated based on the Young-Laplace equation ( p * c = σ /R ). The error associated with predicting the Laplace pressure can be calculated based on ( p c − p * c ) / p * c and its temporal variation is plotted in Fig. 13 . Although there is large initial discrepancy, the absolute error quickly reduces to values lower than 0.1, in Fig. 13 .

Growth of a rising bubble in a supersaturated solution
The computational domain for modelling the growth and rise of a two-dimensional bubble is 3 mm × 9 mm. The preexisting bubble of radius equal to 250 μm is initialized such that its center is at a distance of 1mm from the bottom boundary and equidistant from the side boundaries. The dissolved carbon dioxide in the water is set at 200.64 mol/m 3 . Both gravity and surface tension are treated in this simulation. The four boundaries are described using a zero gradient conditions for α 1 , C i , U and ψ whereas a Dirchlet condition (equal to 101,325 Pa) is used for p rgh . The mesh used in the simulation has a grid resolution ( x ) and R 0 / x equal to 7.5 μm and 33.33, corresponding to the M4 mesh in Table 6 . The simulation is terminated at 0.04s but the maximum time step, calculated based on Eq. (27) , and maximum Courant number permitted are set to 7.2 ×10 −7 s and 0.05 respectively. The driving force for bubble growth is described using Fick's 1 st law as it provides a more realistic growth rate than Sh based correlation as shown in Fig. 9 . The imbalance in the phases and dissolved gas in the simulation of the growth of the rising bubble is lower than 0.1% of the amount present initially in the system.
The changes in the bubble morphology along with its position in the computational domain as it evolves is illustrated in Fig. 14 . The distribution of the dissolved gas around the bubble and velocity in the domain due to the rising bubble at t = 0.04s is illustrated in Fig. 15 (a) and (b), respectively. As the growing bubble grows and rises, the dissolved gas at rear of the rising bubble gets depleted before the mass transfer by convection and diffusion can replenish it. On the other hand, the incoming supersaturated liquid always replenishes the depleted dissolved gas. This variation in the concentration of the Fig. 9. Comparison between the simulated bubble radii using different driving forces ( Sh based Local driving force and Fick's 1 st law) and theoretical models ( Eqs. (34) and (35) ). dissolved gas around the interface leads to larger local mass transfer rate, calculated based on Eq. (22) , in front of the bubble in comparison to its rear as shown in Fig. 15 (c).
In order to understand the effect of bubble growth on rising, its rise velocity is compared to a bubble that is just rising without any growth. The rise velocity of the bubble is computed as the bubble volume averaged vertical component of the velocity vector [54] . The bubble rising without any growth is implemented by setting the dissolved gas concentration to zero in the simulation. Fig. 16 (a) shows that, for 0.04s simulated, there is no substantial change in the rise velocity due to growth of the bubble by mass transfer across the interface. The corresponding growth of the rising bubble is illustrated through the increase in area of the bubble with time in Fig. 16 (b). The insignificant change in the rise velocity indicates that the change in buoyancy force experienced by the bubble does not change as a result of the growth, which agrees with a previous experimental study [20] .

Influence of user defined parameters
In this section, the effect of user defined parameters, like He i , D t and size of pre-existing bubble, are investigated using the the setup described in the validation studies i.e. neglecting surface tension and gravity while using Fick's 1 st law as driving force for bubble growth.
Smearing of ψ 0 to obtain ψ, using Eq. (23 ), relies on a user defined D t and ideally the solution should be independent of the effect of this parameter. As ψ is used to compute the source term required for bubble growth (i.e. ˙ m which is defined  The parameter He i , which is used to model concentration jump across the interface, should theoretically be set equal to zero to model the transport of the dissolved gas and describe the interface as saturated. When He i is set equal to zero the denominator in B (in Eq. (16 )) becomes infinity for α 1 = 0 . Although He i cannot be set equal to zero, the saturated condition of the interface can be reasonably reproduced by using a low enough value of He i . So He i is set to a non-zero number  (i.e. 10 −4 , 10 −3 , 10 −2 , 10 −1 and one) and its influence on the concentration distribution of dissolved gas and associated bubble growth is compared in Figs. 18 and 19 , respectively. When He i is equal to 10 −3 or 10 −4 , both bubble growth and concentration of dissolved gas in bubble becomes nearly independent of He i which indicates that the saturation condition of the interface has been nearly reproduced. It is also worth pointing out that setting He i to 10 −3 or 10 −4 , introduces some temporal unboundedness in the beginning of the simulation in the value of C i . As the mass balance of the dissolved gas is not affected due to this temporal unboundedness of C i , we consider effect of this error to be rather negligible.
The effect of the size of pre-existing bubbles (radii equal to 250 μ m and 500 μm) on the bubble growth is investigated in Fig. 20 . The computational domain and the mesh resolution used in these simulations has been changed to make the   Fig. 20 , shows that smaller the bubble higher the growth rate.
Although size of the pre-existing bubble (250 μm) used in the simulations is much larger than the critical radius, calculated using Eq. (1) , it is in the same order as the radii of the pre-existing bubbles used in the theoretical work like [14] but larger than the cavity size reported in experiments, which is typically around 50-200μm [4] . These micrometer sized bubble are very difficult to model due the presence of spurious currents that can sometimes be large enough enough to render simulations inaccurate [63] . As lower limit to the size of the pre-existing bubble is dictated by the spurious currents, implementation of more advanced surface tension models are required to model micrometer size bubbles. It is also worth pointing out that boiling studies performed using OpenFOAM R [46,64] and condensation studies [50] typically use pre-existing bubbles radius in the same range as the ones used in this paper but other solvers, which use more advanced interface reconstruction algorithms, like Piecewise-Linear Interface Calculation (PLIC) scheme [44] , and/or surface tension approach, using height functions [65] , enables modelling even smaller pre-existing bubbles (in the order of few micrometers) due to the lower spurious current [44,66] .

Conclusions
In this paper, a new VOF based solver to model bubble evolution in supersaturated systems is implemented in Open-FOAM 6. The proposed solved is created by adding dissolved gas transport equation (C-CST model [53] ), surface tension (SSF [54] ) model, driving force for bubble growth (Fick's 1 st law and a Sh based correlation) and the relevant source terms are implemented, by extending the work of [44,46,47] , in interFOAM.
The mass transfer models applied in the proposed solver are validated based on theoretical models like Epstein-Plesset [14] which do not consider the effect of convection, and Scriven [15] as well as Extended Scriven [16] which accounts for the effect of radial bubble growth. The proposed solver utilizing a driving force based on Eqs. (28) and (33) , matches well with the approximate solution of Epstein-Plesset [14] as shown in Fig. 6 . A driving force based on Fick's 1 st law provides a better agreement to Scriven [15] and Extended Scriven [16] models than a Sh correlation as shown in Figs. 9 and 10 . The modelling of the growth of the rising bubble using Fick's 1 st law shows that the proposed solver, for the duration of 0.04 s, predicted the negligible effect of the bubble growth on rise velocity. Further remarks about the proposed solver are: • Spurious currents, generated while modelling surface tension, introduces numerical convection near the interface that reduces the growth rate of the bubble by advecting away the dissolved gas at the interface. • As expected, the solver is able to predict higher growth rate of smaller bubbles compared to larger bubbles and the increase in local mass transfer rate at the front of the bubble than its rear when a bubble is rising. • The simulations also shows that the bubble growth is sensitive to the value of D t but the solution becomes nearly independent of the parameter at larger values. • The use of He i equal to 10 −4 in the C-CST has been shown to describe the transport of dissolved gas as well as the saturation condition at the interface reasonably well.
Based on the derivation of expression to determine β 3 D , i.e. Eq. (36) which is described in [19] , we use the same approach to obtain a formulation for β 2 D that describes a 2D bubble growing in a supersaturated solution. The growth rate of the bubble can be calculated as where first term describes the rate of increase of the bubble volume, h describes a unit grid thickness which is set equal to 10 −6 m which is used in the simulations, and J 1 + J 2 corresponds to the sum of mass transfer across the interface ( Eq. (A.3) ) and effect of convection generated by bubble growth ( Eq. (A.4b) ) [19] . The theortical distribution of dissolved gas ( C ), once concentration boundary layer is developed, can be expressed as a sigmoid function, based on simulations as shown in Fig. 8 (c) that satisfies the following boundary conditions r = R, C = C sat , where m is a constant that is defined in Eq. (A.6) . The diffusive mass transfer across the interface can be described as 3) The effect of convection established by the increase in bubble radius on the growth rate is treated as where A is the interface surface area that is equal to 2 π Rh which can be used to reduce the above equation as When J 2 is neglected, the above equation must reduce to Eq. (31) , which shows that m = 2 . (A. 6) As bubble radius in 2D also evolve as described by Eqs. (34) and (35) , bubble radius evolution can be described as R =