We propose novel third-order less oscillatory and less diffusive compact stencil-based upwind schemes for the approximation of the continuity equation. The proposed schemes are based on the constrained interpolation profile-conservative semi-Lagrangian schemes. An important feature of the proposed schemes is that the interpolation functions are constructed using only variables within one upwind cell (a cell average and two boundary values). Furthermore, the proposed schemes have third-order accuracy and are also less oscillatory, less diffusive, and fully conservative. The proposed schemes are validated through various benchmark problems and comparisons with experiments of two droplets collision/separation and droplet splashing. The numerical results have shown that the proposed schemes have third-order accuracy for smooth solution, and capture discontinuities and smooth solutions simultaneously without numerical oscillations. The proposed schemes can capture the secondary vorticity of lid-driven cavity flow of Re = 7500 with a Cartesian grid of 64 × 64. The numerical results of two droplets collision/separation of We = 40 show that the proposed schemes can reproduce droplets collision/separation with quite coarse grids. These numerical results of droplet splashing have demonstrated that proposed schemes can reduce numerical diffusions well against existing schemes and robust.

In this paper, we propose novel third-order less oscillatory and less diffusive compact stencil-based upwind schemes for the approximation of the continuity equation
(1)
where ϕ is the scalar and u is the velocity. The proposed method is based on the constrained interpolation profile-conservative semi-Lagrangian (CIP-CSL) schemes.15,16,18,27,34–36

The CIP-CSL schemes are solvers of the conservation equation based on a multi-moment concept, which uses both cell average and boundary value as variables. In CIP-CSL schemes, the order of accuracy can be improved by using both cell average and boundary value without using wider stencils (e.g., the CIP-CSL schemes can achieve third-order accuracy only with variables within a cell). On the other hand, most of other schemes such as the essentially non-oscillatory schemes (ENO),9,22,23 weighted essentially non-oscillatory schemes (WENO),13,17 and the piecewise parabolic method (PPM)6 use only single-moment (a type of variables, for instance, cell average if finite volume methods and point value if finite difference methods). In these single-moment-based schemes, a wider stencil must be used to increase the order of accuracy. Thus, the CIP-CSL schemes are a different type of schemes from finite volume or finite different schemes.

The CIP-CSL schemes have been applied to various fluid problems1,30,31,33 such as droplets collision/separation,44 droplet impact on dry surface,46 droplet splashing on dry surface,41–44 interactions between particles and free surface flow,44 and oceanic flows.37 The fluid solvers30,31,33,44 based on the CIP-CSL schemes are highly efficient especially for incompressible flow simulations. For instance, the fluid solvers can simulate droplet splashing on a hydrophobic substrate at least qualitatively with a coarse grid (192 × 192 × 48) and the numerical simulation was completed within 2 hours on a standard desktop computer (Intel Core i7–3820 3.6 GHz, 8 GB memory).44 As drawbacks of the CIP-CSL schemes against single-moment-based schemes, the schemes require additional memory and calculation time for these additional degrees of freedom (DoF) because CIP-CSL schemes use both cell average and boundary value as variables. Although the drawback on additional memory cannot be resolved, these additional DoF surely improve fluid calculations (i.e., numerical resolution is increased by the additional DoF) and memory is not so expensive recently. If the calculation time is also proportionally increased due to the additional DoF, the use of the CIP-CSL schemes would not have many advantages against single-moment-based schemes because it is almost like just increased numerical resolution of single-moment-based schemes. However, if the CIP-CSL schemes are used within semi-implicit formulations (the pressure is solved implicitly and velocity explicitly) such as VSIAM330,31,33 or VSIAM3-TM,44 the computational time is not increased proportionally. Such semi-implicit formulations are mainly used for incompressible flow simulations. In VSIAM3 and VSIAM3-TM, these additional DoF are used only for velocity and not used for pressure. Then, although the calculation cost of velocity is surely increased, the cost of pressure calculation is identical to that of single-moment-based schemes. As it is well known, in incompressible flow simulations based on the semi-implicit formulations, the pressure calculation dominates the total calculation time (typically >80%). Therefore, additional calculation time by the additional DoF is small for the total calculation time. However, the accuracy of velocity is meaningfully improved due to the additional DoF. If the CIP-CSL schemes are used for full explicit formulations (e.g., compressible flow simulations), the calculation cost is proportionally increased as DoF and the advantages are not so clear compared to incompressible flow cases. Therefore, in this paper, we focus only on incompressible flows as applications of the proposed schemes.

There are several compact stencil-based CIP-CSL schemes such as CIP-CSL2 (CIP-CSL with second-order polynomial function)36 and CIP-CSLR (CIP-CSL with rational function).35 In CSL2, which is based on a second-order polynomial interpolation function, three moments within the upwind cell (i.e., a cell average and two boundary values) are used to construct the interpolation function (the details are explained in Subsection II A). CSL2 has third-order accuracy even though the interpolation function is constructed based on moments within one upwind cell. However, as a disadvantage, CSL2 is oscillatory. In CSLR, the same moments with CSL2 are used to construct the interpolation function and oscillations have successfully been minimized by using the rational function. However, CSLR cannot maintain third-order accuracy.

In this paper, we propose the CSL2T (CSL second-order polynomial and hyperbolic tangent functions) and CSL2R (CSL second-order polynomial and rational functions) schemes, which can maintain third-order accuracy of CSL2 but can also minimize oscillations like CSLR. In this paper, we focus on solvers, which can construct interpolation functions using only variables within one upwind cell because such solvers can avoid loss of accuracy near wall boundaries (the proposed schemes can maintain third-order accuracy even near wall boundaries) and practically important. In this paper, we also propose non-oscillatory selector (NOS). NOS is used to switch between an oscillatory high-order schemes such as CSL2 and a non-oscillatory scheme (e.g., CSLT or CSLR). There are some similar concepts with NOS such as limiters,10,21 boundary variation diminishing (BVD),7,24 and multi-dimensional optimal order detection (MOOD).5 As an advantage of NOS against limiters, NOS can maintain third-order accuracy. Although we also tried some limiters, these limiters did not work well for CSL schemes. This will be because CSL schemes use not only cell averages like finite volume methods but also boundary values. In this paper, we have employed a hyperbolic tangent function12,32,39 as a non-oscillatory solution like BVD schemes. However, the selection strategies of high-order interpolation function and the hyperbolic tangent function are totally different. Although NOS uses discrete maximum principle to select a suitable interpolation function like MOOD, there are many differences with MOOD. For instance, although MOOD is for high-order finite volume schemes, NOS is for mixed formulations, which uses both finite volume and finite difference methods. MOOD has not employed the hyperbolic tangent function as a candidate of interpolation functions.

