## 3 Numerical Results and Discussion
### Rayleigh-Taylor Instability  
The Rayleigh-Taylor instability is of great importance for multiphase flow studies and thus it is studied in this section. The computational zone is a rectangular box of \(\Omega=[-d/2,d/2]\times[-2d,2d]\) (\(d\) is the width of the box), where the heavier fluid is set in the upper location and the lighter fluid is set in the lower location. Initially, there is a small disturbance at the interface between the two fluids. The position disturbance function is \(y(x)=0.1d\ cos(2\pi x/d)\) in the beginning. The periodic boundary condition was adopted  
Figure 13: Development of the interface for bubbles merging. at the left and right walls and the no-slip boundary condition was adopted at the top and bottom boundaries. The study is governed by two non-dimensional parameters: the Reynolds number and the Atwood number, \(\text{Re}=\frac{d\sqrt{gd}}{\
u}\) and \(At=\frac{\rho_{H}-\rho_{L}}{\rho_{H}+\rho_{L}}\), where \(g\) is the gravitational acceleration, \(\
u_{H}=\
u_{L}=\
u\). These parameters were set for the study: \(d=200\), \(\sqrt{gd}=0.01\), \(\rho_{L}=1\), \(\zeta=3\), and \(At=0.5\), and the conditions \(Re=256\) and \(3000\) were studied. The parameters of the coefficient \(k=k_{1}\Gamma\Big{(}\frac{\Delta x}{L}\Big{)}^{2}\) in Equation (11) were set as \(L=d\) and \(k_{1}=1\) or \(k_{1}=10\) to observe the influence of the correlation coefficient. For the case of \(Re=256\), we set \(k_{1}=10\). Figure 14 compares the interface positions and the velocities of the rising bubble of light fluid and the falling spike of heavy fluid with the previous studies. They agree well with the references [10,22]. Figure 15 shows the evolution of the total mass of the lighter phase in the present method and the original DI method. The mass obtained with the original method is not stable compared with the present method. For the case of \(Re=3000\), we set \(k_{1}=1\) or \(k_{1}=10\) to explore the influence of the correlation coefficient on the interface. First, we set \(k_{1}=10\) and compared the interface  
Figure 14: Interface positions and velocities of the bubble and spike for the Rayleigh–Taylor instability at \(\text{Re}=256\) compared with references [10,22]: (**a**) interface positions, (**b**) interface velocities. Figure 15: Variation in the mass for the lighter fluid at \(\text{Re}=256\). positions with the previous studies [10; 45], as shown in Figure 16. The present results agree well with the published data. Figure 17 shows the development of the interface for the Rayleigh-Taylor instability of \(Re=3000\) with the present method at \(k_{1}=10\). The production of the falling spike of the heavier fluid and the rising bubble of the lighter fluid can be clearly observed. At first, the heavier fluid close to the falling spike produces some small droplets, and then the small droplets break away from the spike. As time goes by, the droplets sloughed off will generate some smaller droplets and slough them off as before. The phenomena verify the Kelvin-Helmholtz instability. Finally, Figure 18 shows the interface shapes of the original method and the present method with \(k_{1}=1\) and \(k_{1}=10\) at \(T=5\). For the interface shape obtained with the original DI method, the interfacial dispersion is quite strong. The interface of the droplet sloughed off from the heavy fluid is fuzzy. For the case of \(k_{1}=1\), the interfacial dispersion is weaker than the original method. However, the phenomenon is not obvious, and the bands between the droplets can also be found. For the case of \(k_{1}=10\), the interfacial dispersion is suppressed significantly. The interface shape is much clearer than those of the two other cases and the bands between the droplets almost disappear. It is illustrated that the coefficient of the correction term has a great influence on the interface and an appropriate coefficient determined by the flow parameters can suppress the interfacial dispersion well.