New semi-analytical solution of the problem of vapor bubble growth in superheated liquid

This paper presents a mathematical model of the vapor bubble growth in an initially uniformly superheated liquid. This model takes into account simultaneously the dynamic and thermal effects and includes the well-known classical equations: the Rayleigh equation and the heat conductivity equation, written with consideration of specifics associated with the process of liquid evaporation. We have obtained a semi-analytical solution to the problem, which consists in reducing the initial boundary value problem with a moving boundary to a system of ordinary differential equations of the first order, valid in a wide range of operating parameters of the process at all its stages: from inertial to thermal, including the transitional one. It is shown that at large times this solution is consistent with the known solutions of other authors obtained in the framework of the energy thermal model, in particular, for the high Jacob numbers, it is consistent with the Plesset–Zwick solution.

The growth of a vapor bubble is one of the main processes determining boiling, and it is a complex problem described by the equations of hydrodynamics and heat transfer. Obviously, the detailed study of the mechanism of this process is of great importance for both theory and practice. The history of numerous attempts to obtain the law of vapor bubble growth dates back more than a century, beginning with classical works [1][2][3][4][5][6] and ending with modern research [7][8][9][10][11] . They were based either on highly simplified mathematical models 12,13 , or on a fairly complete system of classical equations as applied to the problem under consideration, solved numerically [14][15][16][17][18] . Similar studies can be found for gas bubbles [19][20][21][22][23][24] with the only difference that the mass transfer processes are the governing ones instead of the heat transfer processes.
In the most general statement, the growth of a vapor bubble is determined by various factors: dynamic, kinetic, thermal, etc. Depending on the properties of the two-phase system under consideration, at some stage of the process, the factors that determine it change. Therefore, to solve the problem analytically, various assumptions are used, which allow the construction of simplified mathematical models with one main mechanism of the process control. Let us list the the most common models 25 .
Dynamic inertial model In this model, the vapor pressure in the bubble during the entire process is kept constant and equal to the saturation pressure at the initial liquid temperature, and the pressure drop determines the bubble growth. The temperature of liquid at the bubble boundary does not change and equals the temperature at a great distance from the bubble, i.e., the effect of liquid cooling near the interface due to evaporation is neglected. The analytical solution under these assumptions is obtained by simple integration of the Rayleigh equation. In this model, the dependence of the bubble radius on time is linear.
Dynamic viscous model In this model, the bubble growth rate is limited by the viscous forces. Obviously, such a model can be applied only to highly viscous liquid media. The initial equation for this limiting case is obtained from the Rayleigh equation (neglecting the inertial terms), which is easily integrated. The dependence of the bubble radius on time in this model is exponential.
Energy molecular-kinetic model In this model, an attempt to take into account the kinetics of phase transformation is made. The bubble growth rate in this case is determined, among other things, by the rate of liquid evaporation from the superheated interface. At that, the vapor in the bubble during its growth is in the unsaturated state.
Energy thermal model This model is most common in the literature. In this model, the bubble growth is determined by the supply of heat to the interface from external superheated layers of liquid. At that, all supplied heat is spent on evaporation. This requires knowledge of temperature distribution in the liquid volume. The vapor pressure in the bubble throughout the process is kept constant and equal to the pressure of the surrounding liquid. www.nature.com/scientificreports/ The vapor in the bubble is in the saturated state. The dependence of the bubble radius on time in this model has the root nature. Moreover, the proportionality coefficient (which is sometimes called the growth modulus) is a function of the Jacob number and the vapor -liquid density ratio and it is found by different authors in different ways. The most famous is the Plesset-Zwick formula 1 , which is valid for high Jacob numbers. A sufficiently complete analytical solution to the problem in integral form was obtained by Scriven 4 . Along with this, various empirical approximation dependencies are widely used 26,27 . It is obvious that the energy thermal model, despite its widespread use, as well as all other mathematical models, has a limited scope. In particular, a number of approximations made in this model lead to a solution with an infinite rate of the bubble growth at the initial time, and this is physically incorrect. In this regard, the attempts to create various hybrid models that take into account simultaneously the dynamic and thermal effects are being made. For this, the Rayleigh equation should be solved together with the heat conductivity equation, which, as practice shows, is associated with significant difficulties. In this paper, we tried to find a relatively simple semi-analytical solution that would simultaneously describe both inertial and thermal effects and which would become an alternative to direct numerical calculations. problem statement Main equations. Let us consider the growth of a single supercritical vapor bubble, formed at the initial time moment in uniformly superheated liquid. The system of equations describing the behavior of such a bubble includes well-known classical equations written in relation to the problem under consideration, taking into account the specifics associated with the process of liquid evaporation. When stating the problem, we will use the following assumptions. The temperature and vapor pressure in the bubble are uniform. The vapor in the bubble is stationary and it stays saturated throughout the process. The liquid around the bubble is viscous and incompressible. The flow is spherically symmetric.
The velocity field around the bubble is obtained from the continuity equation and has the following form: where v l is the radial velocity of liquid, r is the radial coordinate with the beginning in the center of the bubble and R is bubble radius. Hereinafter, the subscripts "l" and "v" refer to the liquid and vapor phases, respectively and subscript "R" refers to the value at the interface. The equation of liquid motion can be written as where t is the time, p is the pressure and ρ is the density. The last equation can be integrated along the radial coordinate, taking into account (1): Hereinafter, superscripts "i" and "f" correspond to the initial and final states, respectively. Note, that the liquid velocity at the interface v lR differs from the rate of bubble growth Ṙ due to the phase transition: v lR =Ṙ − j/ρ l , where j is the density of the mass flux at the interface and ˙ ≡ d( )/dt is the time derivative of arbitrary function . In the case of no phase transition, Eq. (2) becomes well-known Rayleigh equation. The boundary conditions representing the mass, momentum and energy conservation laws can be written as follows 28 : where σ is the surface tension, η is the dynamic viscosity, is the thermal conductivity, T is the temperature and L is the specific heat of phase transition.
The presented system of equations is closed by a thermal boundary-value problem solved inside the liquid surrounding bubble. The temperature field dynamics is described by the heat conductivity equation (volumetric heat release caused by viscous dissipation is neglected): where c is the heat capacity. Thermophysical properties of liquid will be considered constant.
At the initial moment of time, the temperature of liquid is assumed to be uniform and higher than the vapor saturation temperature at pressure p li (the liquid is superheated): (T l ) t=0 = T li > T s (p li ) . At a great distance from the bubble, the temperature field remains unperturbed during its growth: (T l ) r→∞ = T li . On the interface, www.nature.com/scientificreports/ the condition of local thermodynamic equilibrium is satisfied: (T l ) r=R = T s (p v ) (the superscript "s" refers to the value of the function on the saturation line). It should be noted that the temperature of liquid at the interface changes over time.
The problem should be supplemented by the equation of vapor state in the bubble and temperature dependence of the saturation pressure in the following general form: Model equations and more accurate empirical dependencies can be used, which are numerous in the literature. We should note that dependencies (7), (8) uniquely determine functions T s (ρ v ) and p s (ρ v ) , which will be used below.
According to the literature on the subject of vapor bubble growth the vapor can be considered an ideal gas and the Mendeleev-Clapeyron equation can be used as the equation of state. Moreover, to find the temperature dependence of the saturated vapor pressure, the Clapeyron-Clausius equation can be applied. It must be said that this somewhat narrows the generality of both the statement of the problem itself and the solutions obtained. However, in most cases this is enough to ensure the acceptable accuracy of the desired solution. An example of using these equations in a problem is given below.
In the framework of assumptions made, the formulated system of equations is closed and fully describes the dynamics of a vapor bubble in uniformly superheated liquid.
nondimensionalization. Further, we will use the following dimensionless variables: = is the initial superheating of liquid and �p i = p s (T li ) − p li is the initial overpressure caused by superheating; χ = r/R is the dimensionless radial coordinate; R = R/R 0 ; v = v/v 0 ; τ = t/t 0 are dimensionless radius, speed and time, respectively; ρ v = ρ v /ρ l is the dimensionless vapor density in the bubble. We define the characteristic size R 0 , the speed v 0 , and the process time t 0 in such a way that the characteristic time of the dynamic stage and the characteristic time of reaching the thermal stage of the bubble growth are In this case, R 0 = a l ρ l /�p i ; v 0 = �p i /ρ l ; t 0 = a l ρ l /�p i . The dimensionless variables introduced in this way allow us to abstract as much as possible from the properties of liquid under consideration and to focus on the behavior of the studied dynamic system. As for the choice of characteristic values, they are most convenient for the analysis of the transition stage of the process. All similarity criteria that arise in dimensionless equations and determine the process under consideration are introduced below.

