Laser-induced vapour bubble as a means for crystal nucleation in supersaturated solutions—Formulation of a numerical framework

We use in this work numerical simulations to investigate the evolution of a laser-induced vapour bubble with a special focus on the resolution of a thin layer of liquid around the bubble. The application of interest is laser-induced crystallization, where the bubble acts as a nucleation site for crystals. Experimental results indicate the extreme dynamics of these bubbles where the interface during the period of 200 us, from nucleation to collapse, reaches a maximum radius of roughly 700 µm and attains a velocity of well above 20 m/s. To fully resolve the dynamics of the bubble, the volume of fluid (VOF) numerical framework is used. Inertia, thermal effects, and phase-change phenomena are identified as the governing phenomena for the bubble dynamics. We develop and implement into our numerical framework an interface phase-change model that takes into account both evaporation and condensation. The performed simulations produce qualitatively promising results that are in fair agreement with both experiments and analytical solutions from the literature. The reasons behind the observed differences are discussed and suggestions are made for future improvements of the framework.


Introduction
Crystals of proteins, salts, and other chemicals are needed in a variety of fields, from chemists and biologists analyzing the crystal atomic structure using X-rays to industrial processes that rely on controlling crystal formation at the correct time and place in the manufacturing process for drugs and other useful compounds (Tatalovic, 2009). In the production of such compounds, crystallization is a common separation or purification process and has been used for a long time in chemical, pharmaceutical, food, and material industries. The main goal of the crystallization process is to produce a set of crystals with desired and controlled properties, such as crystal size, morphology, or purity. However, controlling the properties of the crystals during the crystallisation procedure is not trivial and is a subject of much ongoing research and development. A big concern is to control where, when, and at what rate crystals are formed in a solution in order to achieve their desired properties. At present, there are a number of potentially interesting technologies, for example, sonocrystallization (Ruecroft et al., 2005) and non-photochemical laser-induced nucleation (Garetz et al., 1996). Among these, laser-induced cavitation is a promising example in order to achieve good control of the crystallisation process (Nakamura et al., 2007;Yoshikawa et al., 2014).
When a focused laser beam is applied to a supersaturated solution, a vapour bubble is formed that rapidly increases in size due to evaporation of the superheated liquid. The bubble expands until the superheated liquid is evaporated, at which point the vapour starts to condense back to liquid form, due to heat losses to the surroundings, and the bubble collapses. These laser-induced cavities have been observed to induce crystal nucleation in supersaturated solutions, but the mechanism behind the nucleation is not entirely clear. Most suggestions as to why the crystals are formed are based on assumptions of physical properties in the liquid around the cavitation bubble (Ruecroft et al., 2005;Nakamura et al., 2007;Iefuji et al., 2011;Knott et al., 2011;Yoshikawa et al., 2014;Mirsaleh-Kohan et al., 2017).
One hypothesis about the reason for the crystal nucleation is that there is an increased concentration of solute and a simultaneous cooling of the liquid at the bubble interface, due to evaporation, that leads to crystal nucleation (Soare, 2014). The superheated liquid near the bubble interface is Vol. 1, No. 4, 2019, 242-254 Experimental and Computational Multiphase Flow https://doi.org/10.1007/s42757-019-0024-z Laser-induced vapour bubble as a means for crystal nucleation in supersaturated solutions-Formulation of a numerical framework 243 cooled by evaporation and by the surrounding liquid. The evaporation of the solvent increases the concentration of the solute near the interface and the cooling lowers the liquid solubility. The combination of those two effects, and the fact the crystals have a higher probability to nucleate with increased liquid saturation, may be the prerequisites for crystal nucleation just after the bubble starts to grow. To test this hypothesis, a better understanding of the actual conditions in the vicinity of the bubble is important, but the small scales and fast dynamics of the problem make it very difficult to monitor using conventional measuring techniques. Experimental results by Soare et al. (2011), carried out using an aqueous solution of an inorganic salt, (NH 4 ) 2 SO 4 , give an appreciation of the extreme dynamics of a laserinduced vapour bubble. During the bubble lifetime of about 200 μs, from nucleation to collapse, the bubble reaches a maximum radius of roughly 700 μm and attains an interface velocity of well above 20 m/s. Crystals are observed in a ring around the bubble about 1 s after the laser pulse, but, in the same location as the crystals, small optical disturbances are visible already after 30 μs. Since crystal nucleation occurs at such small scales, it is not possible to observe the process visually, and the nucleation is thought to occur well before the optical disturbances are visible. In such a case, the nuclei have grown large enough to generate the optical disturbances. In the experimental study, it is suggested that the nucleation takes place at the point of maximum rate of vaporization, just after the start of the bubble formation during the laser pulse. Therefore it is of special interest to investigate this early phase of bubble growth.
Numerical simulations that are able to resolve all relevant physical phenomena may provide detailed information about the conditions in the thin layer of liquid around the bubble during its early growth period. According to the generalized Rayleigh-Plesset equation (Rayleigh, 1917;Plesset, 1949), the governing factors for the growth of a vapour bubble in a superheated liquid are initially inertial effects and, after some critical time, thermal effects (Brennen, 1995). During the initial growth phase of the bubble, the growth rate is limited by the inertia of the surrounding liquid. After this period, the growth rate becomes limited by the mass-transfer rate across the bubble-liquid interface which is in turn restricted by the ability of the liquid to transport heat to the interface. In the case of the bubble described in Soare et al. (2011), the critical time that determines the transition from the inertia-controlled to the thermalcontrolled regime is estimated as ~3 ns. Therefore, it seems that the thermal effects are governing the bubble growth rate already after a few nanoseconds and those effects need to be considered in numerical simulations. Zein et al. (2013) and Magaletti et al. (2015) studied the collapse phase of spherical nanobubbles using diffusive interface numerical methods that account for the effects of phase change and obtained results that were in good agreement with experiments. It is, however, not certain how these models are able to capture the growth phase of a bubble and, especially, how well the conditions in the liquid around the bubble can be determined using the diffusive interface approach. Sagar et al. (2018) used a sharp interface method in combination with a cavitation model to model the collapse phase of a laser-induced cavitation bubble near a solid boundary. The obtained results were in favourable agreement to experiments, but the cavitation model that they used neglects the effects of surface tension and inertia, which are important in the case of a small and rapidly expanding bubble that is studied in the present work. Koch et al. (2016) also used a sharp interface approach to model laser-generated bubbles. Their model included the effects of inertia and was used to simulate both the growth and collapse phases of the bubbles with results in good agreement with experiments. The dynamics of the bubbles in their simulations were, however, considered dominated by inertia and compressibility effect and phase change was neglected.
In this work, we use a sharp interface numerical approach that takes into account all the phenomena identified above for the dynamics of a laser-induced vapour bubble. For that purpose, the volume of fluid (VOF) numerical framework is used in combination with an interface mass transfer model outlined in Section 2. The interface mass transfer model is based on a special distribution of source terms in the continuity and energy conservation equations close to the interface. The interface mass transfer rate is computed using the model by Tanasawa (1991), and the source terms are distributed based on the method by Hardt et al. (2008) and Kunkelmann (2011). The aim of the simulations is to capture the dynamics of a laser-induced vapour bubble with the focus on its early growth phase and the conditions of the liquid around the bubble. To simplify the problem, pure water is used for both the liquid and the vapour phases and therefore neither solute concentrations, nor crystal nucleation, are taken into account. Validation results are presented in Section 3 and the results from a laser-induced vapour bubble simulation in Section 4. Finally, in Section 5, the results are discussed and conclusions are drawn.

