A Modified Normalized Weighting Factor method for improving the efficiency of the blended high-resolution advection schemes in the context of multiphase flows

This work deals with a new methodology for the implementation of high-resolution (HR) schemes employed to advect the volume fraction in the volume of fluid (VOF) method, in which the numerical stability and convergence depend heavily on the numerical advection scheme and implementation method. The proposed method is based on the normalized weighting factor (NWF) method, which linearizes the normalized interpolation profile and rewrites the face value directly using the donor, acceptor, and upwind nodes. However, unlike the NWF, which is fully implicit and results in pentadiagonal linear systems, the new modified normalized weighting factor (MNWF) method only forms the implicit terms with the contribution of the donor and acceptor nodes, while the contribution of the upwind node explicitly forms part of the source term. Therefore, the method results in a tridiagonal linear system. The comparison of the new method with the deferred correction (DC), downwind weighting factor (DWF), and the RNWF methods shows that the MNWF requires about 5%–25% fewer iterations than DC and RNWF, and around 10%–85% less than DWF. Thus, a similar order of accuracy of the results can be o btained with less computational time.


Introduction
The volume of fluid (VOF) method of Hirt and Nichols (1981) is a well-established conservative method to solve multiphase flow problems. The VOF introduces an additional transport equation to advent a marker function called volume fraction to define the position of the interface between the fluids. The volume fraction must be updated every time that the fluids move, and the boundary between the different fluids changes position. Nevertheless, updating the marker function is critical for the success of the simulation of multiphase flows and also is not a trivial job due to the purely convective nature of the transport equation (Tryggvason et al., 2001).
An option to deal with this problem is using the blended High Resolution (HR) schemes also known as interface-capture schemes which combine a high order (HO) diffusive scheme, a compressive scheme, and the Convection Boundedness Criterion (CBC) ensuring that no oscillatory behavior is experienced in the solution and have relatively low numerical diffusion (Tryggvason et al., 2001;Moukalled et al., 2016). An example of such schemes is CICSAM (Ubbink and Issa, 1999). However, the direct introduction of the blended HR schemes into the discretized equation is not suitable because of their composite nature. Thus, some techniques initially created to implement the HR schemes in the momentum equation have been used to overcome this difficulty.
For instance, Meyer et al. (2016) implemented the blended interface capture scheme BICS (Wackers et al., 2011) with the deferred correction (DC) method of Rubin and Khosla (1977) to develop a new code for simulating free-surface flows around modern sailing yachts. For the DC method, the implicit terms of the discretized equation are based on the upwind scheme, whereas the difference between the BICS schemes and the upwind scheme is Vo Vol. 3, No. 3, 2021, 208-225 Experimental and Computational Multiphase  , L factors for the NWF formulation n+1, n, superscript to indicate the next, current, and last n−1 time-step f subscript that denotes variable approximated at the face of the control volume ζ convergence criterion considered as a source term. Although according to Darwish and Moukalled (1996), DC suffers from low convergence rates, whereas the general approach of Meyer et al. showed better performance than other codes. However, the influence of DC on the general approach was not studied. On the other hand, CICSAM and its modifications such as THOR (Hogg et al., 2006), MCICSAM-W (Wacławczyk et al., 2007), and MCICSAM-Z (Zhang et al., 2014) employ the Downwind Factor Method (DWF) method of Leonard and Mokhtari (1990). The DWF introduces an auxiliary factor that implicitly contains higher-order wide-stencil information, but its implementation involves only the adjacent upwind and downwind node values. So, this method is suitable for tridiagonal solvers. However, the coefficients obtained from a DWF implementation sometimes are not diagonally dominant; thus, the formulation is not stable for many flow configurations and requires substantial relaxation to achieve convergence. Despite the described problem, this method is still used commonly in the multiphase community.
Another technique that overcomes the shortcomings of the DWF method, but which is rarely applied in the context of multiphase flows is the full implicit Normalized Weighing Factor (NWF) method (Darwish and Moukalled, 1996). The NWF linearizes the normalized interpolation profiles and rewrites the face value directly using the central, upwind, and downwind nodes so that the method uses a pentadiagonal stencil, and the diagonal coefficient results are always positive. Consequently, the NWF is much more robust than the DWF and faster than DC methods (Darwish and Moukalled, 1996). Nevertheless, NWF is not frequently used because it requires the pentadiagonal matrix algorithm (PDMA) to solve the system of equations.
In 2018, a revision of the described normalized weighting factor (RNWF) method was presented by Chourushi (2018), which is applicable for tridiagonal equation solvers. This method relies on the final discretization of the normalized weighting factor method and removes the contribution of far-off nodal values from the diagonal coefficient. These terms are later added as a source term. According to the author, the RNWF is four times faster than DC and 1.3 times faster than NWF.
Because of the stability advantages and efficiency of the NWF formulation compared to DC and DWF, and the new possibility of using it with tridiagonal equation solvers, we tested the RNWF method in the context of multiphase fluids and found that the convergence rate of the RNWF is similar to the DC method in the case of multiphase flows and that the RNWF method tends to degenerate the interface slightly. We supposed that the problem lies in the introduction of two explicit terms in the source term, the value of the center point and the upwind point instead of only the upwind point as is suggested in the original NWF method. Our new idea is only to introduce the upwind value as a source term. This paper presents this new idea that we call Modified Normalized Weighting Factor (MNWF) method which we apply for the numerical implementation of six blended HR schemes: CICSAM, MCICSAM-W, MCICSAM-Z, HRIC (Muzaferija et al., 1998), FBICS (Tsui et al., 2009), and CUIBS (Patel and Natarajan, 2015). The implementations are realized on the in-house finite-volume flow solver FASTEST, based on a block-structured collocated grid arrangement. For the investigation, we consider four test cases: the slotted circle, the circle in shear fluid, the rising bubble, and the break-dam with an obstacle. The convergence rate is given by the total number of iterations required for convergence, and the accuracy of the results is analysed for each test case and compared with the obtained one using the DC, DWF, and RNWF methods.
In the following sections, the mathematics and results of the new approach are presented. First, the governing conservation equations of the multiphase dynamics are shown, as well as the discretization process of the volume fraction equation. Then, the DC, DWF, NWF, and RNWF methods are reviewed, and their disadvantages are pointed out. Finally, the modified normalized weighing factor (MNWF) method is described in detail and used to solve the test problems.