problem solving and analysis of results
Semi-analytical solution. Let us turn from the variables t and r to the variables τ and χ , thereby reducing the formulated heat problem to the problem with fixed boundaries (similar to how it was done in a number of publications 21,22 , when solving the problem of diffusion growth of a gas bubble).
The heat conductivity equation will take the following form: where α =v lRR and β =ŔR represent the functions of time; ´ ≡ d( )/dτ is the function derivative by dimensionless time τ. Initial and boundary conditions: In the general case, the boundary-value problem (9), (10) (even without consideration of inertial effects) can be solved only numerically. However, there is a way that allows us to find a very good approximate solution. To do this, we find a steady-state solution of Eq. (9), taking into account the boundary conditions (10): where I(χ, α, β) = Numerical calculations show that the time of establishing a quasi-stationary state (in variables τ and χ ) is extremely short on the scale of characteristic time of the whole process. This means that solution (11) describes with good accuracy the dynamics of a temperature field forming around a bubble throughout the entire process.
Making the initial system of Eqs. (2)-(5) dimensionless and using the found temperature profile (11) in the dimensionless Eq. (5), we obtain a system of differential equations of the first order: If vapor is considered an ideal gas, then the Mendeleev-Clapeyron equation should be used as the equation of state: , where R g = R /M is the reduced gas constant; R is the universal gas constant; M is the molar mass of vapor. In the dimensionless form it is written as follows: The Clapeyron-Clausius equation takes form: In the dimensionless from it is written as: where ε = L /(R g T li ).
Thus, the problem is reduced to solving the system of three ordinary differential equations of the form ý = f (y) , where y = (α,ρ v ,R) is the desired vector function, and function f (y) is set by the right side of Eq. (12).
We should note that integral function I * (α, β) has obvious asymptotic approximations, and in some cases they can be used to simplify the analysis of the resulting solution: • at α, β ≪ 1 : I * (α, β) ≈ 1; • at α, β 1 : where erfc(z) = ∞ z e −ζ 2 dζ. Fig. 1 illustrates the the dynamics of the processes involved in bubble growth in water under the atmospheric pressure and initial superheating 100 K. It can be seen from the Fig. 1a that the temperature field with time tends asymptotically to the stationary distribution (in variables τ , χ ), when the vapor pressure in the bubble is almost equal to the pressure of surrounding liquid.
Time dependencies of the bubble radius and vapor density in the bubble, plotted on the basis of the found semi-analytical solution (12), (13), are presented in Fig. 1b. These dependencies clearly illustrate the features of the growth mechanism of a vapor bubble in superheated liquid, combining the dynamic and thermal effects. We should note that the difference in the calculation results obtained by direct numerical simulation and on the basis of Eqs. (12), (13) turned out to be extremely small, which indicates the consistency of the found semianalytical solution. Dependence of the bubble radius on time, constructed within the framework of the thermal energy (see the next paragraph) is also shown in Fig. 1b. The form of the plotted curves emphasizes once again that this model is essentially asymptotic and does not describe the considered process in the entire time range (as it was mentioned in Introduction). Figure 2 shows the applicability of the found semi-analytical solution for a wide range of initial superheating. The chosen values of initial superheating T i = 10 K, T i = 50 K and T i = 100 K approximately correspond to the values of Jacob number Ja of 30, 160 and 320 respectively. It can be seen that the results of direct numerical calculation of the boundary problem (2)-(6) differ insignificantly from the results obtained using found semianalytical solution (12), (13) for all values of initial superheating under consideration. This means that obtained semi-anlytical solution can be used as an alternative to direct numerical calculations in solving the problems involving vapor bubble growth. Thermal stage: self-similar solution. According to the simple analysis of solution (12), (13) obtained above, as the bubble grows, the pressure in it decreases gradually and tends asymptotically to the pressure of (12) Scientific RepoRtS | (2020) 10:16526 | https://doi.org/10.1038/s41598-020-73596-x www.nature.com/scientificreports/ the surrounding liquid ( v → 0 ), and the vapor density in the bubble and the liquid temperature at the bubble boundary tend to the certain constant values: ρ v →ρ vf ; (� l ) χ =1 → −1 , where ρ vf is the vapor density at pressure p li (calculated on the saturation line). It should be noted that if we do not take into account any exotic cases, then ρ vf ≪ 1 . At this stage of the process, called the thermal stage, which is essentially asymptotic and can be described in terms of the energy thermal model, the bubble growth is determined exclusively by heat supply to the interface, the temperature field around the bubble becomes stationary, functions α(τ ) and β(τ ) become constant, and the solution to the boundary value problem becomes self-similar (the self-similar variable here is variable χ , and in this case it is proportional to dimensionless complex r/ √ a l t): where κ = 1 −ρ vf . The similar solution was obtained by Scriven 4 .
The coefficient β f occurring in Eq. (16) at this stage is a function of only Jacob number and is found from the implicit equation (which follows directly from Eq. (13)):   (16), (17), the temperature field does not depend explicitly on time; therefore, an assumption about the quasi-stationary nature of the process (in variables τ , χ ) made above is satisfied here automatically.
In the case of low and high superheating, in the approximation of ρ vf ≈ 0 , which is often used in the literature, there are obvious asymptotic approximations. For Ja ≪ 1 : β f ≈ Ja . For Ja ≫ 1 : β f ≈ (6/π)(Ja + 4/9) 2 . The latter coincides with the well-known Plesset-Zwick solution 1 (if we do not take into account correction 4/9 to the Jacob number, which significantly improves this asymptotic approximation). Figure 3 illustrates the dependence of coefficient β f on Jacob number Ja calculated by formula (17) and using asymptotic approximations. Let us note that there is some discrepancy between solution (17) and asymptotic approximations at high Jacob numbers, since the latter was obtained with the assumption of ρ vf ≈ 0 , which is not entirely allowable for β f 1/ρ vf . From the above analysis, we can conclude that due to the good correspondence of direct numerical and approximation solutions, the thermal stage of the vapor bubble growth has been studied in sufficient detail in the literature.