Computational method
The volume of fluid numerical framework is used to treat the two-phase problem. To account for mass transfer between the phases, a continuity equation for each phase is expressed as l l l l l where C is the color (phase indicator) function, ρ is the fluid density, u is the fluid velocity, l represents the liquid phase, v the vapour phase, and S C is the mass source term. The momentum and energy equations are applied to the combined phases and are defined, respectively, as where the pressure field is denoted as p, μ is the kinematic viscosity, σ is the surface tension coefficient, κ is the interface curvature, E is the energy per unit mass, and S h is the energy source term. The mass and energy transfer across the liquid-vapour interface is accounted for by using a special distribution of the mass and energy source terms in a narrow region close to the interface. Solving the momentum conservation equation in a combined form for both phases makes a momentum source term due to phase change unnecessary because the recoil pressure is the only manifestation of the phase change and that is a direct result of the source terms in the continuity equation (Kunkelmann, 2011;Kharangate et al., 2017).

Phase change model
To determine the values of the mass and energy source terms close to the interface, a model that can accurately estimate the local mass transfer rate is needed. Two common approaches used for this purpose in CFD simulations are the Energy jump condition (Gibou et al., 2007) and variations of the Schrage model (Schrage, 1953). In the Energy jump condition, the saturation temperature is assumed at the interface, and the mass transfer rate is based on the net energy flux across the interface. This method requires a non-trivial computation of the temperature gradients on either side of the interface. The Schrage model is based on the Hertz-Knudsen equation (Knudsen, 1934) and makes use of the Kinetic Theory of Gases to relate the flux of molecules crossing the interface to the temperature and pressure of the phases. Tanasawa (1991) further simplified the Schrage model by assuming that the vapour temperature T v and the interfacial liquid temperature T i satisfy (T v -T i )/T v << 1 and that the corresponding saturation pressures P v and P i fulfil (P v -P i )/P v >> (T v -T i )/ (2T v ). Those assumptions are valid for water vapour except when the vapour temperature and the corresponding saturation pressure are very low. In the cases considered in this study, the vapour temperature is at standard conditions or higher. Therefore, and due to its relatively simple implementation, the Tanasawa model is chosen in this work. There are, however, assumptions of saturation conditions in both phases and non-significant effects of higher vapour phase pressure, made in the derivation of the Schrage model that makes the applicability of that model uncertain in the case of a laser-induced vapour bubble. Nonetheless, the Tanasawa model has been successfully used in somewhat similar cases such as droplet evaporation (Hardt et al., 2008) and vapour bubble growth in a superheated liquid (Magnini et al., 2011). The Tanasawa model is used to compute the interfacial mass flux according to where j is the mass flux due to phase change at the interface, T is the local temperature, T sat is the saturation temperature determined by the local pressure and the Clausius-Clapeyron relation, L is the latent heat, M is the molar mass of the fluid, R is the universal gas constant, and χ is the phase change coefficient. The phase change coefficient χ is used to define the fraction of molecules that change phase and transfer across the interface and 1-χ is the fraction of molecules that is reflected back. In general, different phase change coefficients may be defined for the cases of evaporation and condensation but these coefficients are usually, and in this work also, considered equal (Kharangate et al., 2017). The major uncertainty in the above relation is that the value of χ may depend on both material and flow properties and has to be obtained from experiments. When implementing the Tanasawa model at a local cell level, it has to be ensured that the same amount of mass disappearing on one side of the interface reappears on the other side in the form of the other phase. In addition, the liquid should be cooled due to the absorption of energy during evaporation and the vapour heated due to the release of energy during condensation (latent heat). Hardt et al. (2008) proposed a method that accomplishes this using a special distribution of both mass and energy source terms. First, the evaporation rate on the liquid side of the interface is computed, and then the condensation rate on the vapour side is determined in a similar way. The initial evaporation rate field is defined as where 0 φ is in the form of generation rate per volume and N is a normalization factor representing the ratio between the total interface area to the respective phase part of the interface. The latter factor is determined from the equation: where Ω denotes the full computational domain. By using the liquid volume fraction C l , the evaporation rate is only non-zero on the liquid side of the interface. Also, since the gradient of the volume fraction field is non-zero only at the interface, 0 φ will be localized in a very narrow region at the liquid side of the interface. The normalization factor ensures that the correct interface area is used. To avoid numerical instabilities by having potentially large source terms at the interface, the initial evaporation rate field is smeared using a diffusion equation expressed as where D is a diffusion coefficient. This smearing procedure has to be done for each time step in the simulation and its aim is to smear the evaporation field a certain number a of computational cells. An imaginary diffusion time T * is introduced, representing the duration of the diffusion that smears 0 φ . This parameter, together with the diffusion coefficient D, determines the length scale over which the original source term 0 φ is smeared and it is given by (DT * ) 1/2 . A stability criterion for a conservative explicit Euler discretization in 2D is that the time step size Δt * must fulfil: which in the case of a uniform mesh dx = dy yields: The length scale of the smearing should be in the order of a few computational cells, a, given by To achieve an efficient computation of the smeared source term field, the number of time steps, n T t * * = /D , should be minimized while still fulfilling the before mentioned constraints. Solving this minimization problem gives a value of 2 4 n a = time steps and a diffusion constant of for an arbitrary value of T * . In practice, a value of n slightly higher than the minimum is found to give stable and accurate results. The new, smeared evaporation rate field is obtained as 1 0 ( ) φ φ T * = .
To avoid potential instabilities in the interface region by adding source terms, Kunkelmann (2011) implemented a cropping step where the source term of 1 φ is set to zero in all cells that do not contain pure liquid or pure vapour ( l ). In the simulations, a value of cut 0 001 C = . is chosen as a cut-off limit. To conserve the global evaporation rate of 0 φ and to ensure global mass conservation, the remaining evaporation rate field has to be rescaled. This is done by relating the volume integral of 0 φ to 1 φ on each side of the interface by introducing the scaling coefficients N l and N g as where tot m  represents the total phase change rate in the whole domain and H is the Heaviside function. The final evaporation rate field, 2 φ , can then be expressed as where the first term accounts for creation of mass on the vapour side and the second term for the removal of liquid on the liquid side. An illustration of the different steps in this method is given in Fig. 1. To include the effect of condensation of vapour into liquid, the same phase change model as described above is used but with Cv instead of C l in all of the equations. It is the final evaporation and condensation rate fields that are introduced as mass source terms in the continuity equations (1) and (2) respectively as The energy source term field needs to account for two things considering the phase change process. First, it should include a cooling or heating term that reflects the latent heat of phase change coupled to the mass source terms. Second, when mass is added or removed in a computational cell heat also needs to be added or removed in order to conserve energy. The contribution from the latent heat of phase change can be determined from so that the fluid at the interface is cooled or heated at a rate corresponding to the mass transfer rate. The other contribution, introduced as an energy correction term where mass is added or removed, is defined as where p c is the specific heat capacity of the fluid in the cell and T the temperature. Both contributions are added to the source term in the energy equation as

