Improved hybrid approach of monotonic upstream-centered scheme for conservation laws and discontinuity sharpening technique for steady and unsteady flows

Fukushima, Gaku; Kitamura, Keiichi

doi:10.1063/5.0198163

In finite-volume methods, monotonic upstream-centered schemes for conservation laws (MUSCL) offer second-order spatial accuracy but tend to produce highly dissipative solutions for density discontinuity and weak shock waves. To address this limitation within a second-order framework, a novel strategy for hybridizing MUSCL with the tangent of the hyperbola interface capturing technique for both steady and unsteady compressible flows is presented. This hybridization optimizes the process based on the degree of nonlinearity and discontinuity around the target cells, providing a novel method to sharply resolve weak shock waves and robustly compute strong shock waves within the hybrid scheme. The proposed scheme sharply captures exceedingly weak shock waves that conventional MUSCL fails to resolve accurately due to excessive numerical dissipation. Furthermore, for resolving small vortices induced by instability at slip lines, computational results demonstrate high-resolution surpassing fifth-order spatial accuracy schemes within this second-order spatial accuracy framework with less computational cost. Moreover, the scheme exhibits commendable convergence and robustness when applied to steady-state problems featuring strong shock waves. This scheme offers a more precise and high-resolution alternative to conventional MUSCL for compressible flow computations, as it requires no additional stencil for reconstruction, unlike conventional fifth-order schemes.

I. INTRODUCTION

Compressible fluid computations employing shock-capturing schemes within finite-volume methods (FVM) have found widespread use across various applications. The most commonly employed shock-capturing method is the monotonic upstream-centered scheme for conservation laws (MUSCL) with slope limiters,¹ typically providing second-order spatial accuracy. MUSCL is particularly favored in aerospace engineering^2–4 due to its straightforward mathematical expression and accuracy in fundamental flows featuring shock waves. However, it often yields dissipative solutions for complex flows. Even in cases involving shock wave computations, when the shock wave is extremely weak (e.g., shock Mach number ≈ 1.01), it is resolved as a continuous wave-like shape, losing its discontinuous profile.⁵ In multiphase flow simulations, gas–liquid interfaces suffer from excessive dissipation.⁶ There is a lack of widely accepted second-order schemes capable of effectively and accurately handling density discontinuities and weak shock waves.

One approach to achieve accurate solutions is through the utilization of higher-order methods that incorporate multiple surrounding cells for reconstructing target cells. For instance, the weighted essentially non-oscillatory (WENO)⁷ scheme provides less dissipative numerical solutions, achieving a maximum fifth-order spatial accuracy by assessing stencil smoothness, typically considering five cell values. Recent developments have seen various methods based on the WENO scheme emerge.^8–11 In an alternative approach, a fifth-order spatial accuracy shock-capturing scheme based on the boundary variable diminishing (BVD) principle^12,13 has been developed.¹⁴ This scheme selectively incorporates cell boundary values through a fourth-degree polynomial and the tangential hyperbolic interface capturing (THINC) function.^15–17 Numerical tests using these higher-order schemes have yielded promising results; however, they entail greater complexity and expense compared to the widely used second-order spatial accuracy schemes. We maintain that second-order spatial accuracy remains the most suitable choice for broad acceptance in practical applications due to its simplicity, typically relying on neighboring cells for reconstruction. Furthermore, extending the second-order framework to discretization using unstructured grids requires less effort compared to higher-order schemes. From this perspective, our study introduces a novel second-order scheme based on MUSCL.

To address the limitations of the dissipative nature of MUSCL in multiphase flow computations, THINC has been utilized. THINC is a method used to express a density interface in the volume-of-fluid method.¹⁵ Its distinct feature, which can express discontinuity by the tangent hyperbola function, is also effective for the reconstruction process in the FVM. Nonomura et al.⁶ introduced THINC for multiphase computations using a two-fluid model¹⁸ within a finite-volume MUSCL framework. They used THINC only in the reconstruction of the volume fraction and obtained sharp gas–liquid and gas–gas interfaces. Deng et al. proposed the MUSCL-THINC-BVD scheme,¹⁹ determining cell boundary values based on the BVD principle to minimize numerical dissipation. More recently, Chiu et al.²⁰ introduced a hybrid MUSCL-THINC scheme in which the volume fraction is reconstructed using THINC, while other primitive variables are reconstructed by weighting cell boundary values computed through both MUSCL and THINC. This scheme displayed low-dissipation results for single-phase unsteady flows and two-phase problems.

Given the success of THINC in multiphase flow computations, its application appears promising for single-phase compressible flow simulations aiming to accurately capture density discontinuities and weak shock waves. However, prior research on schemes combining MUSCL and THINC for solving single-phase steady and unsteady compressible flows remains limited, despite the prevalence of applications involving multiphase flow computations, which often involve unsteady flows. The only prior study on this hybrid approach was the original hybrid MUSCL-THINC scheme by Chiu et al.²⁰ While it has been reported that this scheme can provide low-dissipative solutions in single-phase compressible flows, the numerical examples examined in Ref. 20 were limited to a few unsteady flow cases (shock-tube and double-Mach reflection problems). Furthermore, its performance in terms of convergence for steady-state problems, crucial for practical applications, was not discussed. In fact, as we will later report in this paper, the original hybrid MUSCL-THINC scheme exhibited severe oscillations in steady problems with strong shock waves. Additionally, the scheme suffered from shock anomalies in the presence of strong shock waves, known as, the carbuncle phenomenon.^21–24 Therefore, the original approach in Ref. 20 may not be universally suitable for single-phase steady and unsteady compressible flow computations. As such, we believe that the hybridization strategy of MUSCL and THINC for single-phase compressible flow computations needs to be thoroughly investigated and reformulated.

It should be noted that the BVD principle^15–17 is another possible way to effectively incorporate THINC. However, schemes relying on the BVD principle tend to be intricate, as they require a selection procedure after calculating cell boundary values. This complexity was confirmed in prior research,²⁰ where the hybrid MUSCL-THINC scheme was found to be more efficient than the MUSCL-THINC-BVD scheme. Therefore, for the sake of simplicity and efficiency, we did not employ the BVD principle in this study.

In this paper, we propose a novel hybridization strategy of MUSCL and THINC aimed at achieving a high-resolution and robust second-order scheme for single-phase steady and unsteady compressible flows. Specifically, the weights assigned to cell boundary values computed using MUSCL and THINC have been optimized based on the degree of nonlinearity and discontinuity within the stencil. This framework is implemented through the introduction of the “nonlinearity-weighted parameter” and the “slope ratio-weighted parameter.” The former is sensitive to the magnitude of nonlinearity in the phenomenon, reflecting its self-sharpening characteristics. This parameter is newly introduced in this study and plays a crucial role in obtaining stable solutions for strong shock waves. The latter, sensitive to discontinuities, was initially proposed in the original hybrid MUSCL-THINC scheme.²⁰ While the slope ratio-weighted parameter is borrowed from the original hybrid MUSCL-THINC scheme,²⁰ the proposed method primarily relies on the idea that the scheme automatically attains a physically consistent balance between the phenomenon's self-sharpening nature and the reconstruction process. This feature is crucial for achieving accurate computations across a broad spectrum of single-phase steady and unsteady compressible flow conditions. Thus, the proposed scheme can be viewed as distinct from the original hybrid MUSCL-THINC scheme, in which the weight is determined primarily based on the slope ratio, and the primary focus is on unsteady multiphase flows.

The remainder of this paper is structured as follows: Sec. II provides an overview of the calculation method employed in this study, outlining the reconstruction methods and presenting the proposed scheme. Section III details the results of accuracy validation and numerical experiments for both steady and unsteady problems, including a discussion of computational costs. Finally, Sec. IV concludes this study.

II. NUMERICAL SETUP

A. Governing equations

The governing equations are the two-dimensional compressible Euler equations or the Navier–Stokes equations, as follows:

\begin{matrix} \frac{\partial Q}{\partial t} + \frac{\partial F_{k}}{\partial x_{k}} = \frac{\partial F v_{k}}{\partial x_{k}}, \\ Q = [\begin{matrix} ρ \\ ρ u_{l} \\ ρ E \end{matrix}], F_{k} = [\begin{matrix} ρ u_{k} \\ ρ u_{l} u_{k} + p δ_{l k} \\ ρ u_{k} H \end{matrix}], {Fv}_{k} = [\begin{matrix} 0 \\ τ_{l k} \\ u_{m} τ_{m k} + κ \frac{\partial T}{\partial x_{k}} \end{matrix}], τ_{l k} = μ (\frac{\partial u_{l}}{\partial x_{k}} + \frac{\partial u_{k}}{\partial x_{l}}) - \frac{2}{3} μ \frac{\partial u_{n}}{\partial x_{n}} δ_{l k}, \end{matrix}

(1)

