Visualization of Coalescence of Multiple Small Bubbles with Closed B-Spline Curve

The smoothness of surface shape is one of the key issues to simulate coalescence of underwater bubbles. In this paper, B-spline closed curve is used to realize the visual simulation of multibubble coalescence. The main idea of the proposed algorithm is to construct a continuous bubble deformation which is guided by the normal direction of each control point and weighted by the distance from the point to the geometry center of the contour. The advantages of this algorithm include the smoothness of the bubble contour in the processing of deformation and the randomness of dynamic process and coalescence process. Experiment results show that the simulation algorithm works well and can be used in 3D computer games and animations.


Introduction
The bubble phenomenon is a common natural phenomenon that occurs in carbonated drinks, waste water, the breath of fish in the water, etc. It is widely considered in various disciplines and industries. Here are three classical applications and literatures.
In the petrochemical industry, during oil production and transmission, artificial injection bubbles can carry heavy oil from the ground to the surface [4]; the secondarymetallurgical treatment of the liquid steel relies on the injection of purge gas for improving the steel cleanliness [19]. In addition, liquid gas reactors depend on the interaction of bubbles to change the contact area of different phases [16].
In materials science, due to the different mechanical properties of nitrogen bubbles with different diameters on metals, the mechanical effects of nitrogen bubbles and their polymerization behavior on HCP rep-Zr alloy were studied [24]. The stability of bubbles and the potential effects of their interactions are also studied [30]. The interaction between bubbles in liquid affects the mass transfer performance of bubbles [2].
In movies and games, bubble simulation is a common material, as shown in Figure 1. It enhances the user experience in entertainment systems [5,25,29] and the reality visualization effect [33]. Kakuta et al. proposed a mixed reality application for enjoying playing with soap bubbles [17]. Jeong-Mo applied bubble research to animation and game creation based on bubble simulation [15]. Urikovic took a complete computer simulation of a soap bubble from a dynamics perspective and discusses the animation of bubbles [31]. Kunimatsu et al. used computer graphics to animate the phenomenon of throwing the jar into the sink, and the resulting individual bubbles move through the tank [22].
As we often encounter multiple bubbles or groups of bubbles in these practices, it is important to fully understand the interaction and binding of bubbles in a fluid [13,32]. Coalescence is one key behavior of bubbles and simulated in many literatures [6][7][8]14]. These literatures focus on the effects of initial bubble spacing, size and arrangement angle on the coalescence process, the changes in contour [20], the rise velocity and shape of bubbles in pure water [9], the pressure field around the bubble [7,14], and so on.
In general, coalescence of two bubbles takes place in three main steps: bubble collision and deformation, film drainage, and film rupture [19]. Kee and Chhabra [18] observed the bubble shape and coalescence behavior by using high-speed camera technology and investigate the effects of rheological properties and surface tension on bubble properties and coalescence. Li et al. [28] numerically simulates the interaction of bubbles online to obtain the effect of dimensionless distance on the coalescence of bubbles.
For a long time, many scholars have conducted extensive research on bubble behavior [26] and established a relatively sound theoretical basis. However, existing methods mainly study behaviors of one or two bubbles and pay a little attention on the smoothness of the contour of bubbles after coalescence. In the view of computer games, visual effects and the game experience in some cases are more important than physics demands [29]. We in this paper emphasize on the smoothness of the contour and the visualization of multibubble coalescence phenomenon, which results in solving issues as follows: (i) Multibubble coalescence rules and deformation rule are built to guide the implementation of the visualization algorithm (ii) The contour in 2D and the surface in 3D of a bubble should be C1-continuous, so the smoothness of the bubble contour in the processing of deformation is kept (iii) The multibubble coalescence algorithm is developed

Method
For visualizing the process of multibubble coalescence, the surface model, the continuous deformation, and rules of motion should be addressed in advance.

Modeling Surface of a Bubble.
For the case of small size bubble, its shape can be represented by a sphere [3] or an ellipsoid [1,12]. However, in most cases, the contour of a bubble is not a simple sphere or a circle, and it can be approximated by a Fourier series [10,11] or other functions. Given the coalescence of two bubbles from any directions, methods that use x-axis or y-axis as symmetry axis to model contour and surface are not ideal. Given the closed B-spline shows the advantage of the design deformability and flexibility [35], we here use the closed B-spline method to model the contour of a bubble. Let P i , i = 1, 2, ⋯, n, be the control points which denote a sequence of ordered data points on the bubble curve. Let denote a B-Spline curve to be determined to approximate the bubble curve. In formula (1), N i,p ðuÞ are the B-spline basis functions of a certain order and can be expressed as follows: We illustrate the calculation of the contour after coalescence of two bubbles in two-dimensional (2D) case. Let a bubble be displayed by a circle or a set of contour points which are ordered and represent discrete closed curve before coalescence. When two bubbles are too closed, the coalescence occurs due to interbubble attraction, as shown in Figure 2(a). For obtaining the new contour of the bubble, we at first find the intersection points (points A and B, denoted by red points) of two curves with Little's method [23]. Then, those inner points (P 1 , P 2 , Q 1 , and Q 2 ) are detected with Kpalma's method [21] and removed. Finally, we fit those contour points by the closed B-spline method, as shown in Figure 2 For building the 3D surface of bubble, the surface can be obtained as an axis-symmetric surface of rotation. In our experiments, the perpendicular bisector of the line segment AB in Figure 2(a) is used as the axis of symmetry. Figure 3(a) shows another contour with the axis of symmetry. Figure 3(b) is the surface that generated by axissymmetric surface of rotation using the contour points of Figure 3(a).