Validation
To validate the phase change model, two validation cases are performed. The first case, simulating the growth of a vapour bubble, is chosen since an analytical solution by Scriven (1959) is available for comparison and validation. The second validation case is performed to qualitatively assess the ability of the phase change model to handle both evaporation and condensation during the growth and collapse of a vapour bubble.

Validation case 1: Bubble growth in the superheated liquid
The validation case is designed to evaluate the accuracy of the phase change model for the growth of a spherical vapour bubble suspended in an infinitely extended and uniformly superheated liquid. The initial bubble radius is 0.1 mm, and the simulations are performed on a uniform 2D axisymmetric mesh, shown schematically in Fig. 2. An axisymmetric boundary condition is used at the bottom and a symmetry boundary condition at the left boundary, so that only a quarter of the cross-section of the bubble is simulated. At the outer boundaries, a pressure-outlet condition is applied. The fluid is the refrigerant HFE-7100, the initial velocity field of the domain is = 0 u , and the initial pressure and temperature of the surrounding superheated liquid are 50,000 Pa and 317.95 K respectively, corresponding to a liquid superheat level of 5 K. Inside the bubble, the initial conditions are the saturation conditions, where the bubble internal pressure is initialized according to the Young-Laplace equation and the temperature of the bubble is determined as the saturation temperature at that pressure. The surface tension coefficient is set to a constant value of 0.0136 N/m resulting in an initial internal bubble pressure of 50,272 Pa and an initial temperature of 313.1 K. In the liquid close to the interface, a temperature profile from the analytic solution is implemented in order to avoid the sharp temperature jump at the interface and to increase the accuracy of the initial evaporation rate modelling. The solution is only valid in the liquid phase and the vapour phase temperature is assumed uniform. For the HFE-7100 fluid, a value for the evaporation coefficient χ has not been found, but, for water, Marek et al. (2001) have given numerous values for different situations and shown that it may vary in the range of roughly is chosen that gives reasonable agreement to the analytical solution, and a mesh independence study is performed to investigate whether the results converge after successive mesh refinement. The refined mesh resolutions are chosen such that the cell size is reduced by a factor 2 for each consecutive refinement, resulting in mesh resolutions of 25, 50, 100, and 200 cells per initial bubble diameter. Figure 3 shows the bubble radius over time for the simulations performed using different refinement levels, together with the analytical solution by Scriven (1959). The difference between the analytical and simulation results during the initial growth period is to be expected since the analytical solution starts at a zero bubble radius and zero velocity while the simulations start at a non-zero radius and zero velocity. After the initial period, where the growth rate is inertia controlled, the bubble growth rate becomes Fig. 3 Bubble radius during the simulation of validation case 1 and Scriven (1959) analytical solution. The fluid is refrigerant HFE-7100. In the simulations, a bubble is initialized with a 0.1 mm radius and is surrounded by a liquid that is superheated 5 K to a temperature of 317.95 K. The operating pressure is set to 50,000 Pa. For successive mesh refinement, the solutions converge and the asymptotic behaviour tends to the analytical solution. The mesh resolutions are specified as the number of computational cells per initial bubble diameter. Due to different initial conditions the early, inertia controlled, growth phase is not comparable between the simulations and analytical solution. But, for the second, thermally controlled, growth phase, the bubble growth rate should be similar. thermally controlled and the results can be compared. The critical time that determines this transition can be estimated as (Brennen, 1995): where v ( ) p T ¥ is the saturation pressure corresponding to the temperature of the surrounding liquid, p ¥ is the pressure of the surrounding liquid, and l ρ is the density of the liquid. Σ is a parameter related to the thermo-physical properties of the multiphase system: where v ρ is the saturation density of the vapour at the temperature of the surrounding liquid, l p c is the specific heat capacity of the liquid, T ¥ is the temperature of the surrounding liquid, and l α is the thermal diffusivity of the liquid. The resulting critical time for the present validation case is around 1.4 μs. In the present simulation with simulation time of the order of milliseconds, the major part of the bubble evolution is thermally controlled. The slope of the curve as time increases is in fair agreement with the analytical solution. An even better agreement may be achieved by further fine tuning of the evaporation coefficient. Considering mesh independence, it is indicated in Fig. 3 that, as the cell size is reduced, the corresponding simulation results converge. This is a reasonable behaviour, since during a thermally controlled bubble growth, it is the ability of the liquid to transport heat to the interface that is the limiting factor for bubble growth and as the temperature gradient in the liquid becomes more accurately resolved, the heat flux stabilizes to a certain value.