Governing equations
Two viscous, incompressible, and immiscible fluids are modelled as a single effective fluid whose behavior is described by the conservation transport equations of mass, momentum, and volume fraction: where u is the velocity vector, p the pressure field, α the volume fraction field of one of the fluids, and ρ, μ are the physical properties density and dynamic viscosity, respectively. In the momentum equation (Eq. (2)), the last two terms of the right-hand side represent the gravity forces and the force due to surface tension. Since the multiphase system of the two fluids is treated as a single fluid, the local density and dynamic viscosity can be evaluated as a function of the volume fraction using the following relations: where the subscripts 1 and 2 refer to each fluid and α denotes the presence (α = 1) or absence (α = 0) of the traced fluid. A volume fraction between 0 and 1 indicates the presence of a mixture, and the value of 0.5 defines the interface between the fluids. The three conservation equations are discretized in space with the second-order finite-volume method and in time with the second backward differentiation BDF2 time scheme (Ferziger and Perić, 2012). The discretization of the mass and momentum equations for multiphase flow systems follows the methodology detailed in Sauer (2000) and it is not shown here as it is outside the focus of this paper. However, the following sections will be dedicated to the discretization of the volume fraction transport equation.

Discretization of the volume fraction transport equation
In order to solve the volume fraction governing equation, the domain is subdivided into an arbitrary number of control volumes which form a structured grid domain. Then, Eq. (3) is integrated over each control volume, and the volume integral of the divergence of α is transformed to a surface integral applying the Gauss theorem. Hence, the integral form of the volume fraction transport equation results in where V is the volume and S is the surface of the control volume, and n is the normal vector to the surface. The volume integral is resolved using the BDF2 time scheme, and the surface integral is numerically approximated by the mid-point ruler. For an arbitrary control volume P Eq. (5) becomes where the superscripts n+1, n, and n-1 represented the values at the next, current, and last time-step respectively. The subscript f denotes the variable approximated at the center of each face of the control volume P, and C is the convective flux. After some arithmetic operation, the algebraic form of the discretized equation reads . A P and A F are defined according to the advection scheme used to approximate the face volume fraction 1 n f α + .

