Correlations and Asymptotic Behaviors of the Shape Parameters of Floating Bubbles Using an Improved Numerical Procedure

: An improved numerical procedure is used to present the correlations between the shape parameters and Bond numbers of floating bubbles for a wider range of Bond numbers ( 5 (cid:3400) 10 (cid:2879)(cid:2873) (cid:3407) 𝐵𝑜 (cid:3407) 5000 ) than the previously reported range of Bond numbers ( 0.003 (cid:3407) 𝐵𝑜 (cid:3407) 241 (cid:4667) , and their as ‐ ymptotic relations as Bo  0 and Bo  ∞ . The proposed method is proven to be more precise and robust than the conventional methods in comparison with previous numerical and experimental results. In addition, the profile of floating bubbles and the related parameters are presented for a wide range of bubble sizes. The shape parameters are divided into three distinct Bond number re ‐ gions, and are fitted with a fifth ‐ order polynomial as a function of Bond number on a log ‐ log scale for each region. The parameters show two asymptotes, which can be expressed in a simple power law. In addition, the dimensionless maximum depth of the floating bubble is obtained as 𝐻 (cid:3404) 0.7291015 when 𝐵𝑜 (cid:3404) 4.755563 . These correlations and asymptotic relations are expected to assist in the development of scale models of dynamic bubble ‐ related phenomena such as bubble bursting.


Introduction
Bubbles floating on a liquid can be seen in various places, and an understanding of the related phenomena is important in the natural sciences and industrial fields. For example, the bursting of floating bubbles in the ocean affects atmospheric circulation and marine pollution [1][2][3][4]. In a nuclear power plant accident, the bursting of bubbles in the reactor acts as a source of radionuclides [5,6]. Volcanic eruptions are caused by bubbles with sizes measured in meters [7]. In food engineering, bubbles are involved in the release of foams and flavors, e.g., champagne, espresso, and beer [8][9][10][11]. Recently, boiling heat transfer in a complex nanofluid is closely related to bubble nucleation and growth [12,13].
Bubbles generated in a liquid rise due to buoyancy, and then burst after reaching the liquid surface [14]. Bubbles floating on a liquid can remain in quasi-equilibrium for some time before bursting. During this time, the thin liquid film that defines the surface of the bubble gradually becomes thinner and unstable due to gravity, external disturbance, and impurities [15]. After bursting, the dynamics of the bubble cannot be simulated by the present numerical method, and the partial differential equations of the conservation laws must be solved using the computational fluid dynamics techniques such as the finite element method [16], the arbitrary Lagrangian-Eulerian moving mesh method [17], the twofluid two-phase flow model [18,19], and the adaptive multi-phase method [20]. As shown schematically in Figure 1, the shape of the floating bubbles in the quasi-equilibrium state can be divided into three curves, namely the cap, the cavity, and the meniscus. These three curves meet at a junction point . As the various dynamic phenomena that occur when a floating bubble bursts depend largely on its initial shape, the accurate prediction of this shape is important [21][22][23][24][25][26][27]. If the thickness and drainage of the liquid film can be neglected, then the shape of the bubble in the quasi-equilibrium state is determined by the Young-Laplace equation representing the balance between the forces of hydrostatic pressure and capillary action [28][29][30][31][32]. Floating bubbles vary in shape from spherical to hemispherical, depending on their size. Therefore, the Bond number, , is used as a key parameter. Bashforth and Adams [28] presented a systematic theoretical equation for hydrostatic pressure and capillary forces and obtained a series solution of drops of fluids. Subsequently, Toba [29] calculated the shapes of floating bubbles by using a semi-graphical method. At the same time, Princen [30] used the Table of Bashforth and Adams [28] to find the junction point and calculate the shapes of bubbles of various sizes. Medrow and Chao [31] numerically calculated the bubble shape for a broader range of the bubble sizes (which corresponds to a range of 0.003 241, computed based on the parameter values of [31]), and derived an analytic solution using a perturbation technique for tiny bubbles that are difficult to calculate numerically. Recently, Lhuissier and Villermaux [23] obtained the bubble shape via a method similar to that of Toba [29]; however, there is no detailed information on their numerical method and on the accuracy of their solution. Bartlett [32] used a similar method to Princen [30] and presented the shapes of bubbles for Bond numbers between 0.1 to 100. Cohen et al. [33] developed a model considering the film weight of giant soap bubbles and compared the results of numerical solutions and experiments. Shaw and Deike [34] studied the coalescence of floating bubbles for a wide range of Bond number, and characterized the evolution of the underwater neck and the surface bridge. Miguet et al. [35] investigated the life of a floating bubble by exploring the role of drainage and evaporation on film thinning. Despite a long period of study, the bubble sizes used for analysis have been limited, and there is no rigorous evaluation of the accuracy and robustness of the numerical methods used. As bubble shape parameters play an important role in the theory describing various dynamic phenomena related to floating bubbles [36], accurate prediction of the bubble parameter is crucial, especially in a wide range of Bond numbers. For instance, the lifetimes of floating bubbles are related to the surface area of the bubble film [15]. The jet velocity generated when the bubble bursts is related to the depth of the cavity, and the amount of aerosol emitted depends on the radius and height of the cap [23,27]. Koch et al. [5] used the correlation between the dimensionless cap area and the dimensionless diameter to calculate the number of film droplets generated in the bubble bursting process, and the correlation between the dimensionless base radius and the dimensionless diameter to obtain the critical film thickness. Teixeira et al. [37] employed a perturbation method for large bubbles to obtain analytical expressions for bubble parameters as a function of the Bond number, and compared their results with experimental data. Puthenveettil et al. [36] also derived approximate analytical solutions for bubble parameters over a moderate range of Bond number (0 1) and compared these with their experimental results. Moreover, to the best of the present author's knowledge, no study has yet been published detailing the correlations between the Bond number and various bubble parameters over an extensive range of Bond numbers.
In the present work, an improved numerical procedure is proposed to calculate the shapes of floating bubbles over a much wider range of Bond numbers than that examined in the previous studies. The proposed method is validated against the earlier numerical and experimental results to prove its accuracy and robustness. New correlations between the Bond number and several important bubble parameters, along with their asymptotic relations, are then provided.
The novelty of this study is as follows. A new method is developed by introducing specific functions for accurately finding the location of the connection points on the bubble surface regardless of the ODE solvers and their boundary points. An improved algorithm to find the junction point is developed, based on a novel function to determine the type of meniscus.

Theoretical Model
A slightly modified form of the theoretical model presented by Toba [29] and Medrow and Chao [31] is used in the present work. Thus, the shape of the floating bubble is described by the Young-Laplace equation: where Δ is the pressure difference across the interface, is the surface tension, and and are the principal radii of curvature. A half-section of an axisymmetric bubble floating on a liquid is presented, along with its associated parameters, in Figure 1, where ℎ 0 is the depth of the cavity below the undisturbed liquid surface, is the point where the meniscus meets the undisturbed surface, and and are the tangent angles of the cavity and meniscus curves, respectively.
If the thickness of the cap film is ignored, then the cap has a spherical shape with a radius , and Equation (1) can be written as where Δ denotes the pressure difference across the cap film. Furthermore, if the weight of gas inside the bubble is neglected, then the cavity satisfies where and denote the liquid density and gravitational acceleration, respectively, Similarly, the meniscus shape is given by 1 1 .
As the bubble shape is a surface of revolution about the z-axis, the principal radii of curvature are given as 1 sin , 1 d sin d .
By using the capillary length, 2 / , the bubble parameters are nondimensionalized as given by Substituting Equations (2) and (5)-(7) into Equation (3) yields the dimensionless ordinary differential equations (ODEs) for the cavity shape, given here as where the parameter represents the dimensionless pressure difference across the bottom of the bubble, as defined by 4 2 .
Equations (8) and (9) The second initial condition of Equation (9) is simply given as Similarly, the relation is used to obtain the ODEs for the meniscus shape as follows: The boundary conditions of Equations (15) and (16) at are given as sin 0 and 0 At the junction point , the geometrical relationship given in Equation (18) must be satisfied: From Equations (10) and (18), the relation between and is given as 2 sin 2 .
As the cap is spherical with a radius , its shape is given by cos cos .
The volume of the bubble is computed analytically via integration by parts [38]. The volume of the cap is given by where 1 cos . The volume of the cavity is given by 1 .
The surface area of the cap is obtained as However, there is no analytic form for the surface area of the cavity. It should be numerically integrated as follows: where and respectively denote the X and Y positions on the cavity at /4.

The Bond number is defined by
, where denotes the equivalent spherical radius of the bubble.

Numerical Procedure
Given the parameter , the shape of the floating bubble can be determined. Since the bubble is symmetric about the z-axis, only half of it is needed. First, the shape of the cavity shape between the bottom and the junction point is calculated, then the shape of the meniscus between the point and the boundary point is determined. Finally, the shape of the circular cap is quickly obtained. The junction point between the cavity and the meniscus is not given analytically and should be found via a numerical iterative method. The accuracy of the bubble shape is significantly affected by the accuracy with which is located.
The dimensionless height is used to determine whether the guessed junction point ̃ matches the exact point . Given the position of , is readily calculated from Equation (19). Furthermore, after computing the meniscus shape, another height, , can be obtained. Thus, the bubble shape is considered to have converged when the error falls within a defined level of tolerance that indicates the accuracy of the numerical method used.

Improving the Accuracy of the Cavity Shape
The cavity shape can be obtained by integrating Equations (8) and (9) or Equations (11) and (12) numerically. However, previous studies have shown that it is necessary to divide the cavity into two [29,32,39] or three [31] parts in order to avoid divergence of Equations (8) and (9) at /2. In the present study, when the cavity was divided into two parts, the numerical error was found to increase rapidly as the bubble size decreased (see Figure 6). Hence, the cavity shape was divided into the following three parts or intervals: (I) 0 /4, (II) /4 3 /4, and (III) 3 /4 , and Equations (8) and (9) were applied to intervals (I) and (III), while Equations (11) and (12) were applied to interval (II). However, because is not an independent variable in these equations, the numerical integrations are not completely finished at /4 and 3 /4. Hence, a conventional method is to set an arbitrary far point (or ) as a boundary condition of the ODEs and to finish the calculation when is just greater than /4 (or 3 /4). Unfortunately, this leads to potential errors at the endpoint (or ) of the divided cavity depending on ODE solvers. The conventional method works well on medium-sized bubbles, but produces a severe error in large bubbles, especially at the position of , because the bottom of the cavity becomes substantially flat. This error is not only affected by the boundary points (or ), but also by the tolerance of the ODE solver.
To remedy this problem, an improved method was proposed for accurately finding the location of the connection points and regardless of the ODE solvers and their boundary points. In this approach, a monotonically increasing function ; is defined as follows: function ; Integrate Equations (8) and (9) from 0 to to find sin at the endpoint return sin sin end function An accurate position is then found by solving ; /4 0 in the interval 0 using a root-finding technique. Here, the boundary point is any value greater than ; usually, 100 is appropriate for most of the bubble sizes considered in this study. As shown in Figure 2b, does not affect the accuracy of the solution obtained using this method at all. Once has been determined, Equations (8) and (9) are integrated again from 0 to to obtain the first part of the cavity profile along with the remaining and . Similarly, the second part of the cavity and its endpoint , , can be computed using the following function: function ; Integrate Equations (11) and (12) from to to find cos at the endpoint return cos cos

end function
The computed bubble shapes according to the position of the boundary point for a large bubble ( 10 ) are presented in Figure 2. Here, the bubble shapes are seen to vary significantly depending on when the conventional technique is used (Figure 2a). When the root-finding method is used, however, the bubble shapes not affected by ( Figure 2b). Similarly, the effect of ODE solver tolerances for the same size of bubble are presented in Figure 3. Here, the bubble shape changes slightly according to the tolerance using the improved method (Figure 3b), but the results obtained using the conventional method are unacceptable (Figure 3a). Therefore, it is essential to compute the connection point accurately for large bubbles.

Finding the Junction Point Based on the Type of Meniscus
There are two general approaches to finding the junction point according to the direction of integration of the meniscus. The first approach is to integrate the meniscus from to in the negative X-direction [29,31]. Here, is sought in an iterative manner so that the derivatives of the meniscus and cavity profiles match at the junction point . However, it is difficult to compute the shape of the meniscus using this method because the initial slope becomes zero when the exact initial condition, 0 at , is applied. To circumvent this, Toba [29] and Freud and Freud [40] employed 0.00001 at as an initial condition. Medrow and Chao [31] further improved this method by using an approximate analytical solution of with the Bessel function. Nevertheless, as these methods do not give an exact initial condition at , the accuracy of the solution is still limited.
The second approach is to integrate the meniscus from to in the positive Xdirection [23,30,32]. Although this approach has the advantage of obtaining the position more precisely, there is another difficulty in that the meniscus shape changes dramatically around the junction point , as shown in Figure 4. Hence, a method for effectively locating the junction point is proposed herein by using a specific function that distinguishes the type of meniscus shape. A conventional bisection method for finding narrows the interval by comparing the sign of at both ends of the solution interval. At each iteration step, it is necessary to find the height , which is calculated differently depending on the meniscus type: if the meniscus profile is convex, its minimum value is chosen; if the meniscus profile is monotonically decreasing, an inflection point is selected. Interestingly, without calculating , the bisection method can become simpler by considering the meniscus type only. Hence, the following function is defined to determine whether a meniscus is convex at a guessed junction point:

function is_meniscus_convex
Integrate Equations (15) and (16) (15) and (16) from to to find sin at the endpoint return sin end function This modified bisection method is more efficient and robust than the conventional one because there are no expensive routines to find the minimum or inflection point of the meniscus. Although a similar meniscus type to that used herein was previously used by Princen [30], the method employed therein remained limited in terms of bubble size and accuracy because it relied on Bashforth and Adam's table [28].
The overall procedure for calculating a floating bubble shape is summarized in Algorithm 2.

Algorithm 2
The improved numerical procedure for calculating the shape of a floating bubble Input: Output: Floating bubble shape and parameters Calculate the first part of the cavity profile by integrating Equations (8) and (9) from 0 to /4 and by solving ; /4 0.
Calculate a part of the cavity profile ′ by integrating Equations (11) and (12) (20) and (21) from to . Combine the profiles to obtain the floating bubble shape: ⋃ ⋃ ⋃ .

Results and Discussion
In the present study, the ODE solver employs a variable step, the 4th and 5th-order Runge-Kutta algorithm using Dormand and Prince pairs [41]. A range of the Bond numbers, i.e., 5 10 5000, is applied to the improved numerical method. Accordingly, the results obtained using the conventional and improved methods are compared in the following subsection. Then the present results are compared with those of previous computational and experimental studies.

Comparison between the Conventional and Improved Methods
The conventional and improved methods are compared for large bubbles in Figure  5. Here, it can be seen that a smaller value of corresponds to a larger bubble size ( Figure  5a). Thus, when 1E-7, the bubble shapes obtained using the two different methods are almost identical. However, when 1E-14, there is a slight difference between the two results, and when 1E-22, the two methods show a significant discrepancy in the obtained bubble shapes. The cause of inaccuracy in the conventional method is demonstrated in Figure 5b, where the dimensionless cap radius is seen to increase abnormally as the parameter decreases in the conventional method.
(a) (b) The conventional and improved methods are compared for small bubbles in Figure  6. Here, when 1000, the bubble shapes obtained using both methods are almost identical ( Figure 6a). However, when 2000 and 3000, the conventional method raises the bubble in a non-physical manner because the bisection method does not converge and gives a negative value of . As shown in Figure 6b, when is small, the dimensionless absolute height | | is the same in both methods; however, as becomes larger, | | fluctuates significantly in the conventional method. These results clearly demonstrate that the improved method should be used to calculate large or small bubbles robustly.

Comparison of the Present Results with the Earlier Numerical and Experimental Results
The present computational results are compared with the previously published numerical results in Table 1. Various bubble parameters were previously calculated in the range of 1 20, 0.28 8, and 4 10 50 by Toba [29], Princen [30], and Medrow and Chao [31], respectively. However, Princen [30] used the parameter in place of , where the relationship between these is given by where denotes the dimensionless parameter of the present work, and / denotes the dimensionless parameters of Princen's work. The computed parameters in Table 1 agree well with one another, with the exception of Toba's parameters. This is because Toba employed a semi-graphical technique of limited accuracy. The earlier methods used an approximate analytical solution for tiny bubbles [31] or an interpolation method with a pre-calculated Table [30], so the sizes considered were narrow. Medrow and Chao [31] used the boundary conditions of 0.00001 and sin 0.0001 . Princen [30] searched the junction point until both Δ and Δ were less than 0.5 10 . In contrast, the present calculations maintain an accuracy of |sin | 5 10 and 5 10 , thus providing much more precise results. The experimental results of Teixeira et al. [37] are compared with the present numerical results for selected bubble parameters in Figure 7a, where the parameters are plotted with respect to the Bond number based on the cap radius, i.e., / . Here, the present computations are seen to agree well with the experimental results, except for the * values. For this parameter, the numerical results obtained using Surface Evolver S/W [42] have shown a similar trend to the present computations, as described in [37]. Meanwhile, the experimental results of Puthenveettil et al. [36] are compared with the present numerical results for selected parameters in Figure 7b, where denotes the height of the bubble top above the undisturbed surface. All parameters are normalized with respect to the equivalent bubble radius . The experimental results agree well with the present computation, except that the present computation overestimates the experimental data for / . Notably, the numerical results of Puthenveettil et al. [36] also overestimate / by a similar amount.
(a) (b) Figure 7. Comparison of the present numerical results with the experimental results of (a) Teixeira et al. [37] and (b) Puthenveettil et al. [36] for selected bubble parameters, where * √2 .
The presently computed profiles for bubbles of various sizes (red outlines) are compared with the experimental photographs obtained by Teixeira et al. [37] and Puthenveettil et al. [36] in Figure 8. In addition, the present predictions for the angle of the junction point ( are compared with those measured by Teixeira et al. [37] in Table 2. Teixeira et al. [37] used a soap solution ( 0.0282 N/m), and the images were taken in the range of 1 to 32.5 mm (Figure 8a). Here, the calculated bubble shapes are seen to match the experimental images well, although the results in Table 2 indicate slightly different values for the angle of the junction point ( ). Specifically, the present computational work predicts to be about 1 2 degrees higher than that measured by Teixeira et al. [37] Meanwhile, the experimental investigations of Puthenveettil et al. [36] were performed on various liquids (e.g., water and ethanol), and the photographs were taken for 0.004 to 2.46 (Figure 8b). However, when the Bond numbers measured by Puthenveettil et al. [36] were used in the present computations, the resulting numerical profiles for large bubbles were quite different from the experimental photographs. Hence, the values used in the present work were adjusted in order to provide a better match between the predicted and measured profiles. These adjusted Bond numbers are compared with those measured from Figure 8b in Table 3. Here, similar Bond numbers are predicted and measured for the two tiny bubbles, but the present computations predict much higher Bond numbers than those measured for the two large bubbles. While Teixeira et al. [37] measured the bubble cap radius and the junction angles to obtain an alternative form of Bond number ( ), Puthenveettil et al. [36] calculated by measuring the volume of a bubble from the experimental image. It is expected that obtaining is more complicated than obtaining because the lower part of the floating bubble is immersed in the liquid, whereas the upper part is not. Therefore, the discrepancy in Table 3 is thought to be due to experimental errors or the limitations of the current theoretical model.  Table 2. The angle of the junction point ( ) as measured experimentally by Teixeira et al. [37] and as predicted in the present work. Figure 8a Teixeira et al. [ Table 3. The Bond numbers ( ) measured by Puthenveettil et al. [36] and those of the present computational work. Figure 8b Puthenveettil et al. [

Important Bubble Parameters and Profiles for a Wide Range of the Bond Numbers
The present numerical results for important bubble parameters in a wide range of Bond numbers (10 1000, corresponding to 8 10 5723) are listed in Table 4. The relative tolerance of the ODE solver is set to 2.22 10 . For the range of considered here, the present computations satisfy |sin | 10 and 7 10 . Thus, the present method computes the boundary condition of the meniscus more accurately and converges the position of the junction point more precisely than the earlier methods. As a result, the present method provides highly precise and robust numerical results for a wide range of bubble sizes. The calculated bubble shapes for 10 1000 are presented in Figure 9, where each shape is drawn at equal intervals on a log scale of . In the interval of 10 to 10 , the bubble exhibits an almost hemispherical shape, and the cavity is virtually flat (Figure 9a). However, in the interval 10 0.1, the shape becomes less hemispherical as the bottom of the bubble continues to descend into the liquid (Figure 9b). In the interval 0.1 10, the shape continues to change from hemispherical to that of a bisected rugby ball (Figure 9c). Interestingly, the depth is found to increase to a maximum value and subsequently decrease within this interval. In the interval 10 100, the bubbles have acquired a nearly spherical shape and are seen to protrude only slightly above the undisturbed liquid surface (Figure 9d)

Correlations and Asymptotic Behaviors of the Bubble Parameters
The log-log plots of the various bubble parameters with respect to the Bond number for the interval of 5 10 5000 are presented in Figure 10. Each plot can be divided into the following three regions: (I) 5 10 0.1, (II) 0.1 100, and (III) 100 5000. The parameters for each interval are fitted with a fifth-order polynomial on the log-log scale according to Equation (27): where Π denotes the dimensionless bubble parameters. The coefficients are presented in Table 5.  Interestingly, can be well expressed even with a single "bending function" for 5 10 5000, as given by Equation (28) The curves fitted by the three fifth-order polynomials and by the single bending function are compared in Figure 11, where both correlations are seen to fit the numerical data well. Figure 11. Comparison of the fitting results of the fifth-order polynomials and bending function for .
Two asymptotes appear in Figure 10 as the Bond number goes to zero or infinity. These were obtained numerically via the following simple power law: where the coefficients and are shown in Table 6.