Continuous Deformation.
After coalescence, the contour of two coalesced bubbles will continuously deform over time due to the combined effects of multiple forces, including surface tension and water pressure. Then, a method of contour deformation guided by the normal direction and weighted by the distance from the point to the geometry center of the contour is proposed. Let P i , i = 1, 2, ⋯, n, be the control points located on the contour, n ! i be the unit normal direction of P i , and d i be the Euclidean distance from P i to the geometry center point C of the contour, as shown in Figure 4.
At first, we calculate the average distance d from contour to the bubble center. There are two cases. The first case occurs when two bubbles coalesce. where r 1 and r 2 are the radii of two bubbles before coalescence. Accordingly, if the simulating is experimented in 3D space, the average distance d from contour to the bubble center is calculated with the following formula.
Formulas (3) and (4) are, respectively, used for keeping the area and volume unchanged after coalescence.
The second case is for the deformation of one bubble that changes with time.
Then, the moving step τ i of the control point P i is where η is a positive number specified by user, and The control point P i moves to the new position P i ′, as When all control points move to the new positions, the contour can be updated with the closed B-spline fitting.
Note that: (1) The parameter η is used to adjust the speed of deformation (2) If the contour is a circle, then d max = d min , and the contour keeps the same shape. Therefore, the end shape of the contour is a circle or a sphere (3) In our experiment, the normal direction of control point is away from the center point C. If τ i > 0, i.e., d i > d, the control point P i moves towards to the center point C, and vice versa (4) These proposed formulas can be applied to a threedimensional (3D) coalescence simulation 2.3. Rules of Motion. In this paper, motion features of bubbles include rising speed, small random swing, and coalescence. In still water, the rising speed v t of a bubble is related to its volume V, the height H from the bubble to the water surface.
When both H and the change of bubble shape are small, f ðV, HÞ can be a linear function, which is easy to calculate and real time in computer games. So, the contour points rise according to the formula Adding small random swing to the rising bubbles can improve the effect of coalescence and reality. In our 3 International Journal of Computer Games Technology experiment, the direction of rising of a bubble is with random disturbance.
In summary, according to formulas (10), (11), and (8), the moving of a bubble is As for coalescence of two bubbles, if the bubbles are spherical [14], then we can use the bubble radii and the distance between sphere centers to determine whether the two bubbles are intersected or compatible; if the contour of at least one bubble is represented by discrete points, we use with Little's method [23] to determine whether the two bubbles are intersected and get the intersect points if there are.

Algorithm of Coalescence.
In this subsection, we summarize the algorithm for simulating the process of multiple bubble coalescence.
Let the number of bubbles is N, the coordinates of control points of the k th bubble be P ðkÞ i , where k = 1, 2, ⋯, N; i = 1, 2, ⋯, n k . The simulation time T is represented by the number of iterates I c , the still rising speed is v, and the random factor is σ.
Since the deformation of contour is guided by both the normal direction and the distance of contour point to the center of bubble, the proposed algorithm of bubble coalescence is named NDgBCA (Algorithm 1). According to Algorithm 1, the time complexity is OðN 2 × I c Þ, because we use enumerate method to detect the coalescence between any two bubbles in Step 4.