Validation case 2: Bubble growth and collapse
The second validation case is performed to qualitatively assess the ability of the phase change model to handle both evaporation and condensation during the growth and collapse of a vapour bubble where only a region of the liquid around the bubble is superheated. No suitable analytical solution is found for validation and therefore, a case similar to the evaporation validation case, shown in Fig. 2, is constructed using the same HFE-7100 fluid. The bubble is initialized with the same radius of 0.1 mm, temperature of 313.1 K, and pressure of 50,272 Pa. But, instead of being surrounded by an infinitely extended superheated liquid, the liquid is only superheated in the region between the bubble and R = 0.15 mm with a temperature of 317.95 K. Outside this region there is liquid at a temperature of 310 K, i.e., below the saturation temperature of 312.95 K that corresponds to the initial liquid pressure of 50,000 Pa.
The case is constructed so that the bubble will expand until all superheated liquid either evaporates or gets cooled below the saturation temperature. After that, the bubble will be hotter than the surrounding cold liquid and should consequently start to condense and shrink until it disappears. The laser-induced bubble evolution is governed by the same principles, although at smaller scales and much faster dynamics. If the model is able to handle this validation case, it should therefore, in principle, be able to model the governing physical phenomena of the laser-induced bubble as well.
The evolution of the bubble properties during the simulation is shown in Fig. 4 using a semi-logarithmic scale to better visualize the early bubble growth period. Due to the fast evaporation rate of the superheated liquid, the bubble rapidly increases in size, reaching a maximum radius of about 0.26 mm before starting to contract due to condensation of vapour. Condensation begins once the vapour at the interface gets colder than the saturation temperature, which, in turn, happens when all the superheated liquid either evaporates or is sufficiently cooled down by the surrounding liquid. Figure 4(b) shows the bubble mean pressure during the simulation. During the first few milliseconds, there is a spike in the bubble mean pressure which is an effect of liquid evaporation and the surrounding liquid inertia. Along with the bubble interior pressure from the simulation, a calculated Laplace pressure is plotted. The Laplace pressure is given by the expression: where p D is the pressure difference between the bubble interior and exterior pressures at equilibrium conditions, σ is the surface tension coefficient, and R is the bubble radius. The Laplace pressure can be used as an indication of whether the current interior bubble pressure should result in bubble growth or contraction. If the pressure in the bubble is greater than the Laplace pressure, the bubble will push outwards on the liquid, resulting in growth, and the opposite will take place if it is lower. In Fig. 4(b) the curves of the bubble pressure and the Laplace pressure cross each other at about the same time as the bubble radius evolution in Fig. 4(a) shifts from growth to contraction. There is an initial spike in the mean bubble temperature in Fig. 4(c) as well. This spike can be explained by the combination of an increase in the bubble pressure and that the surrounding superheated liquid, at 317.95 K, heats the vapour. It is the temperature at the interface that determines the phase change rate, so even though the mean temperature of the bubble is above the saturation temperature well after bubble contraction has begun, it takes time for the entire bubble to cool down.
During condensation of vapour, latent heat is released at the vapour side of the interface. This increases the temperature and thereby lowers the condensation rate. Therefore, in contrast to evaporation, the condensation rate is governed by the vapour's ability to transport heat away from the interface. This is mainly done through conduction which is directly proportional to the temperature gradient and the thermal conductivity of the fluid. The temperature difference between the phases for this case is higher during the evaporation period, around 5 K, than during the condensation period, around 3 K, and the thermal conductivity is far higher in the liquid than in the vapour phase. This suggests that the evaporation rate of the liquid during the bubble growth phase should be much higher than the condensation rate of the vapour during the bubble contraction phase. In Fig. 4(a) these phenomena are visible in the form of much faster bubble growth rate ( 130 mm/s) than contraction rate ( -4 mm/s). Previous experimental investigations of the formation and collapse of a laser-induced vapour bubble in a microtube filled with water, and with added red dye, report similar vapour bubble dynamics to those in our validation case, namely, rapid expansion and slower contraction (Quinto-Su et al., 2009;Sun et al., 2009). Although our validation case does not have the same characteristics as the mentioned experiments, the qualitative behaviour of the implemented phase change model indicates that the Fig. 4 Evolution of the bubble radius (a), pressure (b), and temperature (c) during the validation case 2. The fluid is refrigerant HFE-7100. A bubble is initialized with a radius of 0.1 mm and surrounded by a liquid to a radius of 0.15 mm that is superheated 5 K to a temperature of 317.95 K. The rest of the domain contains of liquid below the boiling conditions, at 310 K. The operating pressure is set to 50,000 Pa. model is able to capture the governing physical phenomena for expansion and collapse of vapour bubbles formed by superheating of the liquid.

