Hybrid scheme for complex flows on staggered grids and application to multiphase flows

Highlights

• A hybrid scheme for complex flow simulations with staggered grids is developed.

• Reformulation of convective-terms-operators in a flux-based finite volume form.

• Unphysical oscillations are reduced by introducing a new interpolation technique.

• Robustness is increased by introducing hybrid weights, being continuous in time.

• Application to transient LES of the needle-opening process of a GDI injector.

Abstract

A hybrid scheme is developed that enables advanced simulations, such as multiphase flows, in complex flow setups with compressible flow solvers using a staggered arrangement of flow variables. For this purpose, a weighted essentially non-oscillatory (WENO) scheme is adapted for staggered grids and combined with a finite difference central scheme. Hybrid weights that maintain continuity in time are computed based on a limiter constructed from the deviation of the local weights from the optimal weights. A total variation diminishing (TVD) interpolation is introduced to avoid numerical oscillations from the flow field. Additionally, the finite difference operators are reformulated in a flux-based finite volume form for convective terms in order to avoid spurious oscillations for turbulent configurations. For demonstrating the accuracy of the newly developed scheme and its improvements compared to the existing schemes, a modified wavenumber analysis is shown. Furthermore, the new scheme is applied to single and multiphase test cases. Finally, it is used for simulating the needle-opening process of a gasoline direct injector (GDI) and predicting the amount of gas inside the nozzle.

Introduction

Multiphase flows often occur in complex flow situations, in which they interact with other phenomena such as shocks, turbulence, or high-frequency waves. Thus, numerical methods that are robust enough are required to deal with shocks or high-frequency waves, but conserve and accurately describe the level of turbulence and the interface between two phases, such as gas and liquid. This work presents a numerical algorithm to deal with such complex flow situations using shock-capturing methods as a starting point.

The application of shocks are observed across diverse fields such as aerodynamics, medical applications, or internal combustion engines. Lithotripsy, in which shock waves generated from collapsing microscopic cavitating bubbles break down kidney stones [1], is one of the examples applications in which shocks are deliberately induced. In other applications, such as fuel injectors [2] or propellers of ships [3], shocks are unintended and may cause structural damage. Therefore, accurate methods for predicting the impact of shocks and their common collateral phenomena, such as contact discontinuities and rarefactions, under various conditions are crucial. Experimental investigations are often difficult due to the very short time and length scales and hence, numerical simulations help to gain further understanding of shocks and their effect on relevant applications. However, also from a numerical point of view, shocks are difficult due to the resulting non-smooth flow field. More precisely, widely used numerical schemes such as the central difference scheme cannot handle the emerging discontinuities in the flow field, resulting in oscillations which destabilize a simulation, and thus numerical schemes that are suitable for shock capturing need to be investigated.

The Riemann problem, an initial value problem with piecewise constant data having a single discontinuity, is one of the common ways to numerically describe shocks [4] and has been extensively studied in the context of shock capturing schemes. For finite difference (FD) methods, correcting for upwinding as done by the first-order accurate Lax-Friedrichs scheme is one way to solve it. The Lax-Wendroff method [5] is an example for a second-order accurate method in both space and time. Flux limiting approaches such as the hybrid scheme by Harten and Zwas [6] were developed to improve the oscillating behavior of high-order methods. Other approaches for FD methods are summarized by Toro [4]. For finite volume (FV) methods, exact Riemann solvers such as Gudonov's method use non-linear root finding for every state pair across shocks in order to deal with the discontinuities. This results in computationally expensive schemes and has motivated the development of characteristics-based approximate Riemann solvers. Thus, modern shock-capturing schemes for FV methods rely on accurate estimations of numerical fluxes and approximation of variables at the cell boundaries from cell averaged ones called reconstruction.