The details of proposed CIP-CSL schemes and some existing CIP-CSL schemes are given in Sec. II. In Sec. III, numerical results of sine wave propagation test, Jiang-Shu test, Burgers' equation, Zalesak problem, lid-driven cavity flows, droplet collision and separation, and droplet splashing are given. The summary is given in Sec. IV.

In this section, we review CSL2 and CSLR schemes, and propose third-order less oscillatory and less diffusive compact stencil-based upwind schemes (CSL2T and CSL2R with NOS). In Secs. II A and II B, CSL236 and CSLR35 are explained. In Sec. II C, we propose the CSLT scheme based on the hyperbolic tangent function. In Sec. II D, we propose NOS and the CSL2T scheme based on CSL2 and CSLT schemes. In Sec. II E, we also propose the CSL2R scheme based on CSL2, CSLR, and NOS.

The CIP-CSL2 scheme36 would be the simplest CIP-CSL formulation. The CIP-CSL2 scheme uses three moments in the upwind cell (i.e., one cell average ϕ ¯ i and two boundary values ϕ i 1 / 2 and ϕ i + 1 / 2) to interpolate within the upwind cell (i.e., between x i 1 / 2 and x i + 1 / 2) as shown in Fig. 1. Then, we can interpolate between x i 1 / 2 and x i + 1 / 2 by a quadratic function Φ i CSL 2 ( x ),
(2)
FIG. 1.

Schematic figure of the CIP-CSL2 scheme. u i 1 / 2 < 0 is assumed. The moments, which are indicated by gray color ( ϕ i 1 / 2 , ϕ ¯ i, and ϕ i + 1 / 2), are used to construct the quadratic interpolation function Φ i CSL 2 ( x ).

FIG. 1.

Schematic figure of the CIP-CSL2 scheme. u i 1 / 2 < 0 is assumed. The moments, which are indicated by gray color ( ϕ i 1 / 2 , ϕ ¯ i, and ϕ i + 1 / 2), are used to construct the quadratic interpolation function Φ i CSL 2 ( x ).

Close modal
By using the following constraints,
(3)
(4)
(5)
these coefficients C 2 , i CSL 2 , C 1 , i CSL 2 , and C 0 , i CSL 2 can be calculated as follows:
(6)
(7)
(8)
By using the interpolation function Φ i CSL 2 ( x ), the boundary value ϕ i 1 / 2 is updated by the conservation equation of a differential form
(9)
A semi-Lagrangian approach is used for Eq. (9)
(10)
(11)
The cell average ϕ ¯ i is update by a finite volume formulation
(12)
where F i 1 / 2 CSL 2 is the flux at the cell boundary x i 1 / 2,
(13)
The CIP-CSLR scheme is briefly explained here. In the CIP-CSLR scheme,35 the following function Φ i CSLR ( x ),
(14)
is used to interpolate between x i 1 / 2 and x i + 1 / 2. The coefficients, C 2 , i CSLR , C 1 , i CSLR, and C 0 , i CSLR, are determined as follows:
(15)
(16)
(17)
by using the same constraints with these of CSL2 [(3)–(5)]. Here, ε is a small number to avoid zero division. We used ε = 10 15 for all results in this paper. The rest of the procedure is the same with CSL2.
We propose the CSLT (CSL with tangent hyperbolic function) scheme, which uses a hyperbolic tangent function as the interpolation function (Fig. 2). A similar scheme has been proposed in Ref. 16. In Ref. 16, only cell averages were used to construct the hyperbolic tangent function (like standard finite volume methods) and a wider stencil (three cells) was used. In the proposed scheme, not only cell averages but also boundary values are used (like CIP-CSL schemes) and a shorter stencil (one cell average and two boundary values) is used. The interpolation function is defined as
(18)
where ϕ min = min ( ϕ i 1 / 2 , ϕ i + 1 / 2 ) , ϕ max = max ( ϕ i 1 / 2 , ϕ i + 1 / 2 ), and γ = sgn ( ϕ i + 1 / 2 ϕ i + 1 / 2 ). x i ̃ is the jump location and can be calculated by using (5). β is the parameter to control jump thickness. The rest of the procedure is the same with CSL2. The CIP-CSLT scheme is a non-oscillatory scheme because the interpolation function (the hyperbolic tangent function) is a monotonic function with a minimum of ϕ min and a maximum of ϕ max.
FIG. 2.

Schematic figure of the CIP-CSLT scheme. u i 1 / 2 < 0 is assumed. The moments, which are indicated by gray color ( ϕ i 1 / 2 , ϕ ¯ i, and ϕ i + 1 / 2), are used to construct the step function Φ i CSLT ( x ).

FIG. 2.

Schematic figure of the CIP-CSLT scheme. u i 1 / 2 < 0 is assumed. The moments, which are indicated by gray color ( ϕ i 1 / 2 , ϕ ¯ i, and ϕ i + 1 / 2), are used to construct the step function Φ i CSLT ( x ).

Close modal
We propose non-oscillatory selector (NOS), which can prevent numerical oscillations of a high-order scheme such as CSL2 by combining with a non-oscillatory scheme such as CSLT and CSLR. The basic idea of the proposed selector is that first we solve the equation using a high-order scheme like CSL2. If oscillations are detected for some cell averages and/or boundary values, we recalculate these oscillatory cell averages and/or boundary values using a non-oscillatory scheme like CSLT. If numerical oscillations are not detected, we simply use the solution of the high-order scheme. Some schematic examples are shown in Fig. 3. For instance, if CSL2 is used as the high-order scheme, CSL2 possibly causes numerical oscillations as shown in Figs. 3(a) and 3(b) because CSL2 employs the quadratic function as the interpolation function. If the interpolation functions, which are indicated by the dotted ellipses in Figs. 3(a) and 3(b), are used to update ϕ i + 1 / 2 , ϕ i + 1 / 2 possibly causes some numerical oscillations. For such cases (when an oscillation was detected), we recalculate ϕ i + 1 / 2 using CSLT as shown in Figs. 3(a) and 3(b). CSLT has employed the hyperbolic tangent function (a monotonic function) as the interpolation function, so we can prevent the oscillation. If numerical oscillations are not detected like Fig. 3(c), we simply use the solution by CSL2. We also consider a successive gradient in this selector. The procedure is as follows (here we consider that CSL2 is the high-order scheme and CSLT is the non-oscillatory scheme):
FIG. 3.

Schematic figure of the non-oscillatory selector (NOS) for three different cases. (a) and (b) show typical cases when CSL2 is oscillatory and (c) a typical non-oscillatory case of CSL2.

FIG. 3.

Schematic figure of the non-oscillatory selector (NOS) for three different cases. (a) and (b) show typical cases when CSL2 is oscillatory and (c) a typical non-oscillatory case of CSL2.

Close modal
  • Step 1: Update all cell averages ( ϕ ¯ i CSL 2 , n + 1) and boundary values ( ϕ i 1 / 2 CSL 2 , n + 1) using CSL2.

  • Step 2: If the cell average ϕ ¯ i CSL 2 , n + 1 is oscillatory, the fluxes F i 1 / 2 CSL 2 and F i + 1 / 2 CSL 2 are replaced with F i 1 / 2 CSLT and F i + 1 / 2 CSLT, and ϕ ¯ i is re-updated using CSLT (neighbor cell averages associated with F i 1 / 2 CSLT or F i + 1 / 2 CSLT are also re-updated). Then, we also re-update the boundary values, ϕ i 1 / 2 and ϕ i + 1 / 2 using CSLT. We assume that ϕ ¯ i CSL 2 , n + 1 is oscillatory when ( ϕ ¯ i CSL 2 , n + 1 > ϕ ¯ i max , n ε) or ( ϕ ¯ i CSL 2 , n + 1 < ϕ ¯ i min , n + ε). Here, ϕ ¯ i max , n is the maximum value of the variables (at time step n), which are used to calculate ϕ ¯ i CSL 2 , n + 1, and ϕ ¯ i min , n is the minimum value of these variables.

  • Step 3: If the boundary value ϕ i 1 / 2 CSL 2 , n + 1 is oscillatory, ϕ i 1 / 2 is re-updated using CSLT. Then, F i 1 / 2 CSL 2 is replaced with F i 1 / 2 CSLT (if it was not replaced in step 2) and also re-updates cell averages, which are associated with F i 1 / 2 CSLT. We assume that ϕ i 1 / 2 CSL 2 , n + 1 is oscillatory when ( ϕ i 1 / 2 CSL 2 , n + 1 > ϕ i 1 / 2 max , n ε) or ( ϕ i 1 / 2 CSL 2 , n + 1 < ϕ i 1 / 2 min , n + ε). Here, ϕ i 1 / 2 max , n is the maximum value of the variables (at time step n), which are used to calculate ϕ i 1 / 2 CSL 2 , n + 1, and ϕ i 1 / 2 min , n is the minimum value of these variables.

  • Exception CSLT cannot handle preexisting cusps (when ϕ ¯ i > max ( ϕ i 1 / 2 , ϕ i + 1 / 2 ) ε or ϕ ¯ i < min ( ϕ i 1 / 2 , ϕ i + 1 / 2 ) + ε) because then x ̃ is going to be infinity. For preexisting cusps, alternative solvers are required. In this paper, CSL2 is used for cusps when the successive gradient (ri) is less than α; otherwise, the first-order upwind scheme is used. The successive gradient is defined as
    (19)

    In this paper, α = 4.01 is mainly used.

NOS is intuitive and easy to implement. We call the scheme, which combined CSL2 and CSLT as CSL2T (CSL with second-order polynomial and hyperbolic tangent functions).

By using NOS, we can minimize numerical oscillations of CSL2 with any less-oscillatory scheme. In this paper, we also proposed CSL2R (CSL with second-order polynomial and rational functions) scheme, which is based on CSL2 and CSLR. The procedure is almost same with CSL2T. We can simply replace CSLT with CSLR in the CSL2T scheme.

The conservation equation (1) is solved with the initial condition ϕ ( x , 0 ) = sin ( 2 π x ). The domain [ 0 , 1 ], u(x)=1 and periodic boundary conditions are used. Four different grid sizes (N = 40, 80, 160, and 320) are used with Δ t = 0.4 Δ x and Δ x = 1 / N. Errors are defined as follows:
(20)
(21)
Table I shows the numerical results. CSL2 has third-order accuracy. CSLR cannot maintain third-order accuracy. CSL2T and CSL2R have third-order accuracy when α 4.001. This means that NOS does not disturb the accuracy of the CSL2 scheme when α 4.001. If α 4.00 (for instance, α = 3.999) is used, the order of accuracy cannot be maintained in both CSL2T and CSL2R. In this paper, α = 4.01 is used for the rest of results.
TABLE I.

L1 and L errors in sine wave propagation at t = 1.

Method N L1 error L1 order L error L order
CSL2  40  9.79 × 10 5  ⋯  1.53 × 10 4  ⋯ 
80  1.23 × 10 5  2.99  1.93 × 10 5  2.99 
160  1.53 × 10 6  3.01  2.41 × 10 6  3.00 
320  1.92 × 10 7  2.99  3.01 × 10 7  3.00 
CSLR  40  2.29 × 10 3  ⋯  1.19 × 10 2  ⋯ 
80  5.53 × 10 4  2.05  4.75 × 10 3  1.33 
160  1.29 × 10 4  2.10  1.83 × 10 3  1.38 
320  2.75 × 10 5  2.23  6.91 × 10 4  1.41 
CSL2T (any β, α 4.001 40  9.79 × 10 5  ⋯  1.53 × 10 4  ⋯ 
80  1.23 × 10 5  2.99  1.93 × 10 5  2.99 
160  1.53 × 10 6  3.01  2.41 × 10 6  3.00 
320  1.92 × 10 7  2.99  3.01 × 10 7  3.00 
CSL2T ( β = 4.0, α = 3.999 40  8.50 × 10 4  ⋯  4.12 × 10 3  ⋯ 
80  1.75 × 10 4  2.29  1.61 × 10 3  1.36 
160  3.22 × 10 5  2.44  5.86 × 10 4  1.46 
320  5.58 × 10 6  2.53  2.00 × 10 4  1.55 
CSL2R ( α 4.001 40  9.79 × 10 5  ⋯  1.53 × 10 4  ⋯ 
80  1.23 × 10 5  2.99  1.93 × 10 5  2.99 
160  1.53 × 10 6  3.01  2.41 × 10 6  3.00 
320  1.92 × 10 7  2.99  3.01 × 10 7  3.00 
CSL2R ( α = 3.999 40  8.99 × 10 4  ⋯  4.39 × 10 3  ⋯ 
80  2.00 × 10 4  2.17  1.86 × 10 3  1.24 
160  3.68 × 10 5  2.44  6.73 × 10 4  1.46 
320  7.23 × 10 6  2.35  2.59 × 10 4  1.38 
Method N L1 error L1 order L error L order
CSL2  40  9.79 × 10 5  ⋯  1.53 × 10 4  ⋯ 
80  1.23 × 10 5  2.99  1.93 × 10 5  2.99 
160  1.53 × 10 6  3.01  2.41 × 10 6  3.00 
320  1.92 × 10 7  2.99  3.01 × 10 7  3.00 
CSLR  40  2.29 × 10 3  ⋯  1.19 × 10 2  ⋯ 
80  5.53 × 10 4  2.05  4.75 × 10 3  1.33 
160  1.29 × 10 4  2.10  1.83 × 10 3  1.38 
320  2.75 × 10 5  2.23  6.91 × 10 4  1.41 
CSL2T (any β, α 4.001 40  9.79 × 10 5  ⋯  1.53 × 10 4  ⋯ 
80  1.23 × 10 5  2.99  1.93 × 10 5  2.99 
160  1.53 × 10 6  3.01  2.41 × 10 6  3.00 
320  1.92 × 10 7  2.99  3.01 × 10 7  3.00 
CSL2T ( β = 4.0, α = 3.999 40  8.50 × 10 4  ⋯  4.12 × 10 3  ⋯ 
80  1.75 × 10 4  2.29  1.61 × 10 3  1.36 
160  3.22 × 10 5  2.44  5.86 × 10 4  1.46 
320  5.58 × 10 6  2.53  2.00 × 10 4  1.55 
CSL2R ( α 4.001 40  9.79 × 10 5  ⋯  1.53 × 10 4  ⋯ 
80  1.23 × 10 5  2.99  1.93 × 10 5  2.99 
160  1.53 × 10 6  3.01  2.41 × 10 6  3.00 
320  1.92 × 10 7  2.99  3.01 × 10 7  3.00 
CSL2R ( α = 3.999 40  8.99 × 10 4  ⋯  4.39 × 10 3  ⋯ 
80  2.00 × 10 4  2.17  1.86 × 10 3  1.24 
160  3.68 × 10 5  2.44  6.73 × 10 4  1.46 
320  7.23 × 10 6  2.35  2.59 × 10 4  1.38 
We validate proposed CSL schemes through Jiang-Shu problem13  u(x)=1, the domain [ 1 , 1 ], N =200, Δ t = 0.4 Δ x , Δ x = 2 / N, and periodic boundary conditions are used in this test. The initial condition is given as
(22)
where
(23)
(24)
where a =0.5, z = 0.7 , δ = 0.005 , α J S = 10, and β J S = log ( 2 ) / ( 36 δ 2 ).

Figure 4 shows the results at t = 2 (1 period). CSL2 is oscillatory around discontinuities as shown in Fig. 4(a). Although CSLR can prevent numerical oscillations, the result is diffusive as shown in Fig. 4(b). CSL2T (β = 4) and CSL2R can minimize numerical oscillations of CSL2 without increasing diffusions much [Figs. 4(c) and 4(d)].

FIG. 4.

Numerical results of Jiang-Shu test at t = 2 (1 period). (a)–(d) are numerical results by CSL2, CSLR, CSL2T (β = 4), and CSL2R, respectively. N = 200 and CFL = 0.4 are used.

FIG. 4.

Numerical results of Jiang-Shu test at t = 2 (1 period). (a)–(d) are numerical results by CSL2, CSLR, CSL2T (β = 4), and CSL2R, respectively. N = 200 and CFL = 0.4 are used.

Close modal

Figure 5 shows the results when the parameter β in CSLT was altered. If a smaller β [e.g., Fig. 5(a)] is used, the solution becomes slightly diffusive and errors are increased as shown in Table II. If a larger β [e.g., Figs. 5(b) and 5(c)] is used, the solution becomes less diffusive and errors are decreased as shown in Table II. If a large β (e.g., β = 10) is used, CSL2T can delete numerical oscillations almost without increasing diffusion as shown in Fig. 5(c). Although it is better to use larger β for the 1D test problem, there are some side effects of sharp sigmoid function in multidimensional fluid simulations (an example is shown in lid-driven cavity flows later). Therefore, in this paper, β = 4 is mainly used. CSL2T is superior to CSL2R when β 3 is used in this test.

FIG. 5.

Numerical results of Jiang-Shu test at t = 2 (1 period). (a)–(c) are numerical results by CSL2T with β = 3, β = 5, and β = 10, respectively. N = 200 and CFL = 0.4 are used.

FIG. 5.

Numerical results of Jiang-Shu test at t = 2 (1 period). (a)–(c) are numerical results by CSL2T with β = 3, β = 5, and β = 10, respectively. N = 200 and CFL = 0.4 are used.

Close modal
TABLE II.

L1 and L errors in the complex wave propagation test at t = 1 (1 period). N = 200 and CFL = 0.4 are used.

L1 error L error
CSL2  1.49 × 10 2  2.46 × 10 1 
CSLR  1.90 × 10 2  2.60 × 10 1 
CSL2T (β = 3)  1.27 × 10 2  2.33 × 10 1 
CSL2T (β = 4)  1.22 × 10 2  2.33 × 10 1 
CSL2T (β = 5)  1.19 × 10 2  2.32 × 10 1 
CSL2T (β = 10)  1.12 × 10 2  2.31 × 10 1 
CSL2R  1.28 × 10 2  2.33 × 10 1 
L1 error L error
CSL2  1.49 × 10 2  2.46 × 10 1 
CSLR  1.90 × 10 2  2.60 × 10 1 
CSL2T (β = 3)  1.27 × 10 2  2.33 × 10 1 
CSL2T (β = 4)  1.22 × 10 2  2.33 × 10 1 
CSL2T (β = 5)  1.19 × 10 2  2.32 × 10 1 
CSL2T (β = 10)  1.12 × 10 2  2.31 × 10 1 
CSL2R  1.28 × 10 2  2.33 × 10 1 
In this test, we solve the inviscid Burgers' equation in its conservative formulation
(25)
with the initial condition u ( x , 0 ) = 0.5 + 0.4 cos ( 2 π x ). N =100, CFL = 0.2, and periodic boundary conditions are used. Figure 6 shows the results at t = 1. The reference solution is created by using CSL3CW with N = 10 000. CSL2 has numerical oscillation around the discontinuity. CSLR, CSL2T, and CSL2R can capture the discontinuity well without numerical oscillation. CSL2T and CSL2R are slightly less diffusive than CSLR around the discontinuity.
FIG. 6.

Numerical results of Burgers' equation at t = 1. N = 100 and CFL = 0.2 are used.

FIG. 6.

Numerical results of Burgers' equation at t = 1. N = 100 and CFL = 0.2 are used.

Close modal
Zalesak's test problem47 in which a notched circle is rotated is widely used as a test of scalar advection schemes. The initial condition is given by
(26)
(27)
The one revolution is completed with 628 time steps.

Figure 7 shows numerical results (top view of 0.5-contour) by CSL2, CSLR, CSL2T (β = 4), and CSL2R after one revolution. Both results show good agreement with the exact solution. Figure 8 shows the perspective views of the numerical results by CSL2, CSLR, CSL2T (β = 4), and CSL2R. As shown in Fig. 8, the numerical result by CSL2 is oscillatory (a), CSLR can minimize numerical oscillations but the result is diffusive panel (b). CSL2T and CSL2R can minimize oscillations without increasing numerical diffusions much as shown in Figs. 8(c) and 8(d), respectively. Table III shows L1 errors after one revolution and four revolutions. The result indicates that CSL2T is superior to CSL2R in this test and the trend is unchanged after four revolutions.

FIG. 7.

Top views of numerical results of Zalesak problem (after one revolution). (a)–(d) are numerical results by CSL2, CSLR, CSL2T (β = 4), and CSL2R, respectively. The dotted and solid lines represent 0.5 contour lines of numerical results and the exact solution, respectively. A Cartesian grid of 100 × 100 was used.

FIG. 7.

Top views of numerical results of Zalesak problem (after one revolution). (a)–(d) are numerical results by CSL2, CSLR, CSL2T (β = 4), and CSL2R, respectively. The dotted and solid lines represent 0.5 contour lines of numerical results and the exact solution, respectively. A Cartesian grid of 100 × 100 was used.

Close modal
FIG. 8.

Numerical results of Zalesak problem (after one revolution). (a)–(d) are numerical results by CSL2, CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 100 × 100 was used.

FIG. 8.

Numerical results of Zalesak problem (after one revolution). (a)–(d) are numerical results by CSL2, CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 100 × 100 was used.

Close modal
TABLE III.

L1 errors in Zalesak problem (after one revolution and four revolutions). A Cartesian grid of 100 × 100 was used.

Schemes L1 Error after one revolution L1 error after four revolutions
CSL2  1.55 × 10 2  2.21 × 10 2 
CSLR  1.73 × 10 2  2.75 × 10 2 
CSL2T (β = 4)  1.28 × 10 2  1.97 × 10 2 
CSL2R  1.30 × 10 2  1.99 × 10 2 
Schemes L1 Error after one revolution L1 error after four revolutions
CSL2  1.55 × 10 2  2.21 × 10 2 
CSLR  1.73 × 10 2  2.75 × 10 2 
CSL2T (β = 4)  1.28 × 10 2  1.97 × 10 2 
CSL2R  1.30 × 10 2  1.99 × 10 2 

The proposed CIP-CSL schemes were applied to lid-driven cavity flow problems, and these numerical results are compared with numerical results by Ghia et al.8  Figure 9 shows numerical results of lid-driven cavity flows of Re = 7500 by CSLR, CSL2T (β = 4), and CSL2R (CSL2 is not considered hereafter because such oscillatory schemes cause some issues such as unnatural oscillations in solutions, and instability of calculations). A Cartesian grid of 128 × 128 was used. All numerical results have shown good agreements with the reference solution by Ghia (a grid of 256 × 256 was used in the reference solution). However, as shown in Fig. 9, we cannot observe differences of these results clearly with the resolution 128 × 128. Therefore, we compare these results with a lower numerical resolution. Figure 10 shows the numerical results when a Cartesian grid of 64 × 64 was used. In this resolution, it is clear that CSL2T and CSL2R are superior to CSLR. To observe the differences among CSL2T, CSL2R, and CSLR more clearly, we plotted numerical results of CSLR, CSL2T, and CSL2R together, and magnified the area 0.5 x 0 and 1 y 0.5. As shown in Fig. 11, CSL2T is slightly superior to CSL2R. We also study the effect of β of CSL2T in this test. As shown in Fig. 12(a), if β = 3 is used, CSL2T is even less accurate than CSLR. If β is increased from β = 3 up to around β = 4.5, the results are improved and are superior to results by CSL2R as well as CSLR. However, as shown in Fig. 12(b), if β is increased further around from β = 4.5, the results are getting worse. Based on these results in Fig. 12, optimal β would be around 4 β 5. If around 4 β 5 is used, CSL2T is superior to CSL2R. In this paper, β = 4 is used hereafter.

FIG. 9.

Numerical results of lid-driven cavity flows of Re = 7500. (a)–(c) are numerical results by CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 128 × 128 was used. The solid and dot lines represent x- and y-components of the velocity fields on the lines x =0 and y =0, respectively. Dots represent numerical results by Ghia et al.8 

FIG. 9.

Numerical results of lid-driven cavity flows of Re = 7500. (a)–(c) are numerical results by CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 128 × 128 was used. The solid and dot lines represent x- and y-components of the velocity fields on the lines x =0 and y =0, respectively. Dots represent numerical results by Ghia et al.8 

Close modal
FIG. 10.

Numerical results of lid-driven cavity flows of Re = 7500. (a)–(c) are numerical results by CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 64 × 64 was used.

FIG. 10.

Numerical results of lid-driven cavity flows of Re = 7500. (a)–(c) are numerical results by CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 64 × 64 was used.

Close modal
FIG. 11.

Numerical results of lid-driven cavity flows of Re = 7500. These lines are numerical results by CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 64 × 64 was used.

FIG. 11.

Numerical results of lid-driven cavity flows of Re = 7500. These lines are numerical results by CSLR, CSL2T (β = 4), and CSL2R, respectively. A Cartesian grid of 64 × 64 was used.

Close modal
FIG. 12.

Numerical results of lid-driven cavity flows of Re = 7500. (a) Represents numerical results by CSL2T with β = 3, β = 3.5, β = 4, and β = 4.5. (b) Represents numerical results by CSL2T with β = 4, β = 5, β = 6, and β = 10. A Cartesian grid of 64 × 64 was used.

FIG. 12.

Numerical results of lid-driven cavity flows of Re = 7500. (a) Represents numerical results by CSL2T with β = 3, β = 3.5, β = 4, and β = 4.5. (b) Represents numerical results by CSL2T with β = 4, β = 5, β = 6, and β = 10. A Cartesian grid of 64 × 64 was used.

Close modal

Figures 13(a) and 13(b) show numerical results of Re = 1000 and 5000 by CSL2T (β = 4), respectively. CSL2T can also simulate cavity flows of other Reynolds numbers well. We have also confirmed that CSLR and CSL2R schemes can also simulate cavity flows of these Reynolds number well. Figure 14 shows the velocity field of Re = 7500 by CSL2T when a Cartesian grid of 64 × 64 was used. Secondary vorticity at bottom right corner has been captured with the resolution of 64 × 64.

FIG. 13.

Numerical results of lid-driven cavity flows of Re = 1000 (a) and 5000 (b) when a Cartesian grid of 128 × 128 was used.

FIG. 13.

Numerical results of lid-driven cavity flows of Re = 1000 (a) and 5000 (b) when a Cartesian grid of 128 × 128 was used.

Close modal
FIG. 14.

Numerical results of lid-driven cavity flows of Re = 7500 by CSL2T (β = 4). A Cartesian grid of 64 × 64 was used.

FIG. 14.

Numerical results of lid-driven cavity flows of Re = 7500 by CSL2T (β = 4). A Cartesian grid of 64 × 64 was used.

Close modal

We performed three-dimensional numerical simulations of two droplets collision and separation, which involves topology change of liquid interfaces.3 The numerical method to simulate free surface flows is based on VSIAM3-TM,44 the CLSVOF (coupled level set19,26 and volume-of-fluid11,14) method,25,38 the THINC/WLIC (tangent of hyperbola for interface capturing/weighted line interface calculation) scheme,12,32,39 and the density scaled balanced continuum surface force model43,45 with level set curvature correction.42 For the full details of the implementation of the free surface flow solver, see Refs. 40 and 42–44.

Figures 15(a)–15(c) show snapshots of the numerical results of two droplets collision and separation (We = 40) by CSLR, CSL2R, and CSL2R, respectively. The corresponding experimental result can be found in Ref. 3. In these numerical simulations, quantitative parameters were used. The density ratio is 1.25:1000 (air:liquid). The mesh size Δ = D / 7.4 was used, where D is the diameter of initial droplets. In this numerical resolution ( Δ = D / 7.4), only CSL2T could reproduce the separation, and CSL2R and CSLR failed to do it. However, if the resolution was increased, CSLR and CSL2R could also reproduce the separation with Δ = D / 7.5 and Δ = D / 8.0, respectively. The numerical results have demonstrated that CSL2T and CSL2R can reduce numerical diffusion in the free surface flow application compared to CSLR, and CSL2T is slightly better than CSL2R in this test. Although there are some previous numerical work of droplet collision and separations,2,4,20,28 compared to these previous work, our simulation has captured the phenomena with a quite course grid.

FIG. 15.

Numerical results of two droplet collisions and separation by CSLR (a), CSL2T (b), and CSL2R (c) with the corresponding experimental result (d).3 The time evolution is from right to left. The mesh size is Δ = D / 7.4. Reproduced with permission from N. Ashgriz and J. Y. Poo, “Coalescence and separation in binary collisions of liquid drops,” J. Fluid Mech. 221, 183 (1990). Copyright 1990 Cambridge University Press.

FIG. 15.

Numerical results of two droplet collisions and separation by CSLR (a), CSL2T (b), and CSL2R (c) with the corresponding experimental result (d).3 The time evolution is from right to left. The mesh size is Δ = D / 7.4. Reproduced with permission from N. Ashgriz and J. Y. Poo, “Coalescence and separation in binary collisions of liquid drops,” J. Fluid Mech. 221, 183 (1990). Copyright 1990 Cambridge University Press.

Close modal

We performed numerical simulations of prompt splashing using CSLR, CSL2T, and CSL2R with VSIAM3-TM.44 The corresponding experimental result can be found in Ref. 29. In the numerical simulations, quantitative parameters, the densities ρ liquid = 1000, ρ air = 1.25 kg/m3, viscosities μ liquid = 1.0 × 10 3, μ air = 1.82 × 10 5 Pa s, surface tension σ = 7.2 × 10 2 N/m, gravity 9.8 m/s2, initial droplet diameter D =1.86 mm, impact speed 2.98 m/s, and the equilibrium contact angle 163° are used. A Cartesian grid of 224 × 224 × 48 is used. The mesh size Δ = D / 45 was used. Figure 16 shows the result. CSLR, CSL2T, and CSL2R can capture the physics of droplet splashing including satellite droplets and spikes. It seems that CSL2T and CSL2R capture satellite droplets and spikes more than CSLR and that the results by CSLR are more diffusive than these by CSL2T and CSL2R. The result has also shown the CSL2T and CSL2R as well as CSLR are robust.

FIG. 16.

A comparison among the numerical results by CSLR (a), CSL2T (b), and CSL2R (c). A distilled water droplet of 1.86 (mm) impacts onto a superhydrophobic substrate (the equilibrium angle is 163°). The droplet impact speed is 2.980 (m/s). A Cartesian grid of 224 × 224 × 48 and α = 1.5 Δ x is used. The mesh size is Δ = D / 45. The corresponding experimental result can be found in Ref. 29.

FIG. 16.

A comparison among the numerical results by CSLR (a), CSL2T (b), and CSL2R (c). A distilled water droplet of 1.86 (mm) impacts onto a superhydrophobic substrate (the equilibrium angle is 163°). The droplet impact speed is 2.980 (m/s). A Cartesian grid of 224 × 224 × 48 and α = 1.5 Δ x is used. The mesh size is Δ = D / 45. The corresponding experimental result can be found in Ref. 29.

Close modal

We proposed novel third-order less oscillatory and less diffusive compact stencil-based upwind schemes, which are called the CIP-CSL2T and CIP-CSL2R schemes for the approximation of hyperbolic conservation laws. We also proposed NOS, which can select a suitable CSL scheme between CSL2 and CSLT in CSL2T, and CSL2 and CSLR in CSL2R. Both CSL2T and CSL2R have third-order accuracy for the sine wave test and capture discontinuities and smooth solutions without numerical oscillations. The proposed schemes could capture the secondary vorticity at bottom right corner of lid-driven cavity flow of Re = 7500 with a Cartesian grid of 64 × 64. The numerical results of two droplets collision/separation of We = 40 show that the proposed schemes can reproduce droplets collision/separation with low numerical resolutions of 7.4–7.5 meshes for the diameter of initial droplets. The numerical results of droplet splashing have demonstrated that proposed schemes can reduce numerical diffusions well against CSLR and robust.

This research was supported in part by FLEXIS and IROHMS, which are part-funded by the European Regional Development Fund (ERDF), through the Welsh Government. The numerical simulations were partially conducted on computers at Yukawa Institute of Theoretical Physics in Kyoto University and at HPC Wales. We acknowledge the support of the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via Welsh Government.

The author declares that there are no conflicts to disclose.

Kensuke Yokoi: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.

1.
M.
Al-Mosallam
and
K.
Yokoi
, “
Efficient implementation of volume/surface integrated average-based multi-moment method
,”
Int. J. Comput. Methods
14
(
02
),
1750010
(
2017
).
2.
A.
Amani
,
N.
Balcázar
,
E.
Gutiérrez
, and
A.
Oliva
, “
Numerical study of binary droplets collision in the main collision regimes
,”
Chem. Eng. J.
370
,
477
498
(
2019
).
3.
N.
Ashgriz
and
J. Y.
Poo
, “
Coalescence and separation in binary collisions of liquid drops
,”
J. Fluid Mech.
221
,
183
204
(
1990
).
4.
X.
Chen
and
V.
Yang
, “
Direct numerical simulation of multiscale flow physics of binary droplet collision
,”
Phys. Fluids
32
(
6
),
062103
(
2020
).
5.
S.
Clain
,
S.
Diot
, and
R.
Loubère
, “
A high-order finite volume method for systems of conservation laws–multi-dimensional optimal order detection mood
,”
J. Comput. Phys.
230
(
10
),
4028
4050
(
2011
).
6.
P.
Colella
and
P. R.
Woodward
, “
The piecewise parabolic method (PPM) for gas-dynamical simulations
,”
J. Comput. Phys.
54
(
1
),
174
201
(
1984
).
7.
X.
Deng
,
Y.
Shimizu
,
B.
Xie
, and
F.
Xiao
, “
Constructing higher order discontinuity-capturing schemes with upwind-biased interpolations and boundary variation diminishing algorithm
,”
Comput. Fluids
200
,
104433
(
2020
).
8.
U.
Ghia
,
K.
Ghia
, and
C.
Shin
, “
High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method
,”
J. Comput. Phys.
48
(
3
),
387
411
(
1982
).
9.
A.
Harten
,
B.
Engquist
,
S.
Osher
, and
S. R.
Chakravarthy
, “
Uniformly high order accurate essentially non-oscillatory schemes, III
,”
J. Comput. Phys.
71
(
2
),
231
303
(
1987
).
10.
A.
Harten
and
S.
Osher
, “
Uniformly high-order accurate nonoscillatory schemes. I
,”
SIAM J. Numer. Anal.
24
(
2
),
279
309
(
1987
).
11.
C.
Hirt
and
B.
Nichols
, “
Volume of fluid (VOF) method for the dynamics of free boundaries
,”
J. Comput. Phys.
39
(
1
),
201
225
(
1981
).
12.
S.
Ii
,
K.
Sugiyama
,
S.
Takeuchi
,
S.
Takagi
,
Y.
Matsumoto
, and
F.
Xiao
, “
An interface capturing method with a continuous function: The THINC method with multi-dimensional reconstruction
,”
J. Comput. Phys.
231
(
5
),
2328
2358
(
2012
).
13.
G.-S.
Jiang
and
C.-W.
Shu
, “
Efficient implementation of weighted ENO schemes
,”
J. Comput. Phys.
126
(
1
),
202
228
(
1996
).
14.
J.
Li
, “
Calcul d'interface affine par morceaux
,”
C. R. Acad. Sci., Ser. II: Méc., Phys., Chim., Astron.
320
(
8
),
391
396
(
1995
).
15.
Q.
Li
,
S.
Omar
,
X.
Deng
, and
K.
Yokoi
, “
Constrained interpolation profile conservative semi-Lagrangian scheme based on third-order polynomial functions and essentially non-oscillatory (CIP-CSL3ENO) scheme
,”
Commun. Comput. Phys.
22
(
3
),
765
788
(
2017
).
16.
Q.
Li
,
J.
Xia
,
K.
Yokoi
, and
S.
Omar
, “
Boundary variation diminished conservative semi-Lagrangian method for both compressible and incompressible flows
,”
Phys. Fluids
33
(
11
),
117114
(
2021
).
17.
X.-D.
Liu
,
S.
Osher
, and
T.
Chan
, “
Weighted essentially non-oscillatory schemes
,”
J. Comput. Phys.
115
(
1
),
200
212
(
1994
).
18.
N.
Onodera
,
T.
Aoki
, and
K.
Yokoi
, “
A fully conservative high-order upwind multi-moment method using moments in both upwind and downwind cells
,”
Int. J. Numer. Methods Fluids
82
(
8
),
493
511
(
2016
).
19.
S.
Osher
and
J. A.
Sethian
, “
Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations
,”
J. Comput. Phys.
79
(
1
),
12
49
(
1988
).
20.
Y.
Pan
and
K.
Suga
, “
Numerical simulation of binary liquid droplet collision
,”
Phys. Fluids
17
(
8
),
082105
(
2005
).
21.
P. L.
Roe
, “
Characteristic-based schemes for the Euler equations
,”
Annu. Rev. Fluid Mech.
18
(
1
),
337
365
(
1986
).
22.
C.-W.
Shu
and
S.
Osher
, “
Efficient implementation of essentially non-oscillatory shock-capturing schemes
,”
J. Comput. Phys.
77
(
2
),
439
471
(
1988
).
23.
C.-W.
Shu
and
S.
Osher
, “
Efficient implementation of essentially non-oscillatory shock-capturing schemes, II
,”
J. Comput. Phys.
83
(
1
),
32
78
(
1989
).
24.
Z.
Sun
,
S.
Inaba
, and
F.
Xiao
, “
Boundary variation diminishing (BVD) reconstruction: A new approach to improve Godunov schemes
,”
J. Comput. Phys.
322
,
309
325
(
2016
).
25.
M.
Sussman
and
E. G.
Puckett
, “
A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows
,”
J. Comput. Phys.
162
(
2
),
301
337
(
2000
).
26.
M.
Sussman
,
P.
Smereka
, and
S.
Osher
, “
A level set approach for computing solutions to incompressible two-phase flow
,”
J. Comput. Phys.
114
(
1
),
146
159
(
1994
).
27.
R.
Tanaka
,
T.
Nakamura
, and
T.
Yabe
, “
Constructing exactly conservative scheme in a non-conservative form
,”
Comput. Phys. Commun.
126
(
3
),
232
243
(
2000
).
28.
S.
Tanguy
and
A.
Berlemont
, “
Application of a level set method for simulation of droplet collisions
,”
Int. J. Multiphase Flow
31
(
9
),
1015
1035
(
2005
).
29.
P.
Tsai
,
S.
Pacheco
,
C.
Pirat
,
L.
Lefferts
, and
D.
Lohse
, “
Drop impact upon micro- and nanostructured superhydrophobic surfaces
,”
Langmuir
25
(
20
),
12293
12298
(
2009
).
30.
F.
Xiao
, “
Unified formulation for compressible and incompressible flows by using multi-integrated moments. I. One-dimensional inviscid compressible flow
,”
J. Comput. Phys.
195
(
2
),
629
654
(
2004
).
31.
F.
Xiao
,
R.
Akoh
, and
S.
Ii
, “
Unified formulation for compressible and incompressible flows by using multi-integrated moments. II. Multi-dimensional version for compressible and incompressible flows
,”
J. Comput. Phys.
213
(
1
),
31
56
(
2006
).
32.
F.
Xiao
,
Y.
Honma
, and
T.
Kono
, “
A simple algebraic interface capturing scheme using hyperbolic tangent function
,”
Int. J. Numer. Methods Fluids
48
(
9
),
1023
1040
(
2005
).
33.
F.
Xiao
,
A.
Ikebata
, and
T.
Hasegawa
, “
Numerical simulations of free-interface fluids by a multi-integrated moment method
,”
Comput. Struct.
83
(
6–7
),
409
423
(
2005
).
34.
F.
Xiao
and
T.
Yabe
, “
Completely conservative and oscillationless semi-Lagrangian schemes for advection transportation
,”
J. Comput. Phys.
170
(
2
),
498
522
(
2001
).
35.
F.
Xiao
,
T.
Yabe
,
X.
Peng
, and
H.
Kobayashi
, “
Conservative and oscillation-less atmospheric transport schemes based on rational functions
,”
J. Geophys. Res.: Atmos.
107
(
D22
),
ACL 2–1
, https://doi.org/10.1029/2001JD001532 (
2002
).
36.
T.
Yabe
,
R.
Tanaka
,
T.
Nakamura
, and
F.
Xiao
, “
An exactly conservative semi-Lagrangian scheme (CIP-CSL) in one dimension
,”
Mon. Weather Rev.
129
(
2
),
332
344
(
2001
).
37.
S.
Yamashita
,
C.
Chen
,
K.
Takahashi
, and
F.
Xiao
, “
Large scale numerical simulations for multi-phase fluid dynamics with moving interfaces
,”
Int. J. Comput. Fluid Dyn.
22
(
6
),
405
410
(
2008
).
38.
K.
Yokoi
, “
Numerical method for complex moving boundary problems in a Cartesian fixed grid
,”
Phys. Rev. E
65
,
055701
(
May 2002
).
39.
K.
Yokoi
, “
Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm
,”
J. Comput. Phys.
226
(
2
),
1985
2002
(
2007
).
40.
K.
Yokoi
, “
A numerical method for free-surface flows and its application to droplet impact on a thin liquid layer
,”
J. Sci. Comput.
35
(
2
),
372
396
(
2008
).
41.
K.
Yokoi
, “
Numerical studies of droplet splashing on a dry surface: Triggering a splash with the dynamic contact angle
,”
Soft Matter
7
,
5120
5123
(
2011
).
42.
K.
Yokoi
, “
A practical numerical framework for free surface flows based on CLSVOF method, multi-moment methods and density-scaled CSF model: Numerical simulations of droplet splashing
,”
J. Comput. Phys.
232
(
1
),
252
271
(
2013
).
43.
K.
Yokoi
, “
A density-scaled continuum surface force model within a balanced force formulation
,”
J. Comput. Phys.
278
,
221
228
(
2014
).
44.
K.
Yokoi
,
M.
Furuichi
, and
M.
Sakai
, “
An efficient multi-dimensional implementation of VSIAM3 and its applications to free surface flows
,”
Phys. Fluids
29
(
12
),
121611
(
2017
).
45.
K.
Yokoi
,
R.
Onishi
,
X.-L.
Deng
, and
M.
Sussman
, “
Density-scaled balanced continuum surface force model with a level set based curvature interpolation technique
,”
Int. J. Comput. Methods
13
(
04
),
1641004
(
2016
).
46.
K.
Yokoi
,
D.
Vadillo
,
J.
Hinch
, and
I.
Hutchings
, “
Numerical studies of the influence of the dynamic contact angle on a droplet impacting on a dry surface
,”
Phys. Fluids
21
(
7
),
072102
(
2009
).
47.
S. T.
Zalesak
, “
Fully multidimensional flux-corrected transport algorithms for fluids
,”
J. Comput. Phys.
31
(
3
),
335
362
(
1979
).
Published open access through an agreement with Cardiff University School of Engineering