## 3 Numerical Results and Discussion
### The Evolution of a Square Bubble  
In nature, an irregular bubble immersed in another fluid changes to a circular one under the action of the surface tension without considering gravity. In this section, we study the deformation of a static square bubble. To mimic this process, the mesh size of \(301\times 301\) was adopted for simulation. The natural boundary condition was applied to all walls. The following parameters were applied for simulations: \(\rho_{H}=10\), \(\rho_{L}=1\), \(\
u_{H}=\
u_{L}=0.002\), \(\sigma\) as defined as \(0.005\), and \(\xi=3\); in order to measure the size of the bubble, a parameter \(R=D/2\) was set, where \(D\) is the length of the side for the square bubble. Cases of R = 20, 40, and 60 were simulated. The parameters of the coefficient \(k=k_{1}\Gamma\Big{(}\frac{\Delta x}{L}\Big{)}^{2}\) in Equation (11) were set as \(k_{1}=1\) and \(L=2R\).  
Figure 8 shows the variation in the mass for the square bubble obtained with the two methods with different radii for the case of \(\sigma=0.005\). It is obvious that the present method had better behavior in mass conservation. In contrast, the mass obtained with the original method decreases and the mass loss increases as the radius grows, which is consistent with the previous results of verifying the Laplace law. Figure 9 shows the results for the final shape with the two methods and compares them with the theoretical solutions for \(R=20\) at \(\sigma=0.005\). For the original DI method, as shown in the Figure 9, the change of the radius is \(9.2\%\) lower than the value of \(12.5\%\) in Figure 6. Figure 10 shows the contour of the order parameter along the centerline in the y direction of the whole zone of \(R=20\) obtained with the present method, which agrees well with the analytical solution. In Figure 11, the development of the interface for the bubble with \(R=20\) from a square shape to a circular shape is shown. The four corners of the square bubble are shrunk and the shape turns into a circle with the surface tension. All the phenomena verify that the present method has good reliability for mass conservation and the capability of maintaining interfacial shapes.  
parameters were adopted: \(\rho_{H}=10\), \(\rho_{L}=1\), \(\
u_{H}=\
u_{L}=0.002\), \(\sigma=0.005\), and \(\xi=4\). The parameters of the coefficient \(k=k_{1}\Gamma\left(\frac{\Delta x}{L}\right)^{2}\) in Equation (11) were set as \(k_{1}=1\) and \(L=2R\). By calculating the parameters, it is obvious that the bubbles will merge.  
Figure 12 shows the results for the final interface position obtained with the proposed method and with the original DI method and compares them with the theoretical solutions. It can be seen that the mass loss still exists with the calculation, while the result obtained with the present method is consistent with the theoretical solution. Figure 13 shows the development of the interface of two bubbles merging simulated with the proposed method at different times. It can be clearly observed that the two bubbles touch, merge, and produce a stable bubble. During the whole process, the interfacial shape is maintained well without seeing the non-physical phenomenon. The results indicate that the present method has good reliability in mass conservation and is effective in maintaining interfacial shapes.  
Figure 10: Contrast of the contour for the order parameter along the y-centerline with the analytical solution for \(\sigma=0.005\) at R = 20.  
Figure 9: Contrast of the final interface for the bubble with R = 20 at \(\sigma=0.005\).  
Figure 11: Development of the interface for the square bubble of R = 20.  
Figure 12: Comparison of the final interfaces for the merging bubbles.