Simulation of a laser-induced vapour bubble
With a phase change model capable of incorporating both evaporation and condensation effects, the relevant physical phenomena governing the dynamics of a laser-induced vapour bubble are accounted for in our numerical framework. In this section we describe such a simulation performed to evaluate the feasibility of the numerical framework to simulate the early growth period of the bubble. Although uncertainties regarding the accuracy of the phase change model exist for this type of rapid phase change problem, a qualitative result should be attainable for testing the suggested hypothesis, about the crystal nucleation around the bubble.
Several simplifying assumptions are made during the setup of the simulation and pointed out in the description below, together with suggestions for future investigations and improvements.

Problem setup
The problem setup has been chosen from the previously mentioned experiments (Soare et al., 2011;Soare, 2014). A schematic view of the geometrical setup is shown in Fig. 5, where an aqueous solution of an inorganic salt, (NH4)2SO4, is placed between two glass slides, 50 μm apart. Then, a 6 ns laser pulse of 0.27 mJ was focused to the center of the solution resulting in a beam width of about 20 μm. The focused laser pulse superheated the liquid, resulting in a vapour bubble that grew explosively. The bubble rapidly exceeded the distance between the two glass plates, thereby growing in a seemingly two-dimensional manner in the plane and reached a maximum radius of roughly 700 μm with an interface velocity of well above 20 m/s. The whole event, from bubble nucleation to collapse, took about 200 μs. The setup chosen for the simulations is a somewhat simplified representation of these experiments. Pure water at 20 °C is used instead of a salt solution, so that only a single species in each phase needs to be considered. The geometrical setup and laser pulse properties are kept as in the experiment. As previously discussed, it is the early bubble growth period that is of interest in this simulation, since crystal nucleation most likely occurs just after bubble nucleation when the evaporation rate is at its maximum.