The Harten-Lax-van Leer (HLL) solver [7] was the first efficient and robust approximate Riemann solver computing the numerical fluxes by assuming a particular wave structure. However, the assumption of a two-wave configuration is not valid for Euler equations and results in poor resolution of contact surfaces and shear waves. The Harten-Lax-van Leer-contact (HLLC) solver [8] is a more accurate extension of the HLL solver. In the HLLC solver, an intermediate wave corresponding to a contact discontinuity is restored, which leads to practically infinite resolution for stationary contact waves. Roe solvers [9] are the most widely used approximate Riemann solvers and have been used in different applications including compressible and multiphase flows. They are less dissipative but are computationally more intensive. Additionally, contact discontinuities are not resolved sharply, and the solvers may fail for low-density regions as they do not preserve positivity. This is a critical property when computing low-density and low-pressure flows, which, for example, appear in several compressible multi-component flow computations.

An example of reconstruction methods is the monotone upstream-centered scheme for conservation laws (MUSCL), which was introduced by van Leer [10] as the first high-order total variation diminishing (TVD) scheme. Instead of assuming piecewise constant data in the first-order Godunov method, high-order accuracy is obtained by reconstructing data and using limiters to avoid spurious oscillations. Harten et al. [11] extended MUSCL schemes to an arbitrary order of accuracy introducing essentially non-oscillatory schemes (ENO), in which the smoothest stencil is used across a discontinuity. Shu and Osher [12] suggested improvements in the implementation of ENO schemes, such as a procedure based on computing fluxes instead of cell averages, which led to simplification of the scheme for multiple dimensions. Liu et al. [13] introduced weighted essentially non-oscillatory (WENO) schemes. WENO schemes use a convex combination of the candidate stencils for computing interpolating polynomials to achieve ENO behavior rather than choosing the smoothest stencil for the reconstruction. With such an approach, an additional order of accuracy is obtained while reducing the cost of computing on parallel computers for WENO over ENO. In addition, smoother dependence on data may lessen some of the ENO scheme's oscillatory behavior near convergence. Later, Jiang and Shu [14] proposed an efficient implementation of a WENO scheme and further improved accuracy in smooth regions with new smoothness indicators. They also analyzed the computational costs of WENO schemes and found that a WENO scheme can be at least twice as fast as its corresponding ENO scheme on vector supercomputers.

The shock capturing schemes discussed above add artificial dissipation for stability. The added dissipation, however, is detrimental if used for computation of shocks in turbulent flows as it damps the turbulent fluctuations as well. For these cases, hybrid schemes combining central schemes in the smooth regions with upwind schemes across shocks have been introduced to limit the dissipative scheme to regions of discontinuities. Visbal and Gaitonde [15] presented a complex hybrid compact-Roe approach using a compact scheme in smooth regions and a robust shock capturing scheme in regions of shocks. The decoupling of the two schemes across the discontinuity is facilitated by a shock detector and inclusion of a few buffer points on both sides of the detected shock. Adams and Shariff [16] used a high-resolution hybrid compact-ENO scheme for studying shock-turbulence interactions. Pirozzoli [17] extended this non-conservative scheme to a hybrid conservative compact-WENO scheme. Hill and Pullin [18] developed a hybrid tuned center-difference (TCD)-WENO method for large-eddy simulation (LES) containing shocks. They computed the optimal weights of the WENO scheme in such a way that its dispersion properties match the TCD scheme used in smooth regions. This is done to avoid the occurrence of spurious oscillations where two schemes meet. Although the WENO scheme is supposed to provide high-order accuracy in smooth regions, the fluctuating turbulent field prevents the use of optimal weights resulting in unpredictable behavior. This necessitates the use of an explicit switch, which can be constructed based on the smoothness indicators of the WENO scheme. Ward and Pullin [19] described another hybrid strategy for compressible multi-component flows that uses a WENO scheme with a limiter. The limiter is proportional to the norm of the deviation of computed weights from the optimal weights. This approach allows for a smooth transition of a central difference scheme into a WENO scheme across the shock regions without any spurious oscillations.