where Q represents the vector of conservation variables, ρ denotes the density, u_i stands for the velocity component in Cartesian coordinates, E denotes the total energy, p represents the pressure, H denotes the total enthalpy (H = E + p/ρ), and T denotes the temperature. F_k and Fv_k signify the inviscid and viscous numerical fluxes, respectively. If Fv_k = 0, Eq. becomes the Euler equations (inviscid). The working gas was assumed to be air, which is a calorically perfect gas with a specific heat ratio γ = 1.4. The Prandtl number, Pr, was 0.72. The molecular viscosity μ (given by Sutherland's law) and thermal conductivity κ are related by κ = c_pμ/Pr, where c_p denotes the specific heat at constant pressure.

δ_{l k}

represents the Kronecker delta. The subscripts k, l, m = 1 and 2 are in 2D.

The equations were discretized using cell-centered FVM. For time integration, the four-stage Runge–Kutta²⁵ method was employed in this study. SLAU2²⁶ was utilized to compute the inviscid numerical flux at the cell boundary. Viscous flux computations were performed using a second-order central difference.

B. Reconstruction methods

This section provides an overview of the reconstruction methods employed in this study. For simplicity, a one-dimensional equally discretized computational domain is considered. The target cell, denoted as the ith cell, possesses a width defined as [x_i_−1/2, x_i_+1/2], with the mesh size represented as Δx = x_i_+1/2 – x_i_−1/2. Figure 1 schematically illustrates the reconstructed variables within the ith cell. The primitive variable q, comprising density (ρ), velocities (u and v), and pressure (p), within the ith cell is reconstructed by utilizing both its intrinsic information and that of its neighboring cells; specifically, information from three cells is considered. Consequently, a stencil [i – 1, i, i + 1] is employed to reconstruct the ith cell, resulting in the determination of cell boundary values denoted as q_R_,_i_–1/2 and q_L_,_i_+1/2. Here, the subscripts i ± 1/2 denote the cell boundaries between the ith cell and its adjacent cells, with L and R signifying the left and right sides of the boundary, respectively. The superscripts M and T denote the reconstructed cell boundary values by MUSCL and THINC, respectively. In the proposed approach, the cell boundary values computed by the two methods are hybridized depending on flow characteristics. The specific reconstruction functions and cell boundary values will be described later.

FIG. 1.

View large Download slide

Illustration of reconstruction and cell boundary values.

1. MUSCL

First, MUSCL¹ is introduced. MUSCL expresses the intercell property distribution as²⁷^,

q (x) = q_{i} + \frac{1}{Δ x} (x - x_{i}) δ_{i} q + \frac{3 η}{2 Δ x^{2}} [(x - x_{i}) - \frac{Δ x^{2}}{12}] δ_{i}^{2} q,

(2)

where

q_{i} = \int_{x - 1 / 2}^{x + 1 / 2} q (x) dx, δ_{i} q = \frac{q_{i + 1} - q_{i - 1}}{2}, δ_{i}^{2} q = q_{i + 1} - 2 q_{i} + q_{i - 1},

(3)

and

η

denotes the parameter of MUSCL. η = 1/3 was used in this study. In this study, the minmod limiter,²⁸

ψ (r) = max [0, min (1, r)], r = \frac{q_{i + 1} - q_{i}}{q_{i} - q_{i - 1}},

(4)

was employed for slope limiting. Using Eqs. , we obtained the cell boundary value to compute the numerical flux as follows:

\begin{matrix} q_{L, i + 1 / 2}^{M} = q_{i} + \frac{1}{4} [(1 - η) ψ (r) (q_{i} - q_{i - 1}) + (1 + η) ψ (\frac{1}{r}) (q_{i + 1} - q_{i})], \\ q_{R, i - 1 / 2}^{M} = q_{i} - \frac{1}{4} [(1 + η) ψ (r) (q_{i} - q_{i - 1}) - (1 - η) ψ (\frac{1}{r}) (q_{i + 1} - q_{i})] . \end{matrix}

(5)

Hereafter, the aforementioned MUSCL approach (MUSCL parameter η = 1/3 and minmod limited) is denoted as MUSCL for simplicity of notation. This combination (primitive variables reconstructions by η = 1/3 MUSCL with the minmod limiter) has been widely used in previous studies.^2,20

2. THINC

THINC was originally developed to express an interface between two fluids in the volume-of-fluid method^15–17 and was extended to the FVM.^6,19,20 The formula for FVM was used in this study. In the FVM, reconstruction using THINC can be defined if the local profile within the stencil is monotonic. In the non-monotonic profile, the reconstruction by THINC is not considered. In the monotonic case, the profile within the target cell is expressed as follows:

q (x) = q_{min} + \frac{Δ q}{2} \{1 + θ tanh [β (X_{i} - d_{i})]\},

(6)

where

\begin{matrix} \begin{matrix} \begin{matrix} q_{min} = min (q_{i - 1}, q_{i + 1}), \\ Δ q = |q_{i + 1} - q_{i - 1}|, \\ θ = sgn (q_{i + 1} - q_{i - 1}), \end{matrix} \\ X_{i} = \frac{x - x_{i - 1 / 2}}{x_{i + 1 / 2} - x_{i - 1 / 2}}, \\ d_{i} = \frac{1}{2 β} \ln \frac{1 - A}{1 + A}, \end{matrix} \\ A = \frac{B / cosh (β) - 1}{tanh (β)}, \\ B = exp \{θ β [\frac{2 (q_{i} - q_{min} + ε)}{Δ q + ε} - 1]\}, \end{matrix}

(7)

where d_i denotes the center of the jump location, ε represents a small real number to prevent division by zero and ε = 10⁻²⁰ in this paper. β denotes the sharpness parameter in THINC. The THINC function can imitate a steep profile of the variables by taking a large value of β. This is a user-defined parameter in THINC and the proposed hybrid approach. To eliminate this degree of freedom, a suitable value of β for compressible flow simulations is discussed later. The THINC function can only be defined if the local profiles of the variables in the stencil are monotonic. THINC provides the cell boundary value as

\begin{matrix} q_{L, i + 1 / 2}^{T} = q_{min} + \frac{Δ q}{2} [1 + θ \frac{tanh (β) + A}{1 + A tanh (β)}], \\ q_{R, i - 1 / 2}^{T} = q_{min} + \frac{Δ q}{2} (1 + θ A) . \end{matrix}

(8)

3. Hybrid MUSCL-THINC scheme (for single-phase flows)

In the case of two-phase flows, the hybrid MUSCL-THINC scheme²⁰ reconstructs the volume fraction exclusively using THINC, while the other primitive variables are reconstructed by weight-averaging the minmod-limited MUSCL and THINC. However, for single-phase flows where a volume fraction is absent, the reconstruction is solely performed using primitive variables. Subsequently, the single-phase variant of the hybrid MUSCL-THINC scheme will be referred to as the “original hybrid MUSCL-THINC scheme” to clearly distinguish it from the proposed scheme. In the original hybrid MUSCL-THINC scheme, the cell boundary values of MUSCL and THINC are blended when the profile of the primitive variable is monotonic, defined as

(q_{i + 1} - q_{i}) (q_{i} - q_{i - 1}) > ϵ .

(9)

Here,

ϵ

represents a small real number

ϵ

= (1 × 10⁻¹⁵)² = 1 × 10⁻³⁰ to prevent reconstruction by THINC in cases of non-monotonic variable profiles. The cell boundary value in the original hybrid MUSCL-THINC scheme is expressed as

\begin{matrix} q_{L, i + 1 / 2}^{M - T} = (1 - ζ_{i}) q_{L, i + 1 / 2}^{M} + ζ_{i} q_{L, i + 1 / 2}^{T}, \\ q_{R, i - 1 / 2}^{M - T} = (1 - ζ_{i}) q_{R, i - 1 / 2}^{M} + ζ_{i} q_{R, i - 1 / 2}^{T}, \end{matrix}

(10)

where the superscript M–T represents the original hybrid MUSCL-THINC. Subscripts L and R indicate the cell boundary values for each method. The slope ratio-weighted parameter

ζ_{i}

⁠, which determines the ratio of MUSCL to THINC, is defined as follows:

ζ_{i} \equiv 1 - min (\frac{q_{L, i + 1 / 2}^{M} - q_{R, i - 1 / 2}^{M}}{q_{i + 1} - q_{i} + ε}, \frac{q_{L, i + 1 / 2}^{M} - q_{R, i - 1 / 2}^{M}}{q_{i} - q_{i - 1} + ε}) .

(11)

When the local profile of variables within the stencil remains continuous,

ζ_{i}

approaches a value close to 0, and the MUSCL method is applied. Conversely, when the stencil exhibits a discontinuity,

ζ_{i}

approaches a value close to 1, leading to the utilization of THINC reconstruction. The proposed scheme incorporates this slope-ratio-weighted parameter to assess the discontinuity in the variable profile within the stencil. In cases where the profile is non-monotonic, the cell boundary values obtained by MUSCL are employed. This reconstruction process is executed at each sub-time step using the Runge–Kutta method.

4. Novel hybridization approach of MUSCL and THINC

Previous studies have indicated that the selection of the reconstruction process should be contingent on the flow properties. For instance, correctly resolving a very weak shock wave necessitates a reconstruction process that minimizes numerical dissipation because the nonlinearity inherent to a weak shock wave can be easily offset by numerical dissipation.⁵ This requirement is also applicable to density discontinuities.^6,20 Conversely, a strong shock wave, being a strong nonlinear phenomenon, maintains a discontinuous profile without special reconstruction processes. Thus, we advocate avoiding a nonlinear reconstruction process using the discontinuity sharpening technique to mitigate potential numerical oscillations. According to literature,²⁶ a moderate amount of numerical dissipation is preferable to prevent shock anomalies. Consequently, the reconstruction process must be optimized to align with the fluid properties within the stencil. Notably, this approach has not been previously implemented within a second-order FVM framework.

This paper introduces a numerical scheme built upon the concept of adjusting the balance of nonlinearity between the phenomenon itself and the reconstruction function to highly resolve weak nonlinear phenomena and robustly compute strong shock waves. To realize this idea within a second-order framework, we have devised a new hybrid scheme combining MUSCL and THINC, guided by the following principles.

MUSCL is employed for continuous flows, while THINC is reserved for discontinuities.
For weak nonlinear phenomena, such as weak shock waves, density discontinuities, and entropy waves, which can be easily smoothed by numerical dissipations of conventional MUSCL, the reconstruction by THINC is used to oppose the numerical dissipation.
For strong shock waves, only MUSCL is employed.

Criterion (A) is achieved by borrowing the slope-ratio-weighted parameter $ζ_{i}$ ⁠, a feature from the original hybrid MUSCL-THINC scheme described in Sec. II B 3, which is effective in detecting discontinuities.²⁰ To satisfy requirements (B) and (C), novel aspects of the proposed method, we have introduced a straightforward parameter designed to limit the use of THINC near strongly nonlinear regions. This parameter is termed the “nonlinearity-weighted parameter” and is employed to assign additional weight to both MUSCL and THINC. The parameter assumes values close to 1 in the vicinity of regions with weak nonlinearity and approaches 0 only near strongly nonlinear regions, thereby diminishing the weight assigned to THINC.

Here, the calculation procedure of the nonlinearity-weighted parameter at the cell boundary i + 1/2 is described as an example. Pressure and density are the key variables used to compute the nonlinearity-weighted parameter at the cell boundary. Specifically, the relevant following variables are computed:

\begin{matrix} \begin{matrix} d p / d x = (p_{i + 1} - p_{i}) / Δ x, \\ d ρ / d x = (ρ_{i + 1} - ρ_{i}) / Δ x, \\ ϕ_{i + 1 / 2}^{p} = \frac{max (p_{i}, p_{i + 1})}{min (p_{i}, p_{i + 1})}, \end{matrix} \\ ϕ_{i + 1 / 2}^{ρ} = \frac{max (ρ_{i}, ρ_{i + 1})}{min (ρ_{i}, ρ_{i + 1})}, \end{matrix}

(12)

where

d p / d x, d ρ / d x, ϕ_{i + 1 / 2}^{p}

⁠, and

ϕ_{i + 1 / 2}^{ρ}

denote the pressure gradient, density gradient, pressure ratio, and density ratio at the cell boundary of i + 1/2, respectively. These variables are used to ascertain whether the cell boundary represents a shock or a density discontinuity. The relationship between pressure and density gradients, specifically the sign of

(d p / d x) (d ρ / d x)

⁠, is utilized, with multiplication used to prevent division by zero.

If

(d p / d x) (d ρ / d x)

is positive, this indicates that the cell boundary may exhibit various states, including continuous flow, shock waves, or density discontinuities. In Fig. 2, the example of a shock wave case [Fig. 2(a)] and a density discontinuity or entropy wave component case [Fig. 2(b)] are illustrated. In such cases, the value of the pressure ratio over the density ratio,

ϕ_{i + 1 / 2}^{p ρ} = ϕ_{i + 1 / 2}^{p} / ϕ_{i + 1 / 2}^{ρ},

(13)

is calculated. This value,

ϕ_{i + 1 / 2}^{p ρ}

⁠, gives the temperature ratio, and in the case of isentropic flow or shock wave, it takes a value higher than 1. According to the Rankine–Hugoniot relations, the stronger the shock wave, the larger

ϕ_{i + 1 / 2}^{p ρ}

becomes. Thus, if the stencil includes a large

ϕ_{i + 1 / 2}^{p ρ}

⁠, the nonlinearity-weighted parameter should be small to avoid the use of THINC. This fulfills requirement (C). On the other hand, if the density ratio

ϕ_{i + 1 / 2}^{ρ}

is larger than the value expected from isentropic flow or shock waves as shown in Fig. 2(b),

ϕ_{i + 1 / 2}^{p ρ}

takes the value around 1 or smaller. In such cases, the cell boundary is likely dominated by weak nonlinear phenomena, such as weak shock waves, density discontinuities, or entropy wave components. In these stencils, THINC should be employed to achieve a less dissipative solution to fulfill requirement (B). Consequently, the nonlinearity-weighted parameter should approximate one. In instances where

(d p / d x) (d ρ / d x)

is negative, the field represents an unphysical state. However, it is worth noting that this field is observed around density discontinuities and entropy waves, typically arising due to numerically generated transitional states. Thus, in these regions, THINC reconstruction is applied to minimize numerical dissipation.

FIG. 2.

View large Download slide

Schematic explanation of distinguishing strong shock waves from weak ones, density discontinuity, and entropy wave components at cell boundary using pressure and density. Example cases for (a) shock wave and (b) density discontinuity or entropy wave components are shown.

Based on the above considerations, we can conclude that only if

ϕ_{i + 1 / 2}^{p ρ}

takes a higher value, the cell boundary corresponds to the strong shock wave. Thus, we define an exponential function that approaches zero near strong shock waves and approaches approximately one near weak shock waves, isentropic flows, density discontinuities, and entropy wave components by using the variable

ϕ_{i + 1 / 2}^{p ρ}

⁠, as follows:

ξ_{i + 1 / 2} = \{\begin{matrix} 1 & (d p / d x) (d ρ / d x) < 0, \\ exp \{- C [max (1, ϕ_{i + 1 / 2}^{p ρ}) - 1)]\} & (d p / d x) (d ρ / d x) \geq 0, \end{matrix}

(14)

where

ξ_{i + 1 / 2}

is the nonlinearity weighted parameter of the cell boundary at x_i_+1/2.

C

is the constant that controls the nonlinearity-weighted parameter. In this paper,

C

= 25 was adopted: we confirmed that this value showed good performance for shock wave computations through numerical tests. The nonlinearity-weighted parameter value as a function of

ϕ_{i + 1 / 2}^{p ρ}

and the shock Mach number M_s is shown in Fig. 3. The plot against M_s was obtained by assuming that the cell boundary expressed a shock wave without thickness. That is, the following relationship from the Rankine–Hugoniot relations was applied to convert

ϕ_{i + 1 / 2}^{p ρ}

to M_s:^29,30

ϕ_{i + 1 / 2}^{p ρ} = \frac{(2 γ M_{s}^{2} - γ + 1) \{(γ - 1) M_{s}^{2} + 2\}}{{(γ + 1)}^{2} M_{s}^{2}} .

(15)

FIG. 3.

View large Download slide

Output of nonlinearity weighted parameter as a function of (a) $ϕ_{i + 1 / 2}^{p ρ}$ and (b) M_s. (b) is obtained assuming that the variables follow Rankine–Hugoniot relations.

According to Fig. 3(b), we can observe the behavior of the nonlinearity-weighted parameter against M_s, which is a universal indicator of shock wave strength. If the cell boundary expresses a shock wave without internal points, the nonlinearity-weighted parameter becomes almost zero for a shock wave stronger than M_s ≈ 1.3. The shock wave usually has a certain thickness and thus has a smaller $ϕ_{i + 1 / 2}^{p ρ}$ at the internal cell boundaries than the value computed from M_s. We confirmed that stopping usage of THINC for the shock wave stronger than this level is effective through preliminary numerical tests.

There are two choices of the nonlinearity-weighted parameter for the ith cell:

ξ_{i - 1 / 2}

and

ξ_{i + 1 / 2}

⁠. To ensure a safer value, a smaller value was adopted. Thus, the nonlinearity-weighted parameter in the ith cell

ξ_{i}

is expressed as

ξ_{i} = min (ξ_{i - 1 / 2}, ξ_{i + 1 / 2}) .

(16)

A flow chart outlining the proposed scheme is presented in Fig. 4. Prior to the primitive variable reconstruction loop, a loop for calculating the nonlinearity-weighted parameter

ξ_{i}

is implemented. Subsequently, following the computation of

ξ_{i}

⁠, cell boundary values are determined using MUSCL [Eq. ] and THINC (if the variable profile in the stencil is monotonic) [Eq. ], along with the slope ratio-weighted parameter

ζ_{i}

[Eq. ], as inherited from the original hybrid MUSCL-THINC scheme. These values are then hybridized via weighted averaging. In the proposed scheme, the slope ratio-weighted parameter

ζ_{i}

⁠, which is sensitive to discontinuity in the stencil is multiplied by the nonlinearity-weighted parameter

ξ_{i}

⁠, and

ζ_{i} ξ_{i}

serves as the weight. Consequently, if the variable profile within the stencil is monotonic, the cell boundary values can be computed as follows:

\begin{matrix} q_{L, i + 1 / 2}^{M - T - p ρ} = (1 - ζ_{i} ξ_{i}) q_{L, i + 1 / 2}^{M} + ζ_{i} ξ_{i} q_{L, i + 1 / 2}^{T}, \\ q_{R, i - 1 / 2}^{M - T - p ρ} = (1 - ζ_{i} ξ_{i}) q_{R, i - 1 / 2}^{M} + ζ_{i} ξ_{i} q_{R, i - 1 / 2}^{T} . \end{matrix}

(17)

FIG. 4.

View large Download slide

Flowchart of the proposed scheme.

MUSCL alone was employed when the variable profile within the stencil exhibited non-monotonic behavior. The proposed scheme is designated as the “hybrid MUSCL-THINC-pρ” (abbreviated as “T-MUSCL”) scheme, where “pρ” denotes the usage of pressure (p) and density (ρ) for computing the nonlinearity-weighted parameter $ξ_{i}$ ⁠. The scheme name is denoted by the superscript “M-T-pρ” in Eq. (17). According to this formulation, the reconstruction reflecting discontinuity characteristics is achieved: THINC is exclusively applied to discontinuous and weakly nonlinear flow components; otherwise, MUSCL is utilized. As such, the scheme satisfies all the outlined requirements (A)–(C). The functionality of the weight parameter $ξ_{i}$ in the numerical example will be discussed in detail in Sec. III.

III. NUMERICAL TEST

In Secs. III A–III J, we present the numerical results of the benchmark tests. To facilitate a clear distinction between the original hybrid MUSCL-THINC scheme by Chiu et al.²⁰ and the newly introduced hybrid MUSCL-THINC-pρ scheme, we will refer to these schemes as “MUSCL-THINC” and “MUSCL-THINC-pρ,” respectively.

A. Weak shock formation problem (unsteady, accuracy study)

Initially, we need to determine the suitable parameter value for the THINC function. When dealing with a stencil containing a monotonous variable profile, THINC can be used to reconstruct the inner-cell variable assuming a discontinuity expressed as a function of the sharpness parameter β, as represented in Eq. (6). For density and volume fraction (appears in the two-fluid model) at an interface, a large value of β is preferred.^6,20 We believe that this is because the interface did not exhibit self-sharpening. However, in the context of compressible flows and shock waves, it is imperative to carefully select β because shock waves have physically consistent sharpening characteristics. Thus, as a first step, we must identify an adequate variable β in MUSCL-THNC and MUSCL-THINC-pρ. To achieve this, the weak shock wave formation problem was solved. Figure 5 illustrates the shock wave formation problem described in.³⁰ This problem has been used in experimental studies to estimate the shock formation distance;^31,32 however, this is the first time it has been used in numerical simulations. A certain state of the continuous compression waves at t = 0 was considered. We define x_a as the location of the trailing edge of the compression wave, and x_b as the leading edge. Because the wave propagation velocity at the trailing edge of the compression wave is greater than that at the leading edge, the trailing wave catches up at the leading edge. Subsequently, a shock wave is formed at a certain point x_s.^29,30 By successfully addressing this process for a weak compression wave, the scheme can accurately model weak, nonlinear, yet discontinuous fluid properties without introducing excessive numerical dissipation.

FIG. 5.

View large Download slide

Schematics of shock formation problem.

The calculations were conducted within a one-dimensional space ranging from [0, 2500], divided equally into 2500 cells (Δx = 1). Initially, the compression wave was situated between x_a = 135 and x_b = 165, resulting in an initial thickness of δ₀ = x_b − x_a = 30. Because the cell size is unity, the thickness corresponds to the number of cells. In this problem, the conditions behind the compression wave corresponded to those behind the shock wave. The primitive variables within the compression wave were determined by assuming a linear distribution of the shock Mach number, ranging from 1 to M_s. This local shock Mach number within the compression wave, denoted as

M_{s}^{'}

⁠, is defined as

M_{s}^{'} = M_{s} - (M_{s} - 1) \frac{x - x_{a}}{x_{b} - x_{a}} (x_{a} \leq x \leq x_{b}) .

(18)

The primitive variables within the compression wave can be expressed as follows:

(ρ, u, p) = (\frac{(γ + 1) M_{s}^{' 2}}{(γ - 1) M_{s}^{' 2} + 2} ρ_{R}, \frac{2 a_{R}}{γ + 1} (M_{s}^{'} - \frac{1}{M_{s}^{'}}) + u_{R}, [1 + \frac{2 γ}{γ + 1} (M_{s}^{' 2} - 1)] p_{R}) (x_{a} \leq x \leq x_{b}),

(19)

where a represents the speed of sound. In the other domain, the initial conditions were given as follows:

\begin{matrix} (ρ, u, p) = \{\begin{matrix} \begin{matrix} (ρ_{L}, u_{L}, p_{L}) & (0 \leq x < x_{a}) \end{matrix} \\ \begin{matrix} (ρ_{R}, u_{R}, p_{R}) & (x_{b} < x \leq 2500) \end{matrix} \end{matrix}, \\ (ρ_{L}, u_{L}, p_{L}) = (\frac{(γ + 1) M_{s}^{2}}{(γ - 1) M_{s}^{2} + 2} ρ_{R}, \frac{2 a_{R}}{γ + 1} (M_{s} - \frac{1}{M_{s}}) + u_{R}, [1 + \frac{2 γ}{γ + 1} (M_{s}^{2} - 1)] p_{R}), \\ (ρ_{R}, u_{R}, p_{R}) = (1, 0, 1) . \end{matrix}

(20)

The shock Mach number was set as M_s = 1.01. The sharpening parameter in THINC β varied from 1.6 to 3.2 in increments of 0.4. The CFL number was set to 0.5 (Δt ≈ 0.51) for time evolution.

The shock formation distance L_s is defined as the distance from the initial position of the leading edge to the point where the trailing-edge wave overtakes the leading edge, that is L_s = x_s − x_b. We can estimate L_s based on the relationship between the characteristic velocities of the leading and trailing edges of the compression waves. Specifically, L_s is expressed by the following equation:³⁰

L_{s} = \frac{c_{+, b}}{c_{+, a} - c_{+, b}} δ_{0},

(21)

where c_+,_a denotes the right characteristic wave velocity at x_a (c_+,a= u_a + a_a). In this setup, L_s = 1507. Thus, the shock wave was assumed to form at x_s = x_b + L_s = 1672. The thickness of the compression waves was expected to decrease linearly before integration and become zero at x_s.

The space-time (x–t) diagram depicting density contours for MUSCL and MUSCL-THINC-pρ (with β = 2.4) is presented in Fig. 6. In this case, the result for MUSCL-THINC is not shown because, under these weak nonlinear conditions, the nonlinearity-weight parameter approaches 1, making no significant difference between MUSCL-THINC-pρ and MUSCL-THINC evident. The right panels provide enlarged views. Since the contours were generated using results every ten steps, the time resolution of the contours is approximately t = 5.1 (= Δt × 10 steps). The contour plots reveal the propagation of a compression wave from left to right. In the MUSCL case, there is no transition from a compression wave to a shock wave, and the compression wave profile remains unchanged during propagation, as demonstrated in the enlarged views in Fig. 6(a). In contrast, a transition to a shock wave appears to occur in MUSCL-THINC-pρ: the thickness of the compression wave diminishes with propagation, eventually leading to the formation of a discontinuous profile. This process is confirmed by the enlarged view in Fig. 6(b). The specific definition of the formation of a “discontinuous shock front” will be discussed later.

FIG. 6.

View large Download slide

Space (x) – time (t) diagram of shock formation problem solved using (a) MUSCL and (b) MUSCL-THINC-pρ.

Now, turning to the discussion of the density profile for different values of β. Figure 7 shows the density profiles at t = 500 and t = 1520 using MUSCL-THINC-pρ and MUSCL. At t = 500 [Fig. 7(a)], the compression wave reaches approximately x = 645–665, where a shock wave has not yet formed. MUSCL-THINC-pρ with β ≥ 2.0 maintains the initial linear waveform and reduces the thickness of the wave from its initial value. In contrast, MUSCL-THINC-pρ with β = 1.6 and MUSCL exhibit dissipative solutions, with the leading and trailing edges of the compression wave becoming a smooth profile. At t = 1520 [Fig. 7(b)], the compression wave reaches approximately x = 1685 (> x_s), where the formation of a shock wave is expected. The results for β ≥ 2.4 appear to display a discontinuous shock front. Conversely, β = 1.6, 2.0, and MUSCL exhibit dissipative profiles. None of these methods account for the coalescence of weakly compressed waves. In particular, MUSCL exhibits excessive dissipation in this problem. Consequently, it is evident that the use of THINC is necessary for accurately solving the formation of weak shock waves.

FIG. 7.

View large Download slide

Density profile in shock formation problem: (a) t = 500 and (b) t = 1520.

Next, we determine the dependence of β by evaluating the thickness of the compression or shock wave over time variation. The compression or shock wave thickness δ can be defined by using following equation:²⁹

δ = \frac{ρ_{L} - ρ_{R}}{{(d ρ / d x)}_{max}} = \frac{ρ_{L} - ρ_{R}}{{max}_{i}^{} [(ρ_{i} - ρ_{i + 1}) / Δ x]} .

(22)

Because we use the cell size Δx of unity in this problem, the thickness δ corresponds to the number of cells. Also, we confirmed that the number of cells expressing the shock wave is consistent across coarse and fine grids. Thus, the following discussion using the number of internal points of the waves is applicable regardless of the grid resolution.

A plot depicting δ against wave position is presented in Fig. 8. δ was calculated and plotted every ten steps. The wave position is defined as the location of the maximum slope of the wave. The red solid line represents the theoretical value obtained from analysis using the characteristic velocity. It is noteworthy that while δ theoretically reaches zero at the shock formation distance, δ numerically reaches a certain thickness due to the appearance of the internal structure of the shock wave.

FIG. 8.

View large Download slide

Wave thickness and position history in weak shock formation problem.

It was indeed challenging to judge the validity of the final δ value from both physical and numerical viewpoints. From a physical perspective, the shock wave is considered a nearly discontinuous surface across which variables exhibit jumps. From a numerical perspective, on the contrary, a FVM inevitably requires at least one cell (typically a few cells) to express a discontinuity, unless the discontinuity precisely aligns with the gridline.^33,34 In this study, we reconciled these differing viewpoints by defining a “discontinuous shock wave” as “a thin zone (δ < 5 in this paper) established nearly at the theoretical shock formation instant” (in this numerical example, at t = 1520). This criterion was consistently applied throughout all the numerical tests.

The numerical results clearly indicate the existence of a physically consistent value of β. For β = 1.6, the compression wave thickness did not settle at a certain value. At x = 2000, the thickness was 8.6 cells, indicating that β = 1.6 leads to excessive dissipation for shock formation with a shock Mach number of 1.01. With β = 2.0, δ stabilized at 5.7 cells after the wave passed through the expected shock formation point. β = 2.4 and 2.8 showed acceptable solutions: the thickness of the compression waves decreased to a stable, thin value (4.4 and 3.6 cells, respectively). The point at which δ reached a stable value was close to the theoretical value for the shock formation distance. Consequently, under these conditions, the shock wave formation was accurately resolved. β = 3.2 exhibited peculiar behavior: the thickness of the compression waves decreased more rapidly than the theoretical trend, and a shock was formed at approximately x = 1250, much closer to the origin than the theoretical distance. This indicates that the compression wave underwent stronger sharpening than the physically accurate compression under this β.

From the shock formation problem, we recommend β = 2.4 and 2.8 for MUSCL-THINC-pρ (the proposed scheme in this paper) and MUSCL-THINC (the original scheme²⁰). However, as will be demonstrated later in the Appendix, MUSCL-THINC with β = 2.8 exhibits an anomalous shape of vortex in the two-dimensional viscous shock-tube problem (which is not observed in MUSCL-THINC-pρ). To avoid any confusion, we adopt β = 2.4 as the recommended value hereafter.

B. Isotropic vortex advection (unsteady, accuracy study)

The advection of an isotropic vortex³⁵ is solved to verify the accuracy of the proposed method. This is a widely used benchmark test for the Euler equations. The exact solution of an isentropic vortex is given by

\begin{matrix} (\begin{matrix} u \\ v \end{matrix}) = \frac{σ}{2 π r_{c}} exp [1 - \frac{{(r / r_{c})}^{2}}{2}] (\begin{matrix} - y \\ x \end{matrix}), \\ p = ρ T, \\ T = 1 - \frac{(γ - 1) σ^{2}}{8 γ π^{2}} exp [1 - {(\frac{r}{r_{c}})}^{2}] . \end{matrix}

(23)

Here, r represents the distance from the center of the computational domain in this section, while

σ

and r_c are constants that determine the characteristics of the vortex. In this case, we used

σ

= 5 and r_c = 1. The vortex was initially placed at the origin within the square two-dimensional computational domain [−5, 5] × [−5, 5]. The mean flow (u, v) is represented as (u, v) = (1, 1).

To quantitatively assess the accuracy of the numerical scheme for a continuous flow, we compared the numerical results with the exact solution. We obtained the exact solution by applying a Galilean transformation to the initial conditions. Theoretically, when (u, v) = (1, 1), the vortex reaches (x, y) = (2, 2) at t = 2. We quantitatively evaluated the differences between the theoretical and numerical solutions using both the

L_{1} - norm

and the

L_{\infty} - norm

⁠. Each norm is defined as follows:

\begin{matrix} (L_{1} - norm) = \frac{Σ_{i, j} (ρ_{i, j} - ρ_{i, j exact})}{i_{max} \times j_{max}}, \\ (L_{\infty} - norm) = {max}_{i, j}^{} |ρ_{i, j} - ρ_{i, j exact}| . \end{matrix}

(24)

The computational domain was divided equally into grids of size 20 × 20, 40 × 40, 80 × 80, 160 × 160, and 320 × 320. The boundary conditions were set as periodic. The time increments varied from Δt = 8 × 10⁻² (for 20 × 20 grids) to 5 × 10⁻³ (for 320 × 320 grids), while ensuring that the CFL number remained at approximately 0.5.

Graphs illustrating the $L_{1} -$ and $L_{\infty} - norms$ plotted against the mesh size are presented in Fig. 9. The results for MUSCL, MUSCL-THINC, and MUSCL-THINC-pρ are displayed. All methods demonstrated second-order accuracy in this problem. This is because all the methods employ minmod-limited MUSCL for continuous flow, which provides second-order spatial accuracy with linear reconstruction within the cells. MUSCL-THINC and MUSCL-THINC-pρ exhibited $L_{\infty} - norm$ approximately 0.5 orders of magnitude smaller than MUSCL. Consequently, it is found that using THINC reconstruction yielded higher resolutions for continuous flow. MUSCL-THINC showed smaller errors than MUSCL-THINC-pρ on the coarse grid. The former employs THINC depending on the slope ratio within the stencil, and the latter restricts usage of THINC depending on nonlinearity in addition to the slope ratio. Thus, it is implied that the aggressive usage of THINC for continuous flows seems to enhance the accuracy on low-resolution grids. In finer grids, there was not a significant difference in accuracy between MUSCL-THINC and MUSCL-THINC-pρ.

FIG. 9.

Errors of isotropic vortex advection problem: (a) L 1 − norm and (b) L ∞ − norm.

View large Download slide

Errors of isotropic vortex advection problem: (a) $L_{1} - norm$ and (b) $L_{\infty} - norm$ ⁠.

C. Sod shock tube (unsteady, accuracy study)

The next step involves verifying the accuracy of computations for shock waves. To achieve this, we solved Sod's shock-tube problem,³⁶ which features a strong shock wave and a moving density discontinuity. In this study, the computational domain was one-dimensional, ranging from [0, 1], and was divided into 200 equally spaced cells (Δx = 0.005). The diaphragm, initially positioned at the center of the computational domain, separates the driver and driven gases. The initial conditions on the left and right sides of the diaphragm are defined as follows:

(ρ_{L}, u_{L}, p_{L}) = (1.0, 0, 1.0), (ρ_{R}, u_{R}, p_{R}) = (0.125, 0, 0.1) .

(25)

The diaphragm was removed instantaneously at the initiation of the computation. The time increment Δt was 1 × 10⁻³ (CFL number ≈ 0.45). We also prepared the reference solution solved in 2000 cells by MUSCL with the time increment of 1 × 10⁻⁴.

The density profile at t = 0.2 (after 200 steps) is presented in Fig. 10. The top-right panel provides an enlarged view of the moving density discontinuity (contact surface), while the bottom-right panel focuses on the region around the shock wave. Additionally, the weight value of THINC (⁠ $ζ_{i}$ for MUSCL-THINC, and $ζ_{i} ξ_{i}$ in the case of MUSCL-THINC-pρ) is displayed in the same figure. The density profiles generated by MUSCL-THINC and MUSCL-THINC-pρ were nearly identical, with a slight but noteworthy difference. The shock wave (0.84 < x < 0.86) thickness for MUSCL-THINC-pρ was slightly broader than that of MUSCL-THINC, and the waveform closely resembled that of MUSCL. This variation can be attributed to the reduction in the influence of THINC around the shock wave in MUSCL-THINC-pρ, as evidenced by the smaller value of $ζ_{i} ξ_{i}$ compared to $ζ_{i}$ ⁠, as shown in the weight value plot. This slight broadening of the strong shock wave is intended and acceptable because the shock wave is still adequately represented within a few cells. This outcome confirms that THINC reconstruction is unnecessary for strong shock waves. Conversely, at the moving density discontinuity (0.67 < x< 0.7), MUSCL-THINC-pρ accurately captured it, similar to the performance of MUSCL-THINC. The weight of THINC was similar in both cases. Therefore, in line with the intentions outlined in Sec. II B 4, the proposed scheme utilizes MUSCL for strong shock waves and THINC for density discontinuities. As discussed later, this feature contributes to the stability of steady-state flow simulations.

FIG. 10.

View large Download slide

Density profile of Sod's shock tube. Right column shows enlarged view of density and weight of THINC around moving density discontinuity (top) and shock wave (bottom).

Before concluding this section, we discuss the unnecessary activation of THINC behind the shock wave (0.82 < x < 0.84), in front of it (0.86 < x < 0.89), and behind the density discontinuity (0.63 < x < 0.67). Although eliminating this unnecessary activation is challenging, it is crucial to note that the oscillation level in this study was negligible. We also confirmed that such activations did not appear within the smooth region (the expansion wave region).

D. Moving shock (unsteady)

In this section, we present the computational results for a moving shock wave with M_s = 1.01 (very weak), 1.5 (moderate), and 3.0 (strong), where M_s represents the shock Mach number. The computations were conducted in a domain spanning [0, 2000], with 2000 cells evenly distributed in the x-direction (Δx = 1). The initial placement of the moving shock wave was at the interface between the 100th and 101st cells. We focus solely on one-dimensional results here, as our preliminary survey did not reveal any two-dimensional shock anomalies, e.g., the carbuncle phenomenon, in the corresponding test with the computational domain extended in the y-direction of the Cartesian coordinate. The initial conditions are defined by the Rankine–Hugoniot equations for a moving shock wave, given as follows:

\begin{matrix} (ρ_{L}, u_{L}, p_{L}) = (\frac{(γ + 1) M_{s}^{2}}{(γ - 1) M_{s}^{2} + 2} ρ_{R}, \frac{2 a_{R}}{γ + 1} (M_{s} - \frac{1}{M_{s}}) + u_{R}, [1 + \frac{2 γ}{γ + 1} (M_{s}^{2} - 1)] p_{R}), \\ (ρ_{R}, u_{R}, p_{R}) = (1, 0, 1) . \end{matrix}

(26)

The time increment Δt was 0.5, 0.4, and 0.2 (CFL number at the initial state ≈ 0.51, 0.74, and 0.77) for M_s = 1.01, 1.5, and 3.0 respectively.

Figure 11 displays the density profile, computed using MUSCL, MUSCL-THINC, and MUSCL-THINC-pρ. The results at time step 2000 are shown: t = 1000, 800, and 400 for M_s = 1.01, 1.5, and 3.0 respectively. MUSCL failed to sharply resolve the shock wave (M_s = 1.01). In contrast to Sec. III A, the initial conditions were perfectly discontinuous in this case. Thus, MUSCL failed to solve the very weak shock wave, regardless of the initial profile. In contrast, both MUSCL-THINC and MUSCL-THINC-pρ sharply captured the very weak shock wave, signifying a considerable improvement over MUSCL. This outcome underscores the necessity of nonlinear reconstruction via THINC for accurately capturing a very weak, moving shock wave.

FIG. 11.

View large Download slide

Density profile of moving shock wave at time step of 2000: (a) M_s = 1.01, (b) 1.5, and (c) 3.0.

For M_s = 1.5 and 3.0, MUSCL showed a discontinuous profile. MUSCL-THINC and MUSCL-THINC-pρ also showed a discontinuous profile. The only difference was the number of cells expressing the shock wave. MUSCL-THINC solved the shock wave with two cells for M_s = 1.5 and 3.0, whereas MUSCL-THINC-pρ solved it with four cells (M_s = 1.5) or three cells (M_s = 3.0). As discussed in Sod's shock-tube problem (Sec. III C), suppressing THINC at the shock wave using the nonlinearity-weighted parameter resulted in a slight increase in the shock wave thickness.

E. Stationary shock (steady)

This section presents the computational results for the stationary normal shock wave standing in the supersonic flow. The Mach number of the upstream flow M was 1.01, 1.5, and 3.0. The computation was performed in a one-dimensional space of [0, 100], where 100 cells were arrayed in the x-direction (Δx = 1). A normal shock wave was installed at the interface between the 50th and 51st cells in the x-direction. The left and right sides of the shock wave were supersonic and subsonic, respectively. A supersonic inlet was applied to the left boundary (x = 0). The wall-like boundary used in Ref. 24 was applied to the right-side boundary (x = 100); only the mass flux was fixed (ρu = 1) to stabilize the shock wave position, and the other variables were simply extrapolated (q_i_max+1 = q_i_max). For a detailed description of this problem, refer to Ref. 24.

Figure 12 depicts the density profiles obtained using MUSCL, MUSCL-THINC, and MUSCL-THINC-pρ at 100 000 time steps. It is evident that MUSCL failed to sharply resolve the M = 1.01 stationary shock wave due to significant numerical dissipation. In contrast, both MUSCL-THINC and MUSCL-THINC-pρ captured the shock wave sharply. Although a slight post-shock oscillation was observed, the implementation of THINC significantly improved the resolution for a very weak shock wave.

FIG. 12.

View large Download slide

Density profile of stationary shock wave at time step of 100 000: (a) M = 1.01, (b) 1.5, and (c) 3.0.

MUSCL-THINC exhibited periodic post-shock oscillations at M = 1.5 and 3.0. In contrast, MUSCL-THINC-pρ yielded less oscillatory results due to the improved hybrid balance between MUSCL and THINC achieved through the nonlinearity-weighted parameter. MUSCL provided sharp and stable solutions for the M = 1.5 and 3.0 stationary shock waves, indicating that THINC reconstruction is unsuitable for moderate and strong stationary shock computations, a scenario not previously explored in studies such as Ref. 20. Hence, it was verified that reconstruction must be applied based on nonlinearity, primarily determined by the shock wave strength. The proposed scheme effectively adjusts the balance between MUSCL and THINC, resulting in a high-resolution and stable solution. The qualitative convergence performance is evaluated using the 2D blunt-body problem described in Sec. III F.

F. Blunt body problem (steady)

Convergence performance to a steady state has rarely been discussed in the context of developing numerical schemes using THINC reconstruction, despite its importance in application problems. In this section, to evaluate the convergence performance of the proposed scheme for a stationary problem, a two-dimensional bow shock formed in front of a blunt body (so-called “blunt body problem”) was solved. This problem serves as a benchmark to evaluate the performance of computational schemes for strong steady shock waves.^8,24 The focus here is clarifying the effect of THINC usage at strong shock waves on the convergence performance. The uniform-flow Mach number was set to M = 3.0. We utilized the grid used by Jiang and Shu.⁸ The computational plane can be expressed as follows:

\begin{matrix} x = - [R_{x} - (R_{x} - 1) i] \cos [θ (2 j - 1)], \\ y = [R_{y} - (R_{y} - 1) i] \sin [θ (2 j - 1)], \end{matrix}

(27)

where

R_{x}

= 3,

R_{y}

= 6, θ = 5π/12, and the number of cells was 60 in the wall-normal direction (i-direction) and 80 in the tangential direction (j-direction). A supersonic flow inlet (M = 3.0) was applied to the left boundary (i = 0), and a slip wall condition was applied to the right boundary (i = 60). Simple extrapolation was used for the edges of the tangential wall directions (j = 0 and 80). In this grid setup, the bow shock wave did not align perfectly with a grid line but was expected to converge to a stationary state as the calculation progressed. We used a time increment with a CFL number of 0.5.

Figure 13 displays the computational grid [Fig. 13(a)] and density contours at 100 000 time steps for MUSCL [Fig. 13(b)], MUSCL-THINC [Fig. 13(c)], and MUSCL-THINC-pρ [Fig. 13(d)]. MUSCL [Fig. 13(b)] and MUSCL-THINC-pρ [Fig. 13(d)] successfully resolved the bow shock and the downstream flow field without apparent anomalies. However, MUSCL-THINC [Fig. 13(c)] exhibited wrinkled lines behind the shock wave compared to the other cases, roughly around coordinates (x, y) = (−1, ±2). The cause of this anomaly can be understood by examining the THINC weight values. The weight values of THINC (⁠ $ζ_{i}$ for MUSCL-THINC and $ζ_{i} ξ_{i}$ for MUSCL-THINC-pρ) used in the density reconstruction in the i-direction, along with density contour lines (in black), are illustrated in Fig. 14. The key difference between MUSCL-THINC and MUSCL-THINC-pρ in these figures is the weight value of THINC around the shock wave. In MUSCL-THINC [Fig. 14(a)], the weight value of THINC $ζ_{i}$ is close to 1 within and behind the shock wave. In contrast, the weight value of THINC in MUSCL-THINC-pρ $ζ_{i} ξ_{i}$ is suppressed within and behind the shock wave, as shown in Fig. 14(b). (Although a large THINC weight value was applied to two or three cells in front of the shock, this was due to a slight density increase. This phenomenon was observed in both MUSCL-THINC and MUSCL-THINC-pρ.) As observed in previous numerical tests of one-dimensional stationary shock waves (Sec. III E), using THINC leads to numerical oscillations behind strong shocks. Therefore, this test reveals that the weight value of MUSCL-THINC for strong shock waves was inappropriate for this steady, two-dimensional problem.

FIG. 13.

View large Download slide

Two-dimensional blunt body problem (M = 3.0): (a) Computational grid including ghost cells, (b)–(d) density contours ranging from 1.2 to 4.2 with 16 levels at the time step is 100 000. (b) MUSCL, (c) MUSCL-THINC, and (d) MUSCL-THINC-pρ.

FIG. 14.

View large Download slide

Contours showing weight value of THINC in i-direction and density (lines) at the time step of 100 000: (a) MUSCL-THINC and (b) MUSCL-THINC-pρ.

Figure 15 shows the density residuals, which represent the maximum value of the density changes between every time step in all computational domains. In MUSCL-THINC, the residual reached 1 × 10⁻² at a time step of 5000 and then did not decrease further, indicating poor convergence performance for steady computation with a strong shock. Conversely, for MUSCL-THINC-pρ, the density residual converged to 1 × 10⁻¹³ at a time step of approximately 15 000. Compared to MUSCL-THINC, the residual decreased by more than 10 orders of magnitude. As shown in Fig. 14, which illustrates the weight value of THINC, the use of THINC in MUSCL-THINC-pρ was suppressed near the strong shock wave compared to MUSCL-THINC. This suppression contributed to stabilizing the flow field behind the shock wave, leading to convergence to a steady state. MUSCL showed the best convergence performance, with the residual reaching 1 × 10⁻¹⁵. Although the residual of MUSCL-THINC-pρ was approximately two orders larger than that of MUSCL, MUSCL-THINC-pρ exhibited considerably improved convergence over MUSCL-THINC.

FIG. 15.

View large Download slide

Density residual in two-dimensional blunt body problem.

In summary, this test revealed that the proposed hybridization strategy successfully addressed steady-state problems involving strong shock waves. This was achieved owing to successful tuning of the THINC utilization at strong shock waves. The proposed approach, MUSCL-THINC-pρ, represents the first hybrid scheme combining MUSCL and THINC with significantly improved convergence performance for steady solutions of compressible flows with strong shock waves maintaining a high-resolution feature for weak shock waves.

G. Shu–Osher's shock-entropy wave interaction (unsteady)

Here, we solved the Shu–Osher problem,³⁷ where the interaction between shock waves and entropy waves occurs in one dimension. The computational domain was [−5, 5], discretized into N = 400 evenly spaced cells. The initial conditions were specified as follows:

\begin{matrix} (ρ_{L}, u_{L}, p_{L}) = (3.857 143, 2.629 369, 10.333 33) for x \leq - 4, \\ (ρ_{R}, u_{R}, p_{R}) = (1 + 0.2 sin 5 x, 0, 1) for x > - 4 . \end{matrix}

(28)

This configuration results in a moving shock wave with M_s = 3.0 when sin5x = 0. We utilized a time increment of 3.6 × 10⁻³ (yielding a CFL number of approximately 0.68) and conducted computations up to t = 1.8 (500 steps).

Figure 16 presents the density profiles at t = 1.8. For comparison, we also include the reference result obtained using MUSCL on a mesh four times finer (1600 cells) with a time step of 9 × 10⁻⁴. MUSCL exhibited significant numerical dissipation in the entropy waves after their interaction with the shock. In contrast, both MUSCL-THINC and MUSCL-THINC-pρ resolved the entropy waves after the interaction, comparable to the reference result with four times the number of cells. The difference between MUSCL-THINC and MUSCL-THINC-pρ was minimal in this problem, confirming that the nonlinearity-weighted parameter does not obviously compromise the resolution of entropy waves, as intended in Sec. II B 4.

FIG. 16.

View large Download slide

density profile of Shu–Osher's shock/entropy wave interaction problem (t = 1.8): (a) −3 ≤ x ≤ 5 and (b) −1 ≤ x ≤ 3.

H. 2D Riemann problem (unsteady)

In this section, we address the two-dimensional Riemann problem, a widely recognized benchmark test for the Euler equations.³⁸ In this study, we initiated density and tangential velocity discontinuities. As time progressed, the Kelvin–Helmholtz (KH) instability emerged, particularly if the scheme exhibited minimal numerical dissipation. The scheme's resolution can be qualitatively assessed by examining the appearance of small-scale flow structures of the KH instability. The computational domain was defined as [0, 1] × [0, 1], equally divided into 1000 cells in both the x- and y-directions. Initially, the square computational domain was divided into smaller four-square domains as follows:

(ρ_{0}, u_{0}, v_{0}, p_{0}) = \{\begin{matrix} \begin{matrix} (1.0, 0.726, 0.0, 1.0) & x \leq 0.5, y \geq 0.5, \\ (0.8, 0.0, 0.0, 1.0) & x \leq 0.5, y < 0.5, \end{matrix} \\ \begin{matrix} (0.5313, 0.0, 0.0, 0.4) & x > 0.5, y \geq 0.5, \\ (1.0, 0.0, 0.7276, 1.0) & x > 0.5, y < 0.5 . \end{matrix} \end{matrix}

(29)

Time integration was performed using CFL = 0.5.

The density contours obtained by each method at a time step of 2500 are presented in Fig. 17. In Fig. 17(a), MUSCL demonstrates a dissipative density interface within the low-density region where x < 0.5 and y < 0.5. MUSCL struggles to resolve the KH instability due to significant numerical dissipation. In contrast, both MUSCL-THINC [Fig. 17(b)] and MUSCL-THINC-pρ [Fig. 17(c)] exhibit the KH instability at the density discontinuity along the initial slip lines near (x, y) = (0.5, 0.4) and (0.4, 0.5). This is attributed to the THINC reconstruction, which captures the density discontinuity with reduced numerical dissipation. Additionally, instability is observed in the high-density region with x < 0.5 and y < 0.5. Figure 17(d) presents the results obtained using WENO⁸ reconstruction for comparison with higher-order schemes. Reconstruction was carried out in a fifth-order component-wise manner. Notably, MUSCL-THINC-pρ exhibits superior resolution for density discontinuity compared to the WENO reconstruction. This is because, for the discontinuity, the higher-order scheme is not necessarily optimized: the higher-order result is obtained for continuous flow components. The interface capturing by THINC contributed to capturing the density discontinuity with less numerical dissipation. Consequently, despite operating within a second-order framework, MUSCL-THINC-pρ achieves resolution surpassing that of the approximately 26% more expensive higher-order scheme for capturing slip lines, which is computationally more demanding due to its reliance on information from a larger number of cells. The results emphasize the effectiveness of the proposed method (see Sec. III K for the specific numerical costs). This feature is highly advantageous for practical simulations.

FIG. 17.

View large Download slide

Density contour of two-dimensional Riemann problem (2500 steps) solved by (a) MUSCL, (b) MUSCL-THINC, (c) MUSCL-THINC-pρ, and (d) WENO.

I. Double-Mach reflection (unsteady)

In this section, a hypersonic planar shock wave with a shock Mach number M_s = 10 reflected by a 30° ramp is simulated.³⁹ This problem, known as the double-Mach reflection problem, allows for the evaluation of the scheme's performance with regard to strong shock waves and density discontinuities in the recirculation region behind the Mach stem. Chiu et al.²⁰ previously tested the performance of the original hybrid MUSCL-THINC scheme using this problem and reported that the scheme could resolve the small structure of the KH instability in the recirculation region. They also reported that MUSCL-THINC suffered slightly from a carbuncle phenomenon behind the Mach stem. In this study, the proposed scheme, which switches to MUSCL around a strong shock wave, was compared with MUSCL-THINC.

The computation was performed in the [0, 4] × [0, 1] domain and spaced at 1000 and 250 cells in the x- and y-directions, respectively. The initial conditions were as follows:

(ρ, u, v, p) = \{\begin{matrix} (8.0, 7.145, - 4.125, 116.5) \\ (1.4, 0, 0, 1) \end{matrix} \begin{matrix} x < x_{s} (t = 0), \\ otherwise, \end{matrix}

(30)

where

x_{s} (t = 0)

was the initial shock wave position,

x_{s} (t = 0) = \frac{1}{6} + \frac{10}{sin (α)} + \frac{y}{tan (α)},

(31)

where α = 60° was the angle of the shock wave from the bottom boundary. The same boundary conditions as those used in ²⁰ were employed. The time increment was set to 2 × 10⁻⁵ (CFL number ≈ 6.4 × 10⁻²), and the computation was performed until t = 0.2 (10 000 steps).

The density contours obtained from the simulations using MUSCL, MUSCL-THINC, and MUSCL-THINC-pρ are illustrated in Fig. 18. The enlarged views of the recirculation region and Mach stem are presented in the right panels. In the results obtaining using MUSCL [Fig. 18(a)], a dissipative solution is evident, with no vortex observed in the recirculation region. In contrast, MUSCL-THINC [Fig. 18(b)] exhibited the generation of vortices in the recirculation region, consistent with previous findings.²⁰ However, the Mach stem displayed unnatural wrinkled density contour lines extending from (x, y) = (2.8, 0) to (2.74, 0.4), accompanied by a two-dimensional instability behind the shock wave in the vicinity of the bottom wall (around 2.7 < x < 2.8), indicative of an anomalous shock solution. (This anomaly was previously termed the “carbuncle phenomenon” in Ref. 20.) In the case of MUSCL-THINC-pρ [Fig. 18(c)], the density contour of the Mach stem appears straighter without the shock anomaly. In addition, the flow field behind the shock wave exhibited fewer numerical oscillations, while successfully resolving the KH instability in the recirculation region. Thus, the proposed scheme mitigates the shock anomaly by considering the strength of the shock wave while maintaining favorable performance in capturing density discontinuities. This improvement in robustness against strong shock waves is crucial for the practical application of numerical methods, as shock anomalies frequently arise in real-world scenarios.^21–24

FIG. 18.

View large Download slide

Density contours of double-Mach reflection problem at t = 0.2: (a) MUSCL, (b) MUSCL-THINC, and (c) MUSCL-THINC-pρ.

J. Viscous shock-tube (unsteady, viscous)

For the final numerical test, we tackled a two-dimensional shock wave-boundary layer interaction problem within a shock tube (viscous shock-tube problem).⁴⁰ This problem involved solving the Navier–Stokes equations to assess the performance of the scheme in handling viscous effects. Viscous flux computations were carried out using a second-order central difference. We adopted the setup described by Daru et al.,⁴⁰ with the computational domain spanning [0, 1] × [0, 0.5], divided equally into 250 × 125 cells. The initial condition was specified as

\begin{matrix} (ρ_{L}, u_{L}, v_{L}, p_{L}) = (120, 0, 0, 120 / γ) for x \leq 0.5, \\ (ρ_{R}, u_{R}, v_{R}, p_{R}) = (1.2, 0, 0, 1.2 / γ) for x > 0.5 . \end{matrix}

(32)

After evolution, a shock wave propagated toward the right boundary, leading to the emergence of a boundary layer behind it due to viscosity. Upon reflection at the right boundary, an interaction occurred between the reflected shock wave and the boundary layer. Our focus was specifically on evaluating the scheme's capacity to resolve the primary vortex formed after this interaction. The computations were conducted using a time step of 5 × 10⁻⁶ (yielding a CFL number of approximately 0.52) for t = 1 (corresponding to 200 000 time steps), with the Reynolds number set at 200.

The density contours at t = 1 are depicted in Fig. 19. Additionally, we included a reference solution obtained with a grid resolution of 1000 × 500 cells (i.e., 4 × 4 times finer than the base resolution) using MUSCL. High-resolution schemes are known to exhibit longitudinal elongation of the primary vortex at approximately (x, y)= (0.8, 0.05).⁴⁰ The dissipative nature of the MUSCL scheme is evident from the horizontally lengthened vortex shape in Fig. 19(a). In contrast, both MUSCL-THINC [Fig. 19(b)] and MUSCL-THINC-pρ [Fig. 19(c)], as well as the reference solution [Fig. 19(d)], demonstrated characteristics of a high-order scheme, particularly the longitudinal elongation of the primary vortex. Remarkably, the results obtained with MUSCL-THINC-pρ closely resembled those of the reference solution. Furthermore, we observed that the MUSCL-THINC-pρ solution exhibited fewer numerical oscillations compared to MUSCL-THINC. The density profiles along the bottom wall are displayed in Fig. 20. We confirmed that the MUSCL-THINC-pρ solution quantitatively resembled the profile of the reference solution around the primary vortex (x ≈ 0.8). Consequently, our findings confirm that the proposed scheme yields favorable results for handling viscous problems.

FIG. 19.

View large Download slide

Density contours of viscous shock-tube problem at t = 1. Density ranges from 20 to 140 with 41 levels: (a) MUSCL, (b) MUSCL-THINC, (c) MUSCL-THINC-pρ, and (d) reference.

FIG. 20.

View large Download slide

Density profiles of viscous shock tube problem along bottom wall at t = 1.

K. Computational costs

Finally, we discuss the computational costs associated with the proposed method. The CPU time measurement results are summarized in Table I. We conducted measurements for two types of numerical tests: Shu–Osher shock wave-entropy wave interactions (Sec. III G) and the two-dimensional Riemann problem (Sec. III H). The entire process, from the start of computation to completion, was included in the measurement. The CPU time of MUSCL was used as the baseline, and the ratio to MUSCL is presented in parentheses.

TABLE I.

Comparison of CPU time.

	MUSCL	MUSCL-THINC (original)	MUSCL-THINC-pρ (proposed)	MUSCL (four times fine grid)	WENO
Shu–Osher's problem (Sec. III G)	0.300 (1)	0.341 (1.137)	0.357 (1.291)	4.910 (16.36)	⋯
2D Riemann problem (Sec. III H)	2563 (1)	3096 (1.207)	3370 (1.315)	⋯	4256 (1.660)

	MUSCL	MUSCL-THINC (original)	MUSCL-THINC-pρ (proposed)	MUSCL (four times fine grid)	WENO
Shu–Osher's problem (Sec. III G)	0.300 (1)	0.341 (1.137)	0.357 (1.291)	4.910 (16.36)	⋯
2D Riemann problem (Sec. III H)	2563 (1)	3096 (1.207)	3370 (1.315)	⋯	4256 (1.660)

In the case of the Shu–Osher's problem, we observed a 13.7% increase in computational time for MUSCL-THINC and a 29.1% increase for MUSCL-THINC-pρ compared to MUSCL under the same grid conditions. It is important to note that MUSCL provided a suboptimal solution due to excessive dissipation, and to achieve similar high-resolution results with MUSCL, a significantly finer grid and a fourfold reduction in the time increment were required, resulting in an approximately 16-fold increase in computational time.

In the context of the two-dimensional Riemann problem, MUSCL-THINC exhibited a 20.7% increase in computational time, while MUSCL-THINC-pρ showed a 31.5% increase compared to MUSCL. Once again, MUSCL delivered a subpar solution for this problem. Additionally, it is worth noting that the proposed scheme, despite the increase in computational time, proved to be more efficient than WENO, which exhibited a 66.0% increase from MUSCL and presented a more dissipative solution than MUSCL-THINC and MUSCL-THINC-pρ at the density discontinuity.

In summary, MUSCL-THINC-pρ resulted in a computational time increase of 29.1%–31.5% compared to MUSCL and 8.87%–13.7% compared to the original MUSCL-THINC. Considering the enhanced convergence performance and robustness against strong shock waves, this increase in computational time remains acceptable. This represents a significant advantage of the proposed method for practical application problems.

IV. CONCLUSIONS

This paper introduced a novel hybrid scheme of MUSCL and THINC, termed the “MUSCL-THINC-pρ” scheme, or briefly, T-MUSCL, for both steady and unsteady single-phase compressible flow simulations. The design of this scheme revolved around achieving an appropriate balance of nonlinearity between the physical phenomena and the reconstruction process to solve weak shock waves sharply and strong shock waves robustly. This concept was realized through the utilization of two key parameters: a nonlinearity-weighted parameter and a slope-ratio-weighted parameter. The proposed scheme offers the following distinctive features.

Enhanced accuracy for continuous flow simulations, providing second-order spatial accuracy with smaller errors compared to conventional MUSCL.
Precise capturing of extremely weak moving and stationary shock waves (M_s = 1.01 and M = 1.01, respectively). These cases were previously challenging for the original MUSCL method due to excessive dissipation.
Improved convergence behavior in steady two-dimensional blunt body problems. The density residual decreased significantly, surpassing ten orders of magnitude reduction when compared to the original MUSCL-THINC hybrid scheme.
Mitigation of shock anomalies in the presence of strong shock waves while maintaining high-resolution computations for density discontinuities and entropy waves. Remarkably, the proposed scheme achieved results with resolution exceeding that of a more computationally intensive fifth-order reconstruction at slip lines in the two-dimensional Riemann problem, despite being designed within a second-order framework.

It is noteworthy that using THINC reconstruction for strong shock waves resulted in undesirable numerical oscillations and convergence issues, particularly in steady flows with intense shock waves. This finding contributes valuable insights for the continued development of numerical schemes in this field.

Despite these advanced features, the algorithm employed in this scheme remains straightforward: after computing the nonlinearity-weighted parameter based on pressure and density, it is multiplied by the slope ratio-weighted parameter used in the original hybrid MUSCL-THINC scheme. Additionally, the stencil size for reconstruction aligns with that of the second-order spatial accuracy MUSCL approach. The proposed scheme has some drawbacks: higher computational costs than that of MUSCL and the original hybrid approach, slightly insufficient small-scale resolution, and somewhat compromised convergence properties than MUSCL. Nevertheless, the proposed scheme consistently delivers robust, high-resolution computational results in a wide variety of benchmark tests. Consequently, it can be readily adopted in place of MUSCL in various practical applications without introducing unnecessary complexity to existing algorithms. The potential extension of this scheme to unstructured grid systems, which could enhance its applicability in more real-world scenarios, remains a promising avenue for future research.

ACKNOWLEDGMENTS

The authors express their sincere gratitude to Akihiro Sasoh at Nagoya University for engaging in fruitful discussions regarding the proposed scheme and providing valuable insights into weak-shock physics. Special thanks are also extended to Junya Aono, Yoshikatsu Furusawa, and Tomohiro Mamashita of Yokohama National University for their helpful comments. This research was supported by Kakenhi Grant Nos. JP21H04589, JP23H01601, and JP23KJ0981. We would like to thank Editage (www.editage.jp) for English language editing.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Gaku Fukushima: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Resources (equal); Software (equal); Validation (lead); Visualization (lead); Writing – original draft (lead). Keiichi Kitamura: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (equal); Software (equal); Validation (supporting); Writing – original draft (supporting); Writing – review & editing (lead).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX: EFFECT OF β IN VISCOUS SHOCK-TUBE PROBLEM

As discussed in the one-dimensional weak shock formation problem (Sec. III A), β values of 2.4 and 2.8 were found to yield physically consistent compression trends for MUSCL-THINC and MUSCL-THINC-pρ. In this Appendix, we further assess the impact of β on the viscous shock-tube problem (Sec. III J). Figure 21 displays the density contours of the viscous shock-tube problem, solved using MUSCL-THINC with β = 2.4 [Fig. 21(a)] and 2.8 [Fig. 21(b)], MUSCL-THINC-pρ with β = 2.4 [Fig. 21(c)], and β = 2.8 [Fig. 21(d)]. Notably, when β was set to 2.8 for MUSCL-THINC, it resulted in an anomalous deformation of the primary vortex, rendering it excessively thin compared to the other cases. Moreover, numerical oscillations were evident throughout the region behind the shock waves. This behavior stemmed from the strong sharpening effect associated with the larger value of β. Conversely, the results for the other cases did not exhibit anomalous behavior and closely resembled the reference outcomes [Fig. 21(e)]. While MUSCL-THINC-pρ with β = 2.8 did not produce anomalous results, it is recommended to opt for β = 2.4 to simplify the user experience, as this value can be effectively used for MUSCL-THINC computations.

FIG. 21.

View large Download slide

Density contours of viscous shock tube problem at t = 1. Density ranges from 20 to 140 with 41 levels: (a) MUSCL-THINC β = 2.4, (b) MUSCL-THINC β = 2.8, (c) MUSCL-THINC-pρ β = 2.4, (d) MUSCL-THINC-pρ β = 2.8, and (e) reference.

REFERENCES

1.

B.

van Leer

, “

Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method

,”

J. Comput. Phys.

32

,

101

–

136

(

1979

).

https://doi.org/10.1016/0021-9991(79)90145-1

Google Scholar

Crossref

2.

W. K.

Anderson

,

J. L.

Thomas

, and

B.

Van Leer

, “

Comparison of finite volume flux vector splittings for the Euler equations

,”

AIAA J.

24

,

1453

–

1460

(

1986

).

https://doi.org/10.2514/3.9465

Google Scholar

Crossref

3.

A.

Hashimoto

et al, “

Toward the fastest unstructured CFD code ‘FaStar

,’ ” AIAA Paper No. 2012-1075,

2012

.

Google Scholar

4.

K. S.

Abdol-Hamid

,

J.-R.

Carlson

,

C. L.

Rumsey

,

E. M.

Lee-Rausch

, and

M. A.

Park

, “

Sixth drag prediction workshop results using FUN3D with k-kL-MEAH2015 turbulence model

,”

J. Aircr.

55

,

1458

–

1468

(

2018

).

https://doi.org/10.2514/1.C034481

Google Scholar

Crossref

PubMed

5.

G.

Fukushima

,

K.

Kitamura

, and

A.

Sasoh

, “

On peculiar behaviors of captured very-weak moving shockwaves

,” AIAA Paper No. 2022-4127,

2022

.

Google Scholar

6.

T.

Nonomura

,

K.

Kitamura

, and

K.

Fujii

, “

A simple interface sharpening technique with a hyperbolic tangent function applied to compressible two-fluid modeling

,”

J. Comput. Phys.

258

,

95

–

117

(

2014

).

https://doi.org/10.1016/j.jcp.2013.10.021

Google Scholar

Crossref

7.

C.-W.

Shu

and

S.

Osher

, “

Efficient implementation of essentially non-oscillatory shock-capturing schemes

,”

J. Comput. Phys.

77

,

439

–

471

(

1988

).

https://doi.org/10.1016/0021-9991(88)90177-5

Google Scholar

Crossref

8.

G.-S.

Jiang

and

C.-W.

Shu

, “

Efficient implementation of weighted ENO schemes

,”

J. Comput. Phys.

126

,

202

–

228

(

1996

).

https://doi.org/10.1006/jcph.1996.0130

Google Scholar

Crossref

9.

R.

Borges

,

M.

Carmona

,

B.

Costa

, and

W. S.

Don

, “

An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws

,”

J. Comput. Phys.

227

,

3191

–

3211

(

2008

).

https://doi.org/10.1016/j.jcp.2007.11.038

Google Scholar

Crossref

10.

L.

Fu

,

X. Y.

Hu

, and

N. A.

Adams

, “

A family of high-order targeted ENO schemes for compressible-fluid simulations

,”

J. Comput. Phys.

305

,

333

–

359

(

2016

).

https://doi.org/10.1016/j.jcp.2015.10.037

Google Scholar

Crossref

11.

L.

Fu

, “

Review of the high-order TENO schemes for compressible gas dynamics and turbulence

,”

Arch. Comput. Methods Eng.

30

,

2493

–

2526

(

2023

).

https://doi.org/10.1007/s11831-022-09877-7

Google Scholar

Crossref

12.

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

).

https://doi.org/10.1016/j.jcp.2016.06.051

Google Scholar

Crossref

13.

X.

Deng

,

B.

Xie

,

R.

Loubère

,

Y.

Shimizu

, and

F.

Xiao

, “

Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts

,”

Comput. Fluids

171

,

1

–

14

(

2018

).

https://doi.org/10.1016/j.compfluid.2018.05.015

Google Scholar

Crossref

14.

X.

Deng

,

Y.

Shimizu

, and

F.

Xiao

, “

A fifth-order shock capturing scheme with two-stage boundary variation diminishing algorithm

,”

J. Comput. Phys.

386

,

323

–

349

(

2019

).

https://doi.org/10.1016/j.jcp.2019.02.024

Google Scholar

Crossref

15.

F.

Xiao

,

Y.

Honma

, and

T.

Kono

, “

A simple algebraic interface capturing scheme using hyperbolic tangent function

,”

Numer. Methods Fluids

48

,

1023

–

1040

(

2005

).

https://doi.org/10.1002/fld.975

Google Scholar

Crossref

16.

F.

Xiao

,

S.

Ii

, and

C.

Chen

, “

Revisit to the THINC scheme: A simple algebraic VOF algorithm

,”

J. Comput. Phys.

230

,

7086

–

7092

(

2011

).

https://doi.org/10.1016/j.jcp.2011.06.012

Google Scholar

Crossref

17.

S.

Ii

et al, “

An interface capturing method with a continuous function: The THINC method with multi-dimensional reconstruction

,”

J. Comput. Phys.

231

,

2328

–

2358

(

2012

).

https://doi.org/10.1016/j.jcp.2011.11.038

Google Scholar

Crossref

18.

K.

Kitamura

,

M.-S.

Liou

, and

C.-H.

Chang

, “

Extension and comparative study of AUSM-family schemes for compressible multiphase flow simulations

,”

Commun. Comput. Phys.

16

,

632

–

674

(

2014

).

https://doi.org/10.4208/cicp.020813.190214a

Google Scholar

Crossref

19.

X.

Deng

,

S.

Inaba

,

B.

Xie

,

K.-M.

Shyue

, and

F.

Xiao

, “

High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces

,”

J. Comput. Phys.

371

,

945

–

966

(

2018

).

https://doi.org/10.1016/j.jcp.2018.03.036

Google Scholar

Crossref

20.

T.-Y.

Chiu

,

Y.-Y.

Niu

, and

Y.-J.

Chou

, “

Accurate hybrid AUSMD type flux algorithm with generalized discontinuity sharpening reconstruction for two-fluid modeling

,”

J. Comput. Phys.

443

,

110540

(

2021

).

https://doi.org/10.1016/j.jcp.2021.110540

Google Scholar

Crossref

21.

J. J.

Quirk

, “

A contribution to the great Riemann solver debate

,”

Numer. Methods Fluids

18

,

555

–

574

(

1994

).

https://doi.org/10.1002/fld.1650180603

Google Scholar

Crossref

22.

M.

Pandolfi

and

D.

D'Ambrosio

, “

Numerical instabilities in upwind methods: Analysis and cures for the ‘carbuncle’ phenomenon

,”

J. Comput. Phys.

166

,

271

–

301

(

2001

).

https://doi.org/10.1006/jcph.2000.6652

Google Scholar

Crossref

23.

Y.

Chauvat

,

J.-M.

Moschetta

, and

J.

Gressier

, “

Shock wave numerical structure and the carbuncle phenomenon

,”

Numer. Methods Fluids

47

,

903

–

909

(

2005

).

https://doi.org/10.1002/fld.916

Google Scholar

Crossref

24.

K.

Kitamura

,

P.

Roe

, and

F.

Ismail

, “

Evaluation of Euler fluxes for hypersonic flow computations

,”

AIAA J.

47

,

44

–

53

(

2009

).

https://doi.org/10.2514/1.33735

Google Scholar

Crossref

25.

A.

Jameson

and

T.

Baker

, “

Solution of the Euler equations for complex configurations

,” AIAA Paper No. 83-1929,

1983

.

Google Scholar

26.

K.

Kitamura

and

E.

Shima

, “

Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for AUSM-family schemes

,”

J. Comput. Phys.

245

,

62

–

83

(

2013

).

https://doi.org/10.1016/j.jcp.2013.02.046

Google Scholar

Crossref

27.

C.

Hirsch

,

Numerical Computation of Internal and External Flows

, Volume 2 of Wiley Series in Numerical Methods (

John Wiley & Sons

,

1990

).

Google Scholar

28.

P. L.

Roe

, “

Characteristic-based schemes for the Euler equations

,”

Annu. Rev. Fluid Mech.

18

,

337

–

365

(

1986

).

https://doi.org/10.1146/annurev.fl.18.010186.002005

Google Scholar

Crossref

29.

A. H.

Shapiro

,

The Dynamics and Thermodynamics of Compressible Fluid Flow

(

Wiley

,

1953

).

Google Scholar

30.

A.

Sasoh

,

Compressible Fluid Dynamics and Shock Waves

(

Springer Nature

,

Singapore

,

2020

).

Google Scholar

Crossref

31.

G.

Fukushima

,

T.

Tamba

,

A.

Iwakawa

, and

A.

Sasoh

, “

Influence of cellophane diaphragm rupture processes on the shock wave formation in a shock tube

,”

Shock Waves

30

,

545

–

557

(

2020

).

https://doi.org/10.1007/s00193-020-00951-2

Google Scholar

Crossref

32.

G.

Fukushima

et al, “

Losing the shock wave front profile due to interaction with turbulence

,”

Fluid Dyn. Res.

53

,

025504

(

2021

).

https://doi.org/10.1088/1873-7005/abeda4

Google Scholar

Crossref

33.

K.

Kitamura

and

E.

Shima

, “

Numerical experiments on anomalies from stationary, slowly moving, and fast-moving shocks

,”

AIAA J.

57

,

1763

–

1772

(

2019

).

https://doi.org/10.2514/1.J057366

Google Scholar

Crossref

34.

K.

Kitamura

,

Advancement of Shock Capturing Computational Fluid Dynamics Methods Numerical Flux Functions in Finite Volume Method

(

Springer Nature Singapore

,

2020

).

Google Scholar

Crossref

35.

Z. J.

Wang

,

L.

Zhang

, and

Y.

Liu

, “

Spectral (finite) volume method for conservation laws on unstructured grids IV: Extension to two-dimensional systems

,”

J. Comput. Phys.

194

,

716

–

741

(

2004

).

https://doi.org/10.1016/j.jcp.2003.09.012

Google Scholar

Crossref

36.

G. A.

Sod

, “

A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws

,”

J. Comput. Phys.

27

,

1

–

31

(

1978

).

https://doi.org/10.1016/0021-9991(78)90023-2

Google Scholar

Crossref

37.

C.-W.

Shu

and

S.

Osher

, “

Efficient implementation of essentially non-oscillatory shock-capturing schemes, II

,”

J. Comput. Phys.

83

,

32

–

78

(

1989

).

https://doi.org/10.1016/0021-9991(89)90222-2

Google Scholar

Crossref

38.

A.

Kurganov

and

E.

Tadmor

, “

Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers

,”

Numer. Methods Partial

18

,

584

–

608

(

2002

).

https://doi.org/10.1002/num.10025

Google Scholar

Crossref

39.

P.

Woodward

and

P.

Colella

, “

The numerical simulation of two-dimensional fluid flow with strong shocks

,”

J. Comput. Phys.

54

,

115

–

173

(

1984

).

https://doi.org/10.1016/0021-9991(84)90142-6

Google Scholar

Crossref

40.

V.

Daru

and

C.

Tenaud

, “

Evaluation of TVD high resolution schemes for unsteady viscous shocked flows

,”

Comput. Fluids

30

,

89

–

113

(

2000

).

https://doi.org/10.1016/S0045-7930(00)00006-2

Google Scholar

Crossref

Improved hybrid approach of monotonic upstream-centered scheme for conservation laws and discontinuity sharpening technique for steady and unsteady flows

I. INTRODUCTION

II. NUMERICAL SETUP

A. Governing equations

B. Reconstruction methods

1. MUSCL

2. THINC

3. Hybrid MUSCL-THINC scheme (for single-phase flows)

4. Novel hybridization approach of MUSCL and THINC

III. NUMERICAL TEST

A. Weak shock formation problem (unsteady, accuracy study)

B. Isotropic vortex advection (unsteady, accuracy study)

C. Sod shock tube (unsteady, accuracy study)

D. Moving shock (unsteady)

E. Stationary shock (steady)

F. Blunt body problem (steady)

G. Shu–Osher's shock-entropy wave interaction (unsteady)

H. 2D Riemann problem (unsteady)

I. Double-Mach reflection (unsteady)

J. Viscous shock-tube (unsteady, viscous)

K. Computational costs

IV. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: EFFECT OF β IN VISCOUS SHOCK-TUBE PROBLEM

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Improved hybrid approach of monotonic upstream-centered scheme for conservation laws and discontinuity sharpening technique for steady and unsteady flows

I. INTRODUCTION

II. NUMERICAL SETUP

A. Governing equations

B. Reconstruction methods

1. MUSCL

2. THINC

3. Hybrid MUSCL-THINC scheme (for single-phase flows)

4. Novel hybridization approach of MUSCL and THINC

III. NUMERICAL TEST

A. Weak shock formation problem (unsteady, accuracy study)

B. Isotropic vortex advection (unsteady, accuracy study)

C. Sod shock tube (unsteady, accuracy study)

D. Moving shock (unsteady)

E. Stationary shock (steady)

F. Blunt body problem (steady)

G. Shu–Osher's shock-entropy wave interaction (unsteady)

H. 2D Riemann problem (unsteady)

I. Double-Mach reflection (unsteady)

J. Viscous shock-tube (unsteady, viscous)

K. Computational costs

IV. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: EFFECT OF β IN VISCOUS SHOCK-TUBE PROBLEM

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only