Experiments
All experiments are done on a laptop computer, with an Intel Core i7-4710MQ CPU@2.50 GHz processor and 4.0 G memory. Our algorithm is implemented in Matlab (2012 version) language. We test our algorithm (NDgBCA) with several experiments: two bubbles in the horizontal line, coalescence occurs from any direction, multiple bubble coalescence, and summary of the information of bubbles after coalescence.
3.1. Two Bubbles in the Horizontal Line. Many papers focus on the case, coalescence of two 2D bubbles in horizontal direction [11,14]. Bubbles are generated with mathematical method. In this experiment, r 1 = r 2 = 0:6, the centers of two circles are (-0.5, -1.8) and (0.5, -1.8). H = 40 is the final height of rising. v x = 0, which means that the bubble does not move in the horizontal direction. v y = H/f , where f is the number of simulation steps and also represents the number of iterates in simulation process. In this experiment, f = 10, and ten key frames captured from the simulation of the coalescence process are shown in Figure 5  Step 1: The current bubble number N ′ ← N.
Step 4: For each bubble, fit those contour points by the closed B-spline method [35].
Step 5: Detect whether there are coalescence occurring in the new positions.
Step 6: If there are coalescence, Step 7: Find the intersection points using Little's method [23].
where τ abides standard normal distribution N (0, 1), and σ is specified by user. In our experiment, σ = 0:8. The simulation results are shown in Figure 7(a). From bottom to top, 6 bubbles coalesce and become 5 bubbles at last.
For the case of variable rising speed, i.e., dv y /dt ≠ 0, referring that the trailing bubble has the higher acceleration [19], we assume that the rising speed v y is inverse to the distance from current position (Ho) to the top (H) of water surface, which can be calculated with formula (14). where parameter ζ is a positive number and specified by user. Under this assumption, the trailing bubble has a chance to be coalesced to the former bubble. We conduct a simulation experiment using 5 initial bubbles and the assumption. The simulation results are shown in Figure 7(b). From bottom to top, 5 bubbles coalesce and become 3 bubbles at last. The coalescence of two bubbles on the left illustrates the coalescence effect of formula (14).

Visualization with 3D Surface.
In order to enrich a real environment with virtual information, we can use axial symmetry technology to generate 3D surface of bubble ( Figure 3) and merge these 3D bubbles into an underwater image, as shown in Figure 8. From the bottom to the top in the figure, it can be found that after coalescence, the bubble gradually develops into an ellipsoid or a sphere, which is the same as the statement in other literatures [19]. From this example, it is obvious that our method can be used to add bubbles in some computer games and improve the visual effect.

Performance
Analysis. This subsection discusses the performance of the proposed algorithm (NDgBCA). All these data are illustrated with the data from the last experiment, i.e., Figure 8.
The model complexity is represented by the number of points of a curve or the number of polygons of a mesh surface. In this experiment, the number of bubbles is 5. For each bubble, its radius is in [0.40, 0.80], and the number of points that represent the contour is 90. After axial symmetry rotation, the numbers of vertices and polygons of each 3D bubble are 1622 and 3240, respectively.
The time complexity is represented by running time of each step of the coalescence processing. In Figure 9, the blue curve is the running time in each step. The average time of each step is 0.10 seconds. The first step is to generate 5 bubble contour points, and the other steps are corresponding to the proposed algorithm that detect coalescing or not and update the contour if coalescence cost only 5% of whole time. Most of the time is spent in calculating continuous deformation and generating 3D mesh bubbles. The orange curve is the cumulative time, and it shows a liner curve, which means that the algorithm is stable in the simulation iteration processing.

Discussion
The behavior of underwater bubbles is very complex. Our method emphasizes on visual effects of bubble coalescence and pays little attention on the physical theory. As mentioned in [29], creating the underwater feeling is a tricky balance between reality and game experience. So, our method is particularly suited for visualizing bubble coalescence in the underwater computer games.
Although particle system gets good performance when simulating small particle motion, fire and liquid [27], for example, literature [29] uses particle system to create bubbles uniformly; some parameters of particle systems are not open and cannot be dynamically adjusted during operation. Therefore, its performance is limited. Our method is easy to use because there are only a few parameters, including initial radii and the number and locations of bubbles, which are set at the beginning of the NDgBCA. In addition, our method also shows more details of bubble coalescence process and achieves better visualization effect.
Comparing to the method of using fitting curves [34] to visualize the bubble coalescence process in horizontal direction, the proposed method shows an advantage of displaying the coalescence process of two bubbles in any direction. If a bubble pair coalesces in the vertical direction, the coalescence effect visualized with our method is in line with the results from the experimental investigation of bubble coalescence [19], as shown in Figure 10, and numerical simulation [36]. Figure 10(a) displays their sequence of key frames acquired by a high-speed video camera and one of the bubble coalescence processes is marked by a yellow arrow [19]. Figure 10(b) shows our simulation results. In the simulation experiment, we adopt an assumption that for two bubbles on the vertical line, the trailing bubble is accelerated to  approach the leading bubble, resulting in the coalescence of the two bubbles [36]. It is interesting that the phenomenon of droplets has a similar coalescing process [26], which means that the proposed method can be used to visualize droplet coalescence with minor modifications based on the force analysis of droplets.

Conclusions
In this paper, a novel algorithm (NDgBCA) is proposed for simulating the multiple bubble coalescence. In the simulating, the continuous bubble deformation process is guided by the normal direction of each control point and weighted by the distance from the geometry center of the contour. The advantages of this algorithm include two aspects. The one is that the contour surface is uniform and smooth. The other is that it can simulate the coalescence of bubbles of different sizes from any angle and direction.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.

Acknowledgments
This work is supported by the National Natural Science Foundation of China under grant no. 51672028. Time of one step Cumulative time Step Time (s) Figure 9: Time consumption in each step of the coalescence processing.
(a) Image from literatures [19], and coalescence events are marked by yellow arrows