Simulation setup
The equation of state for the vapour phase is the ideal gas law even though this assumption may not be accurate for vapour under high pressure and close to saturation conditions. In future simulations, a more complex but also a more accurate alternative would be to adopt a real gas model.
Due to the exceptionally fast dynamics of the problem, large pressure spikes in the liquid close to the bubble are observed in the simulations. To better resolve the effects of those pressure variations, a compressible liquid equation of state, the Tait equation, is employed. This additional degree of freedom, as compared to the common incompressibility assumption, increases the stability of the simulations. It is not trivial to determine a correct value for the evaporation coefficient, χ, for this type of problem. As previously stated, Marek et al. (2001) have shown that it may vary in the range of 3 10 1 χ -< £ for water in different situations. Therefore, validation against reliable experiments has to be made in order to get a good estimate. In the simulations presented in this section a value of 0 01 χ = . is chosen that seems to produce results that are in fair agreement to experimentally observed growth rates. Surely, a further fine-tuning of this parameter might produce results with even better agreement, but, due to the lack of detailed validation data in the bubble growth period of interest, such a study is not performed in this work.
The choice of advection scheme for the color function is important since obtaining a sharp interface will increase the accuracy of the phase change model. Therefore, we have chosen the PLIC interface reconstruction scheme that reduces the smearing of the interface.
Surface tension is introduced using the Continuum Surface Force model developed by Brackbill et al. (1992) with a constant value for the surface tension coefficient between liquid water and water vapour of σ = 0.0589 N/m. In this simulation case, the temperature of the interface varies significantly and should affect the value of the surface tension coefficient, but this effect is not considered in this study.

Computational domain
To reduce computational cost the same assumptions of a 2D axisymmetric and uniform grid, as in the validation cases, are employed. A schematic representation of the resulting computational domain is shown in Fig. 6. The distance between the center of the bubble and the glass wall is set to 25 μm but the required distance between the bubble center and the outlet is, however, not that obvious. In the experiments the bubble reached a maximum radius of about 700 μm, but since the focus of this simulation is on the early bubble growth phase, a domain of 25 100 μm is chosen with a grid resolution of 100 400 cells.

Boundary conditions
The operating pressure of the system is specified as 101,300 Pa. Fig. 6 Computational domain used in the simulations of the laser-induced vapour bubble. A 2D axisymmetric, uniform mesh is adopted where the bubble is initialized with its center at the lower left corner. The radial axis has a symmetry boundary condition and the axial axis an axisymmetric boundary condition. The top boundary has a no-slip wall boundary condition, representing the glass wall, and the right boundary a pressure-outlet boundary condition, representing liquid far away from the bubble.
The outlet boundary is specified as a pressure-outlet boundary condition with a constant pressure equal to the operating pressure. In the case of backflow through the boundary, backflow conditions are set to a pure liquid phase with a temperature of 293.15 K. The glass wall boundary have a no-slip condition and the temperature is assumed constant due to the short duration time of the simulation. The part of the upper boundary closest to the z-axis is where the laser beam irradiates through the glass. Due to the non-zero fractional absorption of light, this part of the glass will also attain a higher temperature after laser irradiation. According to the study with the experiments this portion of the glass reached about 363 K whereas the rest of the boundary remained at 293.15 K (Soare, 2014). The radius of the focused laser beam was about 10 μm. Thus, in the simulation, the part of the glass wall from the rotational axis and 10 μm to the right is set to a constant temperature of 363 K, and the rest to 293.15 K.

Initial conditions
The initialization of a nucleating bubble induced by a laser is not trivial and may be done in various ways. To study the whole bubble growth phase it would be necessary to initialize and resolve a bubble at the nucleation scale, most probably from one or several existing nano bubbles of non-condensable gas (Soare, 2014). Instead, due to computational limitations, a scale of the order of micro meters is chosen. During the first 6 ns of the simulation a volumetric energy source is added to the left 10 μm of liquid. After this irradiation period, it is assumed that the bubble is formed with an initial radius of 3 μm, and it is introduced with its center at the lower left corner of the domain. The initial pressure and temperature of the bubble are set to the corresponding gauge Laplace pressure of 39,267 Pa, and the temperature to the corresponding saturation temperature, 382.7 K.
The laser pulse energy should supposedly be absorbed in a cone shaped portion of the liquid due to the focusing lens concentrating the laser beam. Also, the intensity of the irradiation should decrease exponentially as the light passes through the absorbing liquid, making the liquid closer to the laser source hotter. However, these conditions and geometry of the laser beam are not known from the experiments. Therefore the laser-induced heating in the simulations is assumed uniform and introduced as a cylinder. The study in Soare (2014) indicates that the liquid reached about 490 K, which is chosen as the temperature achieved after irradiation in the simulations. It would, however, be useful to obtain experimental values of the temperature distribution and the bubble evolution properties during and after laser irradiation, since those values influence the early bubble growth phase. An analytical temperature profile in the liquid close to the interface is calculated using the solution by Scriven (1959). This profile shows that the entire temperature drop from the surrounding superheated liquid to saturation temperature at the interface occurs within the distance of a single computational cell. The validity of the analytical solution in this extreme case is not certain and implementing the temperature profile in the simulation would require an excessively fine computational grid. Therefore, no initial analytic temperature profile is introduced in this case.