Blended high resolution (HR) schemes
The accuracy of the numerical solution of Eq. (7) depends on the proper estimation of the face volume fraction. This demands an advection scheme that should neither produce numerical diffusion nor unbounded values (Muzaferija et al., 1998). Over the last decades, blended advection schemes between a compressive and a diffusive high resolution (HR) scheme have been used to advect the volume fraction. The reason for this is, that the use of just the compressive schemes can cause an alignment of the fluid interface with the grid (Ubbink, 1997), or that the use of only the diffusive HR schemes deteriorates the accuracy when the flow is not orientated along a grid line due to the false diffusion (Darwish and Moukalled, 2006). The switching strategy depends on the angle θ f between the flow direction and the grid lines. This approach has been employed to develop several blending advection schemes, also known as interfacecapture schemes, for example, CICSAM (Ubbink and Issa, 1999), HRIC (Muzaferija et al., 1998), STACS (Darwish and Moukalled, 2006), and FBICS (Tsui et al., 2009). A high resolution scheme is a composite high-order scheme combined with the Convective Boundedness Criterion (CBC) of Gaskell and Lau (1988) to ensure that the interpolation profile at the cell face does not underflow or overflow the cell (Moukalled et al., 2016). Some examples of them are SUPERBEE, MUSCL, SMART, or STOIC. The HR schemes can be formulated in the framework of the Normalized Variable Diagram (NVD) introduced by Leonard (1991), which is illustrated in Fig. 1. The UD line refers to the upwind differencing scheme, DD to the downwind differencing scheme, and the shaded area indicates the part of the NVD that fulfils the CBC. The schemes close to the UD line are linked with numerical diffusion but always produce a bounded solution and are stable, whereas the schemes near the DD line are unstable but introduce a negative numerical diffusion, so they are known as compressive schemes.
For the NVD, the normalized volume fraction α  in the vicinity of the cell P is defined as where the subscripts D, A, and U refer to the donor, acceptor, and upwind cells designated depending on the flux direction (see Fig. 2), and r is the position vector. With this normalization relation U A α α , 0 1 , = =   and the normalized volume fraction at the cell face f α  becomes a function of D α  . Then, a blended HR scheme designed within the NVD framework defines the normalized face volume fraction as (blended HR) (compressive) (HR scheme) (1 ) where the blending function λ = f(θ f ) varies between 0 and 1. Due to its composite nature (blended HR) f α  cannot be directly expressed in terms of the nodal values of the control volume P and neighbors F, which is necessary to determine the A P and A F coefficients and to solve Eq. (7) for the unknown values at the central nodes. So, several methodologies have been developed for the numerical implementation of the blended HR schemes. The following section provides a brief introduction of some commonly used methods, and the next section explains in detail the proposed new alternative.

Deferred corrector approach
The deferred corrector (DC) method of Rubin and Khosla (1977) is a simple technique designed for tridiagonal matrix solvers that define the convective term of the volume fraction equation as The first term, the value obtained by the upwind scheme, is used to form the coefficient matrix A for the nodal algebraic equation, while the second term, the difference between the blended HR and UD schemes is explicitly computed using the last available values and added to the source term. Thus, the resulting coefficient matrix is always diagonally dominant, yielding a numerically stable method. Nevertheless, the convergence rate decreases as the difference between the cell face value estimated by the upwind scheme and the blended HR scheme increases (Moukalled et al., 2016).

