Several numerical methods of two-fluid flow with a moving interface have been posed.
The methods which are used to predict the phenomenon of two-phase flow can be divided into two main categories: Lagrangian and Eulerian schemes (Fig. 1.1) [6].
1.2.1 Lagranian schemes
The first method can keep free surface sharp between the two fluids and present the exact position of the free surface with re-meshing as the calculation proceeds. The mesh of this scheme is deformed and changed all the time. The main procedure of this method is that the position of the free surface at next time step is calculated by using the velocity field which is known. When the new free surface boundary is defined, we can reconstruct grids and update new properties for the new flow field [7]. The position of the interface can be predicted precisely, because the boundary mesh matches the free surface. Although the accurate prediction of the free surface can be carried out by the Lagrangian scheme, this method can not be employed to the quite complex flow field. Many deformations and stretches which result from the breaking, overturning, or gravity wave may cause numerical errors as well as reducing the precision. In [8], there is a Langranian scheme which is presented by Peric et al. This scheme can be employed in some simple and no large interface deforming cases, but its drawback is that the scheme cannot be used while the interface is unfolded.
1.2.2 Eulerian schemes
The second method can reduce errors which result from the deforming grids by using fixed grids that are generated before calculating the movement of the interface. The main
disadvantage of this method is the fact that it is prone to result in numerical diffusion. The diffusion will make the interface to spread over several mesh cells and the interface between the two immiscible fluids is going to be no longer sharp. In reality, the interface remains sharp due to the surface tension and the action of gravity, which separates immiscible fluids of different densities [2]. The interface of the two-phase flow must be tracked by employing some special treatments because its motion can not match the mesh any more as the calculation proceeds. A lot of techniques have been developed to cope with the problems of the multi-fluid flow systems in the past decade. These techniques can be classified to three main categories: 1. tracking the interface by using a set of mass-less particles; 2. using several mass-less maker particles to point out the only one kind of fluids and interface; 3. capturing the interface by a indicator function, such as a level set function or a volume fraction function.
There are several Eulerian schemes will be introduced in the following part.
(A) Front tracking method
This front tracking method [9] (Fig. 1.2) is applied to construct the interface between the liquid and gas by a simple trajectory technique. A lot of mass-less particles are uniformly distributed over the interface in the first instance, but the Navier-Stokes equations are solved in a fixed and Eulerian grid system. The numbers of particles on the interface may be increased or reduced as the calculation proceeds. The new positions of these particles can be obtained by integrating the Eulerian fluid velocity field near the particles for each time step.
This method has been used to deal with the motion of rising bubble, the breaking of water waves, and the collapse of an unsupported water column. This method is quite accurate but rather complex. Its first drawback is that the re-meshing of the Lagranian mesh is needed.
Another difficult is that transforming the mesh data of the Lagranian system into the Eulerian is quite complicated. In the three-dimensional problems, the front tracking method has another problem because the particles on the interface are not a string one any more. This problem will cause the calculational time and the computer storage to increase significantly.
(B) Level set method
A level set method for moving interfaces was proposed in [10]. The interface is identified as the zero level set of a smooth distance function from the front of the interface. This method not only eliminates the problems of the numerical diffusion which will smear the sharp front, but also avoids adding or reducing points to the moving grid. This method presents the interface by solving a scalar convection equation of the level set function. This method is easy to code due to the use of Eulerian grid and can result in more accurate results when the flow motion of the interface coincides with one of the coordinate axis. This method can also be easily generalized to three dimensions. However, the main drawback of this method is that level set methods loses its accuracy because the mass is not conserved when the interface is significantly deformed. Sussman et al. used it to simulate the flows of bubbles and droplets [10], and Li presented the results of Rayleigh-Taylor instability [11].
(C) Marker and cell method
Marker and cell method (MAC) (Fig 1.3) was proposed by Harlow and Welch in [12].
Several mass-less marker particles are distributed over a space which is filled with one particular fluid with a free surface, and these marker particles are used to calculate the motion of the flow field including the free surface. This method is quite accurate and can be used accurately to deal with many complex problems, such as an interface subjected to shearing and vorticity, and wave breaking in two-dimensional system, but it may become expensive of operating in three-dimensional one. More marker particles will be added when treating problems with interface stretching, shrinking, breaking, or merging in three dimensions. The above process results in increase of computational time and computer storage. There have been many studies about this method [13, 14].
(D) Volume-of-fluid method
In volume-of-fluid method, the fluids of two-fluid flow are represented by one scalar indicator function called volume fraction. The value of the volume fraction is bounded
between zero and unity. The value of unity denotes one of the fluids. The volume of zero denotes the other fluid. The volume fraction value between zero and unity indicate the interface. This method is quite popular and easy to code in the finite-volume method. The scalar indicator function is convected through the computational domain by solving a scalar convective equation like other transport equations. The scalar indicator function can not maintain a step function on the interface because most convective schemes result in numerical diffusion and dispersion. There are three categories of this volume-of-fluid method as follows.
Line techniques
This method has been implemented in two-dimensional problems, but the reconstruction of the interface in three-dimensional flows is difficult. The methods are used for interface reconstruction can be classified into three categories as follows.
The first method is SLIC method (Simple Line Interface Calculation) which was proposed by Noh and Woodward in 1976 [15]. The interface is approximated by using lines parallel to one of the coordinate axes. The volume fractions of the left and right cells of the prime cell are used to reconstruct the interface in the prime cell approximately when the sweeping direction coincides with the x-axis. On the other hand, the volume fractions in the cells above and under the prime cell are used in the y-axial sweeping. The second method is the one with improvement on the SLIC method by Chorin [16]. All direct neighbors of the prime cell will be used for interface reconstruction in the prime cell. The third method (PLIC or Youngs’ VOF) which is posed by Youngs [17] is more accurate than the SLIC method. In this method, the interface is approximated by using oblique lines. Unlike the SLIC method, all neighboring cells are used to approximate the slope of the interface in Youngs’ VOF. The above methods are illustrated in figure 1.4.
Donor-acceptor techniques
In this method, the value of the volume fraction transported through a cell face between two cells can be approximated by the volume fraction value of the downwind cell (Figure 1.5).
This method will cause the volume fraction values unbounded, i.e. the values of the volume fraction may become greater than one or less than zero. In order to ensure the boundedness, the method which improves the level of volume fraction value on the face by using the value of the donor cell is proposed by Ramshaw and Trapp [18] , but it will cause the incorrect steeping on the interface due to change any finite gradient into a step. As mention above, the volume fraction on the face correlates closely with the flow and interface direction. Another method was proposed to cope with the above problem in [19]. Hirt and Nichols calculate the volume fraction value on the face by including some information on the slope of the interface into fluxing algorithm.
High-order differencing schemes
In this method, the convective scalar transport equation is discretised by using a high order scheme or a blending scheme to predict the interface of two-fluid system. The main errors of this method are numerical diffusion which smears the front of the fluids and numerical dispersion which causes non-physical oscillation. The first-order upwind scheme has the numerical diffusion. This diffusion becomes significantly strong when the flow direction is normal to the interface direction. In order to reduce the numerical diffusion, the linear upwind scheme (LUDS) [20] and the quadratic upstream interpolation for convective kinematics (QUICK) scheme [21] were proposed. The former is second-order accurate and interpolated by the two upwind values. The latter is third-order accurate and interpolated by the two upwind and one downwind values. These high order schemes can reduce numerical diffusion, but they may cause numerical dispersion, such as oscillations, in the strong gradient regions.
In order to cope with the dispersion problem, the flux-blending and flux-limiter technique have been proposed. The former can be classified into two classes. The first class is based on adding an anti-diffusion flux to a first order upwind scheme [22] and used to resolve sharp gradient without over-/under-shoots. The second method is based on introducing some
smoothing diffusive fluxes into an unbounded high-order scheme, and it can prevent oscillations. The flux-blending technique will become expansive due to their multi-step nature and balancing the two fluxes. In [22], the flux corrected transport (FCT) method has posed.
FCT schemes are non-diffusive in nature, but create unphysical flotsam and jetsam.
The flux-limiter technique can remove non-physical oscillations and is based on the numerical flux on the interface of a cell which can be adjusted by using the flux-limiter function that enforce the boundedness. High resolution schemes (HRs) [23] are the schemes which obey the above criterion. The methods, such as Normalized Variable (NV) and Normalized Variable Diagram (NVD) [24], can be used to employ the flux-limiter technique.
The flux limiter function is presented by Van Leer [25]. Sweby developed the Total Variation Diminishing (TVD) [26] approach for high resolution schemes. In the past decades, many high resolution schemes have been proposed, such as SMART of Gaskell and Lau [27], GAMMA of Jasak [28], SUPERBEE of Roe [29], STOIC of Darwish [30], MUSCL and Van Leer of Van Leer [31].
Numerical diffusion can be classified into two main components, namely cross-stream and stream-wise. These two numerical diffusions can be associated with the angle between the flow and interface direction. The blending strategy was proposed in order to improve the accuracy and less numerical diffusion including cross-stream and stream-wise. The key issue in the composite scheme is not just when to switch, but how to switch [32]. Hence, the best approach must have a continuous switching function whereby the values of compressive and high resolution schemes are blended together with a blending factor. This method has been used in utilized in the HRIC of Muzaferija [33], STACS of Darwish [34], CICSAM of Ubbink [35] and the composite of MUSCL and SUPERBEE [36].