conclusions
A mathematical model of the vapor bubble growth in initially uniformly superheated liquid has been formulated. This model takes into account simultaneously the dynamic and thermal effects and includes the well-known classical equations: the Rayleigh equation and the energy equation, written in relation to the studied problem with consideration of the specifics associated with the process of liquid evaporation. The following assumptions were used: the temperature and pressure of vapor in the bubble are homogeneous; vapor in the bubble is stationary and it is in a saturated state throughout the process.
It is shown that the presented problem can be reduced to solving a system of three ordinary differential equations of the first order: a new semi-analytical solution of the problem, which simplifies the analysis of the process greatly. It is shown that the solution obtained is in good agreement with direct numerical calculations in a wide range of superheating and at all stages of the process, including the transition one, which must be taken into account, especially if bubble growth in a strongly superheated liquid is considered. The dynamics of changes in the temperature field in the volume of liquid and gas density in the bubble is illustrated. The time dependence of the bubble radius is found.
It is shown that at large times the process of bubble growth is determined exclusively by the supply of heat to the interface. The temperature field around the bubble in the variables proposed in the work becomes stationary, and the solution to the thermal problem becomes self-similar. The dependence of the bubble radius on time takes the root form, and the proportionality coefficient depends only on the Jacob number (which corresponds to the well-known solutions of other authors). However, as calculations show, for certain operating parameters of the process (in particular, for high Jacob numbers), this stage is unattainable during the real bubble growth processes.
Thus, a relatively simple semi-analytical solution to the problem of the vapor bubble growth in an initially uniformly superheated liquid is found. It is constructed on the basis of the existence of a quasi-stationary state for the bubble growth process. This allows the original boundary value problem with a moving boundary to be reduced to a system of ordinary differential equations of the first order. It should be noted that this solution accounts simultaneously for both the inertial and thermal effects that control the process. It can become the better alternative to direct numerical calculations, since the use of the found analytical solution removes the need for resource-intensive numerical simulation of the heat problem. License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.