Time step
The governing equations are discretized in time using a first order implicit scheme. A variable time step size is adopted to reduce the computational cost when dealing with large differences in fluid velocities and phase change rates during the simulations. When introducing potentially large mass source terms it is noticed that nonphysical results are produced if the variable time step size is limited only by a Courant number of Co = 0.25. The errors occur in connection to the large mass source terms and therefore additional time step limitations are introduced to reduce the amount of mass that is added or removed during a time step. A limitation for the time step size is that the total flux going out of a cell should not be able to completely empty a cell during a time step. This criterion can be expressed as where f denotes face values and Ω the whole computational domain. Also, a negative mass source term large enough to remove all fluid in a cell during a single time step is probable to cause errors. Therefore, a ratio that further limits the time step size is q q q 1 0 where q represents the fluid phase. This constraint should hold for all cells in the domain and let us define another time step criterion as where a maximum allowable ratio 1 β < is introduced. Different values of β are tested, and, for the laser-induced bubble case, a value of 0 0001 β = .
produces stable simulations.
Too large changes, > 10%, in the time step size are also preferably avoided when using an implicit scheme. Therefore, the smallest of the time step limitations, denoted limit t D , is computed for each time step and the next time step size is chosen as This approach ensures that the time step limitations are fulfilled and, at the same time, make us use an as large as possible time step. The volume fraction field is advected using an explicit formulation with a refined time-step that is always smaller or equal to the time-step used for the governing equations, but also restricted by the capillary time-step constraint discussed in Denner et al. (2015).

Results and discussion
In this section, we present results from a numerical simulation of a laser-induced vapour bubble using a similar setup as in the previously mentioned experiments. It is the evolution of the conditions in the thin layer of liquid around the bubble that is of special interest during this study. It is hypothesized that the crystals are nucleated at the maximum rate of evaporation, which, in this simulation, takes place simultaneously with the maximum cooling rate of the interfacial liquid. When the interfacial liquid is cooled below saturation conditions, the evaporation stops and primary crystal nucleation within that liquid becomes less probable. Therefore it is assumed that the crystal nucleation period is captured in the simulation if the liquid has been cooled significantly below saturation conditions, and once this has occurred, the simulation is stopped.
Qualitative comparisons with the available experimental results can be made in order to evaluate the ability of the numerical framework to capture the governing physics of the problem. Also, the plausibility of the suggested hypothesis about the mechanisms behind crystal nucleation can be discussed by examining the conditions of the liquid around the bubble where the crystals are supposedly nucleated. The simulation results are, however, not intended for definite quantitative conclusions about the conditions in the liquid around a real laser-induced vapour bubble. For such conclusions to be made, more detailed experimental data and validation of the simulation method are necessary.
The temperature of the liquid at the bubble interface and the corresponding saturation temperature are shown in Fig. 7. The interfacial liquid temperature is defined as the temperature at the location where the volume fraction of liquid is equal to 0.99 along the radial axis, and at z = 0 in the computational domain. Along the radial axis there is a minimum distance between the bubble interface and the outer liquid, not heated by the laser pulse. Consequently, it is the place where the cooling of the interfacial liquid should be most rapid, since it is cooled both by evaporation, from one side, at the gas-liquid boundary, and through heat loss to the outer liquid, from the other side. The saturation temperature is computed based on the interfacial liquid pressure and the Clausius-Clapeyron relation. Initially, the pressure in the interfacial liquid is at atmospheric conditions, but as the bubble starts to expand, the pressure rises to around 1.1 MPa and the corresponding saturation temperature increases from 373 to about 460 K. The liquid temperature is at first reduced only a few degrees due to evaporation but then drops rapidly below the saturation temperature at around 0.5 μs, where the maximum cooling Fig. 7 Temperature of the interfacial liquid, taken along the radial axis at z = 0 in the computational domain, during the simulation of a laser-induced vapour bubble. At the same location, the saturation temperature is computed that corresponds to the interface liquid pressure. The superheated liquid that initially surrounds the bubble has an initial temperature of 494 K and the rest of the domain is liquid water at 293.15 K. Once the interfacial liquid temperature is below saturation conditions, the evaporation stops and crystal nucleation becomes less probable. It is therefore assumed that an eventual crystal nucleation period is captured before the end of the simulation. rate is obtained at the rate of about -2 × 10 8 K/s. It is in this phase of the simulation that the super-heated liquid, along the radial axis, is reduced to a very thin layer between the bubble interface and the outer, cool, liquid. Heat losses to the outer liquid and cooling by evaporation both contribute in this area to the intense cooling rate of the interfacial liquid. Based on the proposed hypothesis, that it is the increase in the solute concentration and the simultaneous cooling of the liquid that are the mechanisms behind crystal nucleation, it would therefore seem probable for crystals to form at this location. However, without information about supersaturation levels of solute in the interfacial liquid, it is not possible to determine if these conditions are sufficient for primary nucleation of crystals. Figure 8 shows the evolution of the bubble radius from the simulation and the available experimental data from Soare (2014). Unfortunately, the temporal resolution of the experimental data is 4 μs, which is more than twice the period of interest in this simulation. Therefore, the experimental data cannot be used for a direct comparison with the simulation results, but rather as an indication of whether the trend of the bubble growth rate in the simulation is reasonable. The bubble in the simulation is initialized with Fig. 8 Evolution of the bubble radius during the simulation of a laser-induced vapour bubble. The available experimental data from Soare (2014) have a temporal resolution of 4 μs and can therefore not be used for a direct comparison with the simulation results. Still, the experimental data can give an indication of whether the bubble growth rate obtained by the simulation is reasonable. Note that the first experimental data point in the figure is just an indication and should be at t = 0, R  10 -9 . In the simulation, the bubble has an initial radius of 3 μm and expands rapidly in size, reaching a radius of 44 μm at the end of the simulation, 1.9 μs after the laser pulse. In the experiments, the bubble starts to grow at a nano scale and has a non-zero interface velocity when reaching the initial bubble radius of 3 μm used in the simulation. The bubble growth rate, in the simulation, at these early times should therefore deviate from experiments. The interface velocity, after the initial acceleration period, is around 23 m/s and is in the same range as those observed in experiments. a radius of 3 μm and expands rapidly in size, reaching a radius of around 44 μm at the end of the simulation, 1.9 μs after the laser pulse. After the early acceleration period, the interface velocity is around 23 m/s. In the experiments, the bubble starts to grow at a nano scale and has a non-zero interface velocity when reaching the initial bubble radius of 3 μm used in the simulation. The bubble growth rate in the simulation at these early times should therefore deviate from that in the experiments. As discussed before, the critical time, according to Brennen (1995), when the inertia and thermal effects are in the same order of magnitude is around 3 ns for an identical bubble that is suspended in an infinitely extended and uniformly superheated liquid. In the simulation, the bubble is not suspended in an infinite liquid, but since the total simulation period is several orders of magnitude larger than the critical time, the major part of the bubble growth period should be thermally controlled. The interface velocity, after the initial acceleration period, is in the same range as that observed in the experiments, which indicates that the fundamental features of the bubble expansion phase are captured. We acknowledge that detailed experimental data would be needed for gaining deeper insight into the dynamics of an early bubble growth under such conditions.
During the simulation, the mass conservation in the entire domain is monitored to ensure that the phase change model does not produce nonphysical results. In Fig. 9, the mass in the system and a theoretical mass balance, both normalized by the mass in the system at t = 0, are shown. As the vapour bubble expands, the heavier liquid is flowing out from the domain and the total mass is reduced. The Fig. 9 Mass conservation within the computational domain during the simulation of a laser-induced vapour bubble. The figure shows the actual mass in the system at all time steps and a mass balance that represents the change of mass due to flux across the boundaries. All values are normalized with the actual mass at t = 0. The difference between the two lines corresponds to a nonphysical change of mass inside the computational domain. Throughout the simulation this difference is never greater than 0.2%. actual mass is the total mass in the system at all time steps and the mass change due to flux across the boundaries is shown as a mass balance. The difference between the two quantities corresponds to a nonphysical change of mass inside the computational domain. Since this difference is less than 0.2% at all times, this effect is assumed to not significantly affect the results.

Conclusions
We formulate in this work a numerical method that resolves the conditions of the liquid around a laser-induced vapour bubble. Phase change of the fluid is taken into account using a model for interfacial mass transfer that is incorporated into the VOF framework. Two validation cases are presented that aim at evaluating the framework's ability to resolve the dynamics of vapour bubbles in superheated liquids. Then, a numerical analysis is performed of a laser-induced vapour bubble using a somewhat simplified setup from the experiments (Soare et al., 2011;Soare, 2014). The results show that the dynamics of the bubble are in reasonable agreement with the experimental results. In the same time, we realize that more detailed experimental data would be needed for validation of the method. The mechanisms behind the experimentally observed crystal nucleation around a laserinduced vapour bubble are not entirely clear but a hypothesis is that it is the increase in the solute concentration and simultaneous cooling of the liquid around the bubble, due to evaporation, that are needed for the nucleation process to take place. Our numerical simulations show that the interfacial liquid temperature is rapidly reduced due to the combined effect of evaporation and, especially near the radial axis of the computational domain, due to heat loss to the surrounding liquid. The rapid cooling occurs in the same time frame as the maximum rate of evaporation is obtained, and it would therefore seem probable for crystals to form during this time. Further validation of the framework and incorporation of solute concentrations in the numerical simulations are needed so that the latter indeed become a tool to fully understand the mechanisms that lead to crystal nucleation.