Downwind Weighting Factor method
The Downwind Weighing Factor (DWF) method was developed by Leonard and Mokhtari (1990) to overcome the low convergence rate associated with the DC method. The face value is defined as a weighted average between the donor and acceptor cell written as ( ) where the DWF is the weighting factor that varies between 0 and 1 and is explicitly computed as blended HR Then, the convective term using nodal values and considering the flow direction takes the form which is used to generate the coefficient matrix A. The effect of this weighting formulation is a reduced stencil for the discretized coefficients, which allows the system of equations to be solved with a tridiagonal solver. However, the diagonal coefficient A P becomes negative when 0.5( ) which is a common scenario for all HR schemes when 0.5 D α > . Consequently, this system of equations leads to unphysical results for many flow configurations and requires substantial relaxation to avoid convergence problems (Darwish and Moukalled, 1996).

Normalized Weighting Factor method
The Normalized Weighting Factor (NWF) method of Darwish and Moukalled (1996) describes the normalized face value as a linear function of the normalized donor value write as where  represents the slope and m the intercept of each linear function that is part of the HR scheme employed. Then, this linear relation is rewritten as The convective term using the nodal values is (17) The last expression allows the full implicit treatment of the HR schemes, so the values of the far nodes U α + and U αare actual nodes in the computational domain that can be resolved in the algebraic equation. For the one-dimensional structured grid, as shown in Fig. 2, the NWF form of the algebraic equation becomes A P(time) and b P remain as defined in the last part of Section 3. The above coefficients form a diagonally dominant matrix in which A P is always greater than zero because for almost all HR schemes,  is greater than m. Only for the DD scheme [ , ] [0,1] m =  , A P becomes zero. In this case, [ , ] m  is set to [ where L is the value of  from the previous interval of the composite scheme.
The NWF formulation is more robust than the DWF and faster than the DC method (Moukalled et al., 2016;Chourushi, 2018). However, the discretization above involves a pentagonal stencil that includes the far nodes in each coordinate direction which result in a pentadiagonal coefficient matrix.

Reviewed Normalized Weighting Factor (RNWF) method
The Reviewed Normalized Weighting Factor method developed by Chourushi (2018) is a tridiagonal version of the previously described NWF method. From the discrete equation (18), the change from pentadiagonal to tridiagonal stencil was performed considering the far nodes EE and WW as explicit. Then, these nodal values with their A EE and A WW coefficients are placed in the source term. Also, A EE and A WW are removed from the A P coefficient and added to the source term, resulting in the following modification of the previous A P and b P terms: In addition, to improve the numerical stability, the [ , ] m  factors in the case of the DD scheme are changed to [ ,1 ] The RNWF method was 4 times faster than DC and 1.3 times faster than NWF for the test cases studied in Chourushi (2018). However, for our case, the numerical implementation of the interface-capture schemes in the context of multiphase flows, the RNWF did not show the same excellent efficiency than in the cases studied by Chourushi and sometimes slightly altered the interface geometry.

Modified Normalized Weighting Factor (MNWF) method
Encouraged by the numerical stability and efficiency of the NWF method, we decided to review again its formulation to be used in the context of multiphase flows. Unlike the RNWF, the formulation of our new alternative called Modified Normalized Weighting Factor (MNWF) method starts from the initial formulation of the NWF method and not from its final discretized equation as the RNWF method does. First, the values of  and m for the interface-capture scheme are determined according to Eq. (14) and following the blending concept of Eq. (7). The values of  and m for some popular blended HR schemes used in the context of multiphase flows are listed in Table 1.
Second, to ensure numerical consistency, the m factor for the blended HR scheme is corrected with the CBC condition. For which (blended HR) f α  is explicitly calculated using the blended HR scheme, and then bounded as   Wacławczyk et al., 2007) Compressive scheme: The steady version of HYPER-C  1, 0 else ì ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï í ï ï ï ï ï ( Finally, for a better explanation of the MNWF method, we give the algebraic equation for the one-dimensional structured grid presented in Fig. 2  For all the cases, no under-relaxation factor is used, and the maximum number of iterations per time step is limited to 20, whereas the convergence criterion is different for each case. The solution algorithm employed for all the cases is a sequential one. At each time-step, the volume fraction equation (3) is firstly solved at the beginning of each time-step. Then, the new volume fraction field is used to compute the local density and viscosity using Eq. (4), and finally, the momentum and continuity equations (1) and (2) are solved by the predictor−corrector SIMPLE algorithm.
The first two test cases are the advection of the slotted circle in a rotational flow field introduced by Zalesak (1979) and the advection of a circle in a shear flow presented by Rudman (1997) which have simple exact solutions. These are frequently used in the multiphase community to check the performance of the advection schemes dealing with a non-uniform distribution of the Courant number and a considerable interface deformation, respectively (Ubbink and Issa, 1999;Zhang et al., 2014;Patel and Natarajan, 2015). For both, the accuracy of the simulation results is verified using the root mean square (RMS) error defined as where n i α is the numerical solution, a i α is the exact solution, and N is the total number of control volumes. The initialization of the volume fraction field for these cases are illustrated in Fig. 3.

Advection of a Slotted Circle in a rotational flow
A circle with a diameter of 1 m and a slot of width 0.12 m and depth 0.62 m is centered at (2, 2.65) m of a square 4 m × 4 m domain, and exposed to a clockwise circular velocity given by where (x 0 , y 0 ) = (2, 2) is the center of the rotation and ω = 0.5 rad/s is the constant angular velocity. The time needed for one rotation is 12.57 s. The problem is discretized with a structured grid of 200×200 square control volumes, and it is solved for five different time-step sizes which produce a maximum Courant number (Co) of 0.2, 0.4, 0.6, and 0.8 at point (2, 2.15) m (Darwish and Moukalled, 2006). The convergence criterion is ζ = 5×10 − 3 . The stacked columns shown in Fig. 4 represent the total number of iterations that each HR scheme performed to converge during one rotation of the slotted circle at different maximum Courant numbers. The bottom layer of the column represents the number of iterations performed using the DC method, the second relates to the DWF, the third to RNWF, and the last to MNWF. The results for the six different HR schemes are displayed in (a) to (f), with CICSAM shown in (a), MCICSAM-W in (b), MCICSAM-Z in (c), HRIC in (d), FBICS in (e), and CUIBS in (f). The efficiency between the four methods varies according to the HR scheme. However, for all of them, the new MNWF is the fastest to reach the convergence, whereas the DWF method is the slowest. For example, for CICSAM and MCICSAM-Z, the MNWF is about 20%−35% faster than DC and RNWF, and 35%−80% faster than DWF. Whereas for the HRIC and MCICSAM-W scheme, the MNWF works similarly to DC and RNWF but is noticeably better than DWF with about 20%−85% fewer total iterations. In the cases of FBICS and CUIBS, the total number of iterations used by MNWF is from 12%−25% less than by DC and RNWF, while 10%−85% less than for DWF.
On the one hand, for low Co, the influence of the implementation method is negligible for all schemes. On the other hand, for medium and high Co, the implementation method plays a significant role. For instance, the MNWF is about 10%−30% more efficient than DC and RNWF, and 50%−85% than DWF. Table 2 shows the accuracy of each HR scheme implemented through the four different techniques. For CICSAM, HRIC, and MCICSAM-W, the error varies according to the implementation technic used, and it is more significant for high Co. Nonetheless, the best performance remains in MNWF. To visualize this difference in accuracy, Fig. 5 shows the contour plots at Co = 0.8 of these three schemes. For the other three HR schemes, the error is independent of the implementation method.

Advection of a circle in a shear flow
The volume fraction field is initialized, as is shown in Fig.  3(b). A circle of 0.2π m diameter with its center at (0.5π, 0.2(1+π)) m filled with phase one is in a square domain of phase two. The two-phase configuration is exposed to a shear flow field described by where x, y ∈ [0, π]. The domain discretized with a uniform structured mesh consisting of 160×160 cells, and the time-step is chosen so that the local Courant number is 0.5. For observing the performance of the HR schemes together with the implementation method in the presence of interface deformation, the simulation is firstly run for n time-steps using the velocity defined in Eq. (31), then the flow is reversed, and the simulation is rerun for n time-steps. Hence, the interface should return to its initial shape. For this study, n = 1000 and n = 2000 are investigated. Figure 6 summarizes the total number of iterations required by the HR schemes implemented with the four methods DC, DWF, RNWF, and MNWF in the case of (a) n = 1000 and (b) n = 2000. The convergence criterium for these simulations is ζ = 5×10 − 3 . In case (a), 1000 forward steps followed by 1000 backward steps, the MNWF is 15%−50% faster than DC, 45%−84% faster than DWF, and about 3%−18% faster than RNWF. Only for MCICSAM-Z, the RNWF seems to be a better alternative. In case (b), 2000 steps forward followed by 2000 steps backwards, the high degree of interface deformation reduces the efficiency of the MNWF. Thus, the computational effort between MNWF and RNWF are similar, and the MNWF is now only 5%−27% faster than DC and 45%−58% faster than DWF. Nevertheless, MNWF is still a good option for the implementation of the analyzed schemes. Table 3 contains the errors of these simulations, which are slightly different for all schemes and implementation techniques. For n = 1000, the six advection schemes almost recover the initial shape. Therefore, the errors are small and fluctuate between 1.83×10 − 2 and 4.55×10 − 2 . While for n = 2000, all schemes suffer from numerical diffusion, independently of the implementation technique. Thus, the errors are more significant and range from 8.87×10 − 2 to 11.4×10 − 2 . Figure 7 depicts the results of two of the tested schemes: CUIBS (low error values) and CICSAM (high error values) for (a) n = 1000, and (b) n = 2000 forward and backward steps.

Rising bubble
The third test case, the rising bubble, is one of the multiphase flow benchmark cases of Hysing et al. (2009) which is an example of bubble dynamics with strong surface tension effects. Figure 8 illustrates the initial configuration where Ω 2 denotes the region that the bubble occupies and δV size of the control volume. The above quantities are used to validate the multiphase code and also the quality of the overall solution method. The error of the simulation is quantified using the following relative error norm: where q t is the temporal evolution of the quantity x c or y c , and q ref is the reference solution presented in Hysing et al. (2009).
The computation is performed for three uniform structured grids which consist of 40×80, 80×160, and 160×320 hexahedral control volumes. The time steps selected for each mesh are 0.02 s, 0.01 s, and 0.005 s respectably. This selection results in a maximum local Co of 0.45 for the three grids. The convergence criterion is ζ = 5×10 − 4 for the volume fraction field and ζ u = 10 − 6 for the velocity field. The curvature is calculated with the height-function method (Malik et al., 2007). Figure 9 shows the total number of iterations for the coarse, medium, and fine grid during three seconds of the calculation for the six HR schemes implemented using the DC, DWF, RNWF, and MNWF methods. For the three girds in this test case, the efficiency of the DC, DWF, RNWF and MNWF methods are similar for CICSAM and MCICSAM-Z. Although for the coarse and medium grid, the RNWF seems to have a small advantage. On the other hand, for the three grids and the other HR schemes, the percentage difference between the four methods is similar. The MNWF method is 15%−27% faster than DC, 6%−13% than RNWF, and 20%−50% than DWF. Although the most influential part of this simulation is the calculation of the curvature, the effects of the implementation techniques remained, and again the MNWF shows the best performance and DWF the worst.
Besides that, the relative error of the results computed for the fine grid, presented in Table 4, reveals that the calculations performed using the MNWF method are in the same range of precision than with the other methods. Figure 10(a) shows the final position of the bubble solved   Hysing et al. (2009) with the CUIBS scheme in combination with MNWF for the fine grid, and (b) the evolution of the centre of mass with time, and (c) the rise velocity of the bubble for the six blended HR schemes implemented with the MNWF method.

3D rising bubble
In order to measure the efficiency of the new approach in 3D scenarios, the 2D rising bubble test case was extended to 3D according to the configuration showed in Turek et al. (2019). The physical properties, gravity, and surface tension coefficient are the same as for the 2D case. Only now all walls are set as non-slip boundary conditions. The computational domain is a uniform structured grid of 80×80×160 hexahedral control volumes divided into 32 blocks. The time step is 0.005 s which produces a maximum Co ~0.55. The curvature is calculated with the standard second-order central difference scheme (CDS). The convergence criterion is ζ = 5×10 − 4 for the volume fraction and ζ u = 10 − 7 for the velocity. For this part, the results obtained with CICSAM and MCICSAM-Z are not presented because their low performance shown in the 2D case is almost the same in the 3D case, and these impede the precise observation of the other results. Figure 11 shows the new initial geometric configuration, the simulation results of the position of the bubble at three seconds, and the evolution of the rise velocity. The four schemes MCICSAM-W, HRIC, FBICS, CUIBS, implemented with the MNWF method agree with the reference result published in Turek et al. (2019).
Concerning the efficiency of the MNWF method, the performance shown for the 2D case is almost the same for the 3D case. The MNWF method is the fastest and the DWF method the slowest. Except for the FBICS scheme where the quickest is the DC. See

Dam break flow impacting a rigid structure
Finally, to validate our proposed MNWF method for more realistic applications, the classic dam breaking example is computed. It was experimentally studied by Koshizuka (1995) to describe the collapse of a water column impacting a rigid structure. Figure 13 represents the geometry and physical parameters of the problem. It consists of a box open to the atmosphere that contains a water column which collapses and hits a rigid obstacle. The high-density fluid is water, and the low-density fluid is air, both with standard physical properties. The two-phase flow is considered laminar, and the surface tension effects are neglected. The computation is run for two grids consisting of hexahedral control volumes (CV), a coarse grid of 2384 CV and a fine grid of 9536 CV. Each mesh is formed by five blocks divided according to the contour of the obstacle and following a hyperbolic distribution as is shown in Fig.  13(b). The simulations are computed for 0.9 s with a variable time-step defined for each grid to maintain a Co below one. For these calculations, the adaptative time step is chosen because the kinematic energy of the water increases with the time, which produces a constant increase of the Courant number. If a small constant time step is selected, the calculation is stable but is inefficient, while a large value leads to divergence when the velocity of the water is higher. The convergence criterion for the volume fraction field is ζ = 5×10 − 4 , and ζ u = 10 − 7 for the velocity and pressure. Figure 14 presents the total number of iterations required in each grid for the six HR schemes implemented using DC, DWF, RNWF, and MNWF. For the two grids, the difference in the performance of the methods is less evident than for the other cases. It is because the timing control results in low Courant numbers over a large part of the simulation period, thus decreasing the difference in efficiency that exists between the methods. As shown in the first case, the advection of a Slotted Circle in a rotational flow, for low Co, all four methods act similarly. However, for these grids, one observes that for the CICSAM and MCICSAM-Z schemes, the new MNWF method is considerably better than DC and DWF, and slightly better than RNWF. Although RNWF is the quickest for MCICSAM-Z in the coarse grid domain. MNWF computed results with 31%−44% fewer iterations than DC, 50% fewer than DWF, and 2%−25% fewer than RNWF.
Meanwhile, for MCICSAM-W, HRIC, FBICS, and CUIBS, the advantage in speed of MNWF in comparison with the other techniques is less evident but existent. The MNWF  is approximately 5%−7% faster than DC or RNWF, and 15%−32% than DWF. Furthermore, in this case, we noted that the numerical precision of the implementation method directly affects the calculation of the velocity. Thus, the velocity needs more interactions to converge, especially for DWF as this method on many occasions does not achieve the converge criterion.
Finally, in Fig. 14, a qualitative comparison between the simulation results obtained for the fine grid using the MCICSAM-Z scheme implemented with the proposed method MNWF and the experimental results from Koshizuka (1995) are presented to demonstrate the quality of the simulations.

Conclusions
We presented a new methodology, the MNWF method, to deal with the implementation of blended HR schemes in the context of multiphase flow. The method is based on the NWF methodology and produces a system of tridiagonal equations. The resulting coefficient matrix is always diagonally dominant, with A P coefficients greater than zero, giving a numerically stable method without the necessity of under Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.