Even in the absence of shocks, LES for complex geometries need to damp unresolved wavenumbers in order to remain stable and robust. An accurate numerical solution of the governing equations for LES requires high-order numerical methods so that the truncation errors from the numerical scheme do not interfere with the sub-grid scale model [20]. One way to remove the high-frequency wavenumbers, which cannot be resolved on the grid, is to use high-order compact filters [21]. However, the effect of such filters on turbulent flow structures is not understood clearly, and the interface between two phases is affected. Also, the order of the filter needs to be higher than the sixth order to not affect the smallest scales that are resolved on the grid [20], leading to a large stencil size, which is a problem at the boundaries, especially for complex geometries. Tyliszczak [22] examined the influence of explicit compact filters on the shear layer evolution. It was shown that compact filters remove the high-frequency wavenumbers regardless of their source. The inability of the filter to distinguish between the turbulent fluctuations and oscillations due to numerical errors is detrimental to the accuracy of the computed flow field. Since the numerical errors arising from the unresolved high-frequency wavenumbers are mainly due to the non-linear convective term, the application of a hybrid scheme instead of a compact filter is a promising way to achieve a physically consistent solution. As will be shown, this is helpful not only in the context of multiphase simulations for conserving sharp interfaces, but also for single phase setups, in which the compact filter might change the amplitude of turbulent fluctuations or remove flow features, such as vortices. The hybrid scheme can be used to damp any unphysical oscillations originating from the boundaries of complex geometry, extending the arbitrary high-order numerical framework to arbitrary geometries.

For low-Mach simulations, staggered arrangements of flow variables are popular since they avoid decoupling of pressure known as the checker-board effect and are more accurate for the same spatial resolution as their stencil widths for gradients and divergence are half the size of comparable collocated schemes. Moreover, the more compact stencil size of staggered schemes is advantageous at boundaries, making staggered arrangement suitable for complex geometries. However, most of the shock capturing and hybrid schemes in the literature have been developed for collocated grids. In this work, a hybrid scheme for spatially staggered flow variables is proposed, which has the advantages of staggered grids, especially for flow regimes with moderate Mach numbers but is also stable and robust. The hybrid scheme is a combination of an arbitrary order staggered FD central scheme developed by Morinishi et al. [23] and described in Desjardins et al. [24] and a WENO scheme from Jiang and Shu [14]. The two schemes are coupled across the discontinuity similar to the hybrid strategy of Ward and Pullin [19]. Straight-forward implementation of this hybrid scheme for a staggered arrangement of flow variables was found to give oscillatory solutions in numerical tests since the variables must be interpolated from face to the center and vice versa, typically done by centered interpolation operators. Hence, a MUSCL-type TVD interpolation is constructed and used for interpolation of the flow variables. Additionally, it was found that the continuity of weights of the WENO scheme must be maintained in time to avoid unphysical oscillations behind the discontinuities. Solutions are proposed to alleviate this problem. Furthermore, the FD operators are reformulated in a flux-based FV form for the convective terms of momentum and scalar equations including density and energy. So far, Trisjono et al. [25] showed that this reformulation is important for scalar equation operators to avoid spurious oscillations in the scalar fields in simulations of turbulent reactive flows. In this work, this formulation is extended to operators for vector equations.

The remainder of this article is organized as follows. A mathematical formulation of the shock-capturing problem, which is used for developing the hybrid scheme, and a short description of the used code framework are given in Section 2. Furthermore, the advantages of a staggered arrangement of flow variables for compressible flows are demonstrated. Existing numerical schemes, which are successfully used for shock capturing in the context of collocated flow variables, are reformulated for staggered flow variables and applied to Sod's problem in Section 3. Issues arising from the staggered formulation are emphasized. In Section 4, improvements for the staggered hybrid scheme are introduced, and a comparison to the existing schemes is shown. Then, the new staggered hybrid scheme is analyzed with respect to its wavenumber resolution and its performance with respect to single and multiphase examples in Section 5. In Section 6, an application to a gasoline direct injection (GDI) system is shown with focus on the needle opening process and the resulting multiphase flow. The paper finishes with conclusions in Section 7.