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.

## I. INTRODUCTION

*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 problems^{1,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 solvers^{30,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 VSIAM3^{30,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 function^{12,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.

## II. NUMERICAL METHODS

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, CSL2^{36} and CSLR^{35} 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.

### A. CIP-CSL2 scheme

^{36}would be the simplest CIP-CSL formulation. The CIP-CSL2 scheme uses three moments in the upwind cell (i.e., one cell average $ \varphi \xaf i$ and two boundary values $ \varphi i \u2212 1 / 2$ and $ \varphi i + 1 / 2$) to interpolate within the upwind cell (i.e., between $ x i \u2212 1 / 2$ and $ x i + 1 / 2$) as shown in Fig. 1. Then, we can interpolate between $ x i \u2212 1 / 2$ and $ x i + 1 / 2$ by a quadratic function $ \Phi i CSL 2 ( x )$,

### B. CIP-CSLR scheme

^{35}the following function $ \Phi i CSLR ( x )$,

*ε*is a small number to avoid zero division. We used $ \epsilon = 10 \u2212 15$ for all results in this paper. The rest of the procedure is the same with CSL2.

### C. CIP-CSLT scheme

*β*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 $ \varphi min$ and a maximum of $ \varphi max$.

### D. Non-oscillatory selector (NOS) and CIP-CSL2T scheme

**Step 1:**Update all cell averages ( $ \varphi \xaf i CSL 2 , n + 1$) and boundary values ( $ \varphi i \u2212 1 / 2 CSL 2 , n + 1$) using CSL2.**Step 2:**If the cell average $ \varphi \xaf i CSL 2 , n + 1$ is oscillatory, the fluxes $ F i \u2212 1 / 2 CSL 2$ and $ F i + 1 / 2 CSL 2$ are replaced with $ F i \u2212 1 / 2 CSLT$ and $ F i + 1 / 2 CSLT$, and $ \varphi \xaf i$ is re-updated using CSLT (neighbor cell averages associated with $ F i \u2212 1 / 2 CSLT$ or $ F i + 1 / 2 CSLT$ are also re-updated). Then, we also re-update the boundary values, $ \varphi i \u2212 1 / 2$ and $ \varphi i + 1 / 2$ using CSLT. We assume that $ \varphi \xaf i CSL 2 , n + 1$ is oscillatory when ( $ \varphi \xaf i CSL 2 , n + 1 > \varphi \xaf i max , n \u2212 \epsilon $) or ( $ \varphi \xaf i CSL 2 , n + 1 < \varphi \xaf i min , n + \epsilon $). Here, $ \varphi \xaf i max , n$ is the maximum value of the variables (at time step n), which are used to calculate $ \varphi \xaf i CSL 2 , n + 1$, and $ \varphi \xaf i min , n$ is the minimum value of these variables.**Step 3:**If the boundary value $ \varphi i \u2212 1 / 2 CSL 2 , n + 1$ is oscillatory, $ \varphi i \u2212 1 / 2$ is re-updated using CSLT. Then, $ F i \u2212 1 / 2 CSL 2$ is replaced with $ F i \u2212 1 / 2 CSLT$ (if it was not replaced in step 2) and also re-updates cell averages, which are associated with $ F i \u2212 1 / 2 CSLT$. We assume that $ \varphi i \u2212 1 / 2 CSL 2 , n + 1$ is oscillatory when ( $ \varphi i \u2212 1 / 2 CSL 2 , n + 1 > \varphi i \u2212 1 / 2 max , n \u2212 \epsilon $) or ( $ \varphi i \u2212 1 / 2 CSL 2 , n + 1 < \varphi i \u2212 1 / 2 min , n + \epsilon $). Here, $ \varphi i \u2212 1 / 2 max , n$ is the maximum value of the variables (at time step n), which are used to calculate $ \varphi i \u2212 1 / 2 CSL 2 , n + 1$, and $ \varphi i \u2212 1 / 2 min , n$ is the minimum value of these variables.**Exception**CSLT cannot handle preexisting cusps (when $ \varphi \xaf i > max ( \varphi i \u2212 1 / 2 , \varphi i + 1 / 2 ) \u2212 \epsilon $ or $ \varphi \xaf i < min ( \varphi i \u2212 1 / 2 , \varphi i + 1 / 2 ) + \epsilon $) because then $ x \u0303$ is going to be infinity. For preexisting cusps, alternative solvers are required. In this paper, CSL2 is used for cusps when the successive gradient (*r*) is less than_{i}*α*; otherwise, the first-order upwind scheme is used. The successive gradient is defined as$ r i = max ( | \varphi \xaf i \u2212 \varphi i \u2212 1 / 2 | | \varphi i + 1 / 2 \u2212 \varphi i \xaf | , | \varphi i + 1 / 2 \u2212 \varphi \xaf i | | \varphi \xaf i \u2212 \varphi i \u2212 1 / 2 | ) .$In this paper, $ \alpha = 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).

### E. CIP-CSL2R scheme

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.

## III. NUMERICAL RESULTS

### A. Sine wave propagation

*u*(

*x*)

*=*

*1 and periodic boundary conditions are used. Four different grid sizes (N = 40, 80, 160, and 320) are used with $ \Delta t = 0.4 \Delta x$ and $ \Delta x = 1 / N$. Errors are defined as follows:*

Method . | N . | L_{1} error
. | L_{1} order
. | $ L \u221e$ error . | $ L \u221e$ order . |
---|---|---|---|---|---|

CSL2 | 40 | $ 9.79 \xd7 10 \u2212 5$ | ⋯ | $ 1.53 \xd7 10 \u2212 4$ | ⋯ |

80 | $ 1.23 \xd7 10 \u2212 5$ | 2.99 | $ 1.93 \xd7 10 \u2212 5$ | 2.99 | |

160 | $ 1.53 \xd7 10 \u2212 6$ | 3.01 | $ 2.41 \xd7 10 \u2212 6$ | 3.00 | |

320 | $ 1.92 \xd7 10 \u2212 7$ | 2.99 | $ 3.01 \xd7 10 \u2212 7$ | 3.00 | |

CSLR | 40 | $ 2.29 \xd7 10 \u2212 3$ | ⋯ | $ 1.19 \xd7 10 \u2212 2$ | ⋯ |

80 | $ 5.53 \xd7 10 \u2212 4$ | 2.05 | $ 4.75 \xd7 10 \u2212 3$ | 1.33 | |

160 | $ 1.29 \xd7 10 \u2212 4$ | 2.10 | $ 1.83 \xd7 10 \u2212 3$ | 1.38 | |

320 | $ 2.75 \xd7 10 \u2212 5$ | 2.23 | $ 6.91 \xd7 10 \u2212 4$ | 1.41 | |

CSL2T (any β, $ \alpha \u2273 4.001$) | 40 | $ 9.79 \xd7 10 \u2212 5$ | ⋯ | $ 1.53 \xd7 10 \u2212 4$ | ⋯ |

80 | $ 1.23 \xd7 10 \u2212 5$ | 2.99 | $ 1.93 \xd7 10 \u2212 5$ | 2.99 | |

160 | $ 1.53 \xd7 10 \u2212 6$ | 3.01 | $ 2.41 \xd7 10 \u2212 6$ | 3.00 | |

320 | $ 1.92 \xd7 10 \u2212 7$ | 2.99 | $ 3.01 \xd7 10 \u2212 7$ | 3.00 | |

CSL2T ( $ \beta = 4.0$, $ \alpha = 3.999$) | 40 | $ 8.50 \xd7 10 \u2212 4$ | ⋯ | $ 4.12 \xd7 10 \u2212 3$ | ⋯ |

80 | $ 1.75 \xd7 10 \u2212 4$ | 2.29 | $ 1.61 \xd7 10 \u2212 3$ | 1.36 | |

160 | $ 3.22 \xd7 10 \u2212 5$ | 2.44 | $ 5.86 \xd7 10 \u2212 4$ | 1.46 | |

320 | $ 5.58 \xd7 10 \u2212 6$ | 2.53 | $ 2.00 \xd7 10 \u2212 4$ | 1.55 | |

CSL2R ( $ \alpha \u2273 4.001$) | 40 | $ 9.79 \xd7 10 \u2212 5$ | ⋯ | $ 1.53 \xd7 10 \u2212 4$ | ⋯ |

80 | $ 1.23 \xd7 10 \u2212 5$ | 2.99 | $ 1.93 \xd7 10 \u2212 5$ | 2.99 | |

160 | $ 1.53 \xd7 10 \u2212 6$ | 3.01 | $ 2.41 \xd7 10 \u2212 6$ | 3.00 | |

320 | $ 1.92 \xd7 10 \u2212 7$ | 2.99 | $ 3.01 \xd7 10 \u2212 7$ | 3.00 | |

CSL2R ( $ \alpha = 3.999$) | 40 | $ 8.99 \xd7 10 \u2212 4$ | ⋯ | $ 4.39 \xd7 10 \u2212 3$ | ⋯ |

80 | $ 2.00 \xd7 10 \u2212 4$ | 2.17 | $ 1.86 \xd7 10 \u2212 3$ | 1.24 | |

160 | $ 3.68 \xd7 10 \u2212 5$ | 2.44 | $ 6.73 \xd7 10 \u2212 4$ | 1.46 | |

320 | $ 7.23 \xd7 10 \u2212 6$ | 2.35 | $ 2.59 \xd7 10 \u2212 4$ | 1.38 |

Method . | N . | L_{1} error
. | L_{1} order
. | $ L \u221e$ error . | $ L \u221e$ order . |
---|---|---|---|---|---|

CSL2 | 40 | $ 9.79 \xd7 10 \u2212 5$ | ⋯ | $ 1.53 \xd7 10 \u2212 4$ | ⋯ |

80 | $ 1.23 \xd7 10 \u2212 5$ | 2.99 | $ 1.93 \xd7 10 \u2212 5$ | 2.99 | |

160 | $ 1.53 \xd7 10 \u2212 6$ | 3.01 | $ 2.41 \xd7 10 \u2212 6$ | 3.00 | |

320 | $ 1.92 \xd7 10 \u2212 7$ | 2.99 | $ 3.01 \xd7 10 \u2212 7$ | 3.00 | |

CSLR | 40 | $ 2.29 \xd7 10 \u2212 3$ | ⋯ | $ 1.19 \xd7 10 \u2212 2$ | ⋯ |

80 | $ 5.53 \xd7 10 \u2212 4$ | 2.05 | $ 4.75 \xd7 10 \u2212 3$ | 1.33 | |

160 | $ 1.29 \xd7 10 \u2212 4$ | 2.10 | $ 1.83 \xd7 10 \u2212 3$ | 1.38 | |

320 | $ 2.75 \xd7 10 \u2212 5$ | 2.23 | $ 6.91 \xd7 10 \u2212 4$ | 1.41 | |

CSL2T (any β, $ \alpha \u2273 4.001$) | 40 | $ 9.79 \xd7 10 \u2212 5$ | ⋯ | $ 1.53 \xd7 10 \u2212 4$ | ⋯ |

80 | $ 1.23 \xd7 10 \u2212 5$ | 2.99 | $ 1.93 \xd7 10 \u2212 5$ | 2.99 | |

160 | $ 1.53 \xd7 10 \u2212 6$ | 3.01 | $ 2.41 \xd7 10 \u2212 6$ | 3.00 | |

320 | $ 1.92 \xd7 10 \u2212 7$ | 2.99 | $ 3.01 \xd7 10 \u2212 7$ | 3.00 | |

CSL2T ( $ \beta = 4.0$, $ \alpha = 3.999$) | 40 | $ 8.50 \xd7 10 \u2212 4$ | ⋯ | $ 4.12 \xd7 10 \u2212 3$ | ⋯ |

80 | $ 1.75 \xd7 10 \u2212 4$ | 2.29 | $ 1.61 \xd7 10 \u2212 3$ | 1.36 | |

160 | $ 3.22 \xd7 10 \u2212 5$ | 2.44 | $ 5.86 \xd7 10 \u2212 4$ | 1.46 | |

320 | $ 5.58 \xd7 10 \u2212 6$ | 2.53 | $ 2.00 \xd7 10 \u2212 4$ | 1.55 | |

CSL2R ( $ \alpha \u2273 4.001$) | 40 | $ 9.79 \xd7 10 \u2212 5$ | ⋯ | $ 1.53 \xd7 10 \u2212 4$ | ⋯ |

80 | $ 1.23 \xd7 10 \u2212 5$ | 2.99 | $ 1.93 \xd7 10 \u2212 5$ | 2.99 | |

160 | $ 1.53 \xd7 10 \u2212 6$ | 3.01 | $ 2.41 \xd7 10 \u2212 6$ | 3.00 | |

320 | $ 1.92 \xd7 10 \u2212 7$ | 2.99 | $ 3.01 \xd7 10 \u2212 7$ | 3.00 | |

CSL2R ( $ \alpha = 3.999$) | 40 | $ 8.99 \xd7 10 \u2212 4$ | ⋯ | $ 4.39 \xd7 10 \u2212 3$ | ⋯ |

80 | $ 2.00 \xd7 10 \u2212 4$ | 2.17 | $ 1.86 \xd7 10 \u2212 3$ | 1.24 | |

160 | $ 3.68 \xd7 10 \u2212 5$ | 2.44 | $ 6.73 \xd7 10 \u2212 4$ | 1.46 | |

320 | $ 7.23 \xd7 10 \u2212 6$ | 2.35 | $ 2.59 \xd7 10 \u2212 4$ | 1.38 |

### B. Jiang-Shu test

^{13}

*u*(

*x*)

*=*

*1, the domain $ [ \u2212 1 , 1 ]$,*

*N*=

*200, $ \Delta t = 0.4 \Delta x , \u2009 \Delta x = 2 / N$, and periodic boundary conditions are used in this test. The initial condition is given as*

*a*=

*0.5, $ z = \u2212 0.7 , \u2009 \delta = 0.005 , \u2009 \alpha J S = 10$, and $ \beta J S = log ( 2 ) / ( 36 \delta 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)].

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 $ \beta \u2273 3$ is used in this test.

. | L_{1} error
. | $ L \u221e$ error . |
---|---|---|

CSL2 | $ 1.49 \xd7 10 \u2212 2$ | $ 2.46 \xd7 10 \u2212 1$ |

CSLR | $ 1.90 \xd7 10 \u2212 2$ | $ 2.60 \xd7 10 \u2212 1$ |

CSL2T (β = 3) | $ 1.27 \xd7 10 \u2212 2$ | $ 2.33 \xd7 10 \u2212 1$ |

CSL2T (β = 4) | $ 1.22 \xd7 10 \u2212 2$ | $ 2.33 \xd7 10 \u2212 1$ |

CSL2T (β = 5) | $ 1.19 \xd7 10 \u2212 2$ | $ 2.32 \xd7 10 \u2212 1$ |

CSL2T (β = 10) | $ 1.12 \xd7 10 \u2212 2$ | $ 2.31 \xd7 10 \u2212 1$ |

CSL2R | $ 1.28 \xd7 10 \u2212 2$ | $ 2.33 \xd7 10 \u2212 1$ |

. | L_{1} error
. | $ L \u221e$ error . |
---|---|---|

CSL2 | $ 1.49 \xd7 10 \u2212 2$ | $ 2.46 \xd7 10 \u2212 1$ |

CSLR | $ 1.90 \xd7 10 \u2212 2$ | $ 2.60 \xd7 10 \u2212 1$ |

CSL2T (β = 3) | $ 1.27 \xd7 10 \u2212 2$ | $ 2.33 \xd7 10 \u2212 1$ |

CSL2T (β = 4) | $ 1.22 \xd7 10 \u2212 2$ | $ 2.33 \xd7 10 \u2212 1$ |

CSL2T (β = 5) | $ 1.19 \xd7 10 \u2212 2$ | $ 2.32 \xd7 10 \u2212 1$ |

CSL2T (β = 10) | $ 1.12 \xd7 10 \u2212 2$ | $ 2.31 \xd7 10 \u2212 1$ |

CSL2R | $ 1.28 \xd7 10 \u2212 2$ | $ 2.33 \xd7 10 \u2212 1$ |

### C. Burgers' equation

*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.*

### D. Zalesak problem

^{47}in which a notched circle is rotated is widely used as a test of scalar advection schemes. The initial condition is given by

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 *L*_{1} 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.

Schemes . | L_{1} Error after one revolution
. | L_{1} error after four revolutions
. |
---|---|---|

CSL2 | $ 1.55 \xd7 10 \u2212 2$ | $ 2.21 \xd7 10 \u2212 2$ |

CSLR | $ 1.73 \xd7 10 \u2212 2$ | $ 2.75 \xd7 10 \u2212 2$ |

CSL2T (β = 4) | $ 1.28 \xd7 10 \u2212 2$ | $ 1.97 \xd7 10 \u2212 2$ |

CSL2R | $ 1.30 \xd7 10 \u2212 2$ | $ 1.99 \xd7 10 \u2212 2$ |

Schemes . | L_{1} Error after one revolution
. | L_{1} error after four revolutions
. |
---|---|---|

CSL2 | $ 1.55 \xd7 10 \u2212 2$ | $ 2.21 \xd7 10 \u2212 2$ |

CSLR | $ 1.73 \xd7 10 \u2212 2$ | $ 2.75 \xd7 10 \u2212 2$ |

CSL2T (β = 4) | $ 1.28 \xd7 10 \u2212 2$ | $ 1.97 \xd7 10 \u2212 2$ |

CSL2R | $ 1.30 \xd7 10 \u2212 2$ | $ 1.99 \xd7 10 \u2212 2$ |

### E. Lid-driven cavity flows

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 \xd7 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 $ \u2212 0.5 \u2264 x \u2264 0$ and $ \u2212 1 \u2264 y \u2264 \u2212 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 $ \beta = 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 $ \beta = 4.5$, the results are getting worse. Based on these results in Fig. 12, optimal *β* would be around $ 4 \u2264 \beta \u2264 5$. If around $ 4 \u2264 \beta \u2264 5$ is used, CSL2T is superior to CSL2R. In this paper, *β* = 4 is used hereafter.

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.

### F. Two droplets collisions and separations

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 set^{19,26} and volume-of-fluid^{11,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 model^{43,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 $ \Delta = D / 7.4$ was used, where D is the diameter of initial droplets. In this numerical resolution ( $ \Delta = 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 $ \Delta = D / 7.5$ and $ \Delta = 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.

### G. Droplet splashing

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 $ \rho liquid = 1000$, $ \rho air = 1.25$ kg/m^{3}, viscosities $ \mu liquid = 1.0 \xd7 10 \u2212 3$, $ \mu air = 1.82 \xd7 10 \u2212 5$ Pa s, surface tension $ \sigma = 7.2 \xd7 10 \u2212 2$ N/m, gravity 9.8 m/s^{2}, 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 \xd7 224 \xd7 48$ is used. The mesh size $ \Delta = 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.

## IV. SUMMARY

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author declares that there are no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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