In finitevolume methods, monotonic upstreamcentered schemes for conservation laws (MUSCL) offer secondorder spatial accuracy but tend to produce highly dissipative solutions for density discontinuity and weak shock waves. To address this limitation within a secondorder 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 highresolution surpassing fifthorder spatial accuracy schemes within this secondorder spatial accuracy framework with less computational cost. Moreover, the scheme exhibits commendable convergence and robustness when applied to steadystate problems featuring strong shock waves. This scheme offers a more precise and highresolution alternative to conventional MUSCL for compressible flow computations, as it requires no additional stencil for reconstruction, unlike conventional fifthorder schemes.
I. INTRODUCTION
Compressible fluid computations employing shockcapturing schemes within finitevolume methods (FVM) have found widespread use across various applications. The most commonly employed shockcapturing method is the monotonic upstreamcentered scheme for conservation laws (MUSCL) with slope limiters,^{1} typically providing secondorder 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 wavelike shape, losing its discontinuous profile.^{5} In multiphase flow simulations, gas–liquid interfaces suffer from excessive dissipation.^{6} There is a lack of widely accepted secondorder schemes capable of effectively and accurately handling density discontinuities and weak shock waves.
One approach to achieve accurate solutions is through the utilization of higherorder methods that incorporate multiple surrounding cells for reconstructing target cells. For instance, the weighted essentially nonoscillatory (WENO)^{7} scheme provides less dissipative numerical solutions, achieving a maximum fifthorder 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 fifthorder spatial accuracy shockcapturing scheme based on the boundary variable diminishing (BVD) principle^{12,13} has been developed.^{14} This scheme selectively incorporates cell boundary values through a fourthdegree polynomial and the tangential hyperbolic interface capturing (THINC) function.^{15–17} Numerical tests using these higherorder schemes have yielded promising results; however, they entail greater complexity and expense compared to the widely used secondorder spatial accuracy schemes. We maintain that secondorder 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 secondorder framework to discretization using unstructured grids requires less effort compared to higherorder schemes. From this perspective, our study introduces a novel secondorder 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 volumeoffluid method.^{15} 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.^{6} introduced THINC for multiphase computations using a twofluid model^{18} within a finitevolume 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 MUSCLTHINCBVD scheme,^{19} determining cell boundary values based on the BVD principle to minimize numerical dissipation. More recently, Chiu et al.^{20} introduced a hybrid MUSCLTHINC 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 lowdissipation results for singlephase unsteady flows and twophase problems.
Given the success of THINC in multiphase flow computations, its application appears promising for singlephase compressible flow simulations aiming to accurately capture density discontinuities and weak shock waves. However, prior research on schemes combining MUSCL and THINC for solving singlephase 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 MUSCLTHINC scheme by Chiu et al.^{20} While it has been reported that this scheme can provide lowdissipative solutions in singlephase compressible flows, the numerical examples examined in Ref. 20 were limited to a few unsteady flow cases (shocktube and doubleMach reflection problems). Furthermore, its performance in terms of convergence for steadystate problems, crucial for practical applications, was not discussed. In fact, as we will later report in this paper, the original hybrid MUSCLTHINC 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 singlephase steady and unsteady compressible flow computations. As such, we believe that the hybridization strategy of MUSCL and THINC for singlephase 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,^{20} where the hybrid MUSCLTHINC scheme was found to be more efficient than the MUSCLTHINCBVD 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 highresolution and robust secondorder scheme for singlephase 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 “nonlinearityweighted parameter” and the “slope ratioweighted parameter.” The former is sensitive to the magnitude of nonlinearity in the phenomenon, reflecting its selfsharpening 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 MUSCLTHINC scheme.^{20} While the slope ratioweighted parameter is borrowed from the original hybrid MUSCLTHINC scheme,^{20} the proposed method primarily relies on the idea that the scheme automatically attains a physically consistent balance between the phenomenon's selfsharpening nature and the reconstruction process. This feature is crucial for achieving accurate computations across a broad spectrum of singlephase steady and unsteady compressible flow conditions. Thus, the proposed scheme can be viewed as distinct from the original hybrid MUSCLTHINC 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 equations were discretized using cellcentered FVM. For time integration, the fourstage Runge–Kutta^{25} method was employed in this study. SLAU2^{26} was utilized to compute the inviscid numerical flux at the cell boundary. Viscous flux computations were performed using a secondorder central difference.
B. Reconstruction methods
This section provides an overview of the reconstruction methods employed in this study. For simplicity, a onedimensional 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.
1. MUSCL
2. THINC
3. Hybrid MUSCLTHINC scheme (for singlephase flows)
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.^{5} 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,^{26} 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 secondorder 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 secondorder 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 sloperatioweighted parameter $ \zeta i$, a feature from the original hybrid MUSCLTHINC scheme described in Sec. II B 3, which is effective in detecting discontinuities.^{20} 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 “nonlinearityweighted 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.
According to Fig. 3(b), we can observe the behavior of the nonlinearityweighted 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 nonlinearityweighted 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 $ \varphi i + 1 / 2 p \rho $ 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.
MUSCL alone was employed when the variable profile within the stencil exhibited nonmonotonic behavior. The proposed scheme is designated as the “hybrid MUSCLTHINCpρ” (abbreviated as “TMUSCL”) scheme, where “pρ” denotes the usage of pressure (p) and density (ρ) for computing the nonlinearityweighted parameter $ \xi i$. The scheme name is denoted by the superscript “MTpρ” 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 $ \xi 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 MUSCLTHINC scheme by Chiu et al.^{20} and the newly introduced hybrid MUSCLTHINCpρ scheme, we will refer to these schemes as “MUSCLTHINC” and “MUSCLTHINCpρ,” 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 innercell 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 twofluid model) at an interface, a large value of β is preferred.^{6,20} We believe that this is because the interface did not exhibit selfsharpening. 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 MUSCLTHNC and MUSCLTHINCpρ. To achieve this, the weak shock wave formation problem was solved. Figure 5 illustrates the shock wave formation problem described in.^{30} 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.
The spacetime (x–t) diagram depicting density contours for MUSCL and MUSCLTHINCpρ (with β = 2.4) is presented in Fig. 6. In this case, the result for MUSCLTHINC is not shown because, under these weak nonlinear conditions, the nonlinearityweight parameter approaches 1, making no significant difference between MUSCLTHINCpρ and MUSCLTHINC 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 MUSCLTHINCpρ: 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.
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 MUSCLTHINCpρ and MUSCL. At t = 500 [Fig. 7(a)], the compression wave reaches approximately x = 645–665, where a shock wave has not yet formed. MUSCLTHINCpρ with β ≥ 2.0 maintains the initial linear waveform and reduces the thickness of the wave from its initial value. In contrast, MUSCLTHINCpρ 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.
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.
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 MUSCLTHINCpρ (the proposed scheme in this paper) and MUSCLTHINC (the original scheme^{20}). However, as will be demonstrated later in the Appendix, MUSCLTHINC with β = 2.8 exhibits an anomalous shape of vortex in the twodimensional viscous shocktube problem (which is not observed in MUSCLTHINCpρ). To avoid any confusion, we adopt β = 2.4 as the recommended value hereafter.
B. Isotropic vortex advection (unsteady, accuracy study)
Graphs illustrating the $ L 1\u2212$ and $ L \u221e\u2212norms$ plotted against the mesh size are presented in Fig. 9. The results for MUSCL, MUSCLTHINC, and MUSCLTHINCpρ are displayed. All methods demonstrated secondorder accuracy in this problem. This is because all the methods employ minmodlimited MUSCL for continuous flow, which provides secondorder spatial accuracy with linear reconstruction within the cells. MUSCLTHINC and MUSCLTHINCpρ exhibited $ L \u221e\u2212norm$ approximately 0.5 orders of magnitude smaller than MUSCL. Consequently, it is found that using THINC reconstruction yielded higher resolutions for continuous flow. MUSCLTHINC showed smaller errors than MUSCLTHINCpρ 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 lowresolution grids. In finer grids, there was not a significant difference in accuracy between MUSCLTHINC and MUSCLTHINCpρ.
C. Sod shock tube (unsteady, accuracy study)
The density profile at t = 0.2 (after 200 steps) is presented in Fig. 10. The topright panel provides an enlarged view of the moving density discontinuity (contact surface), while the bottomright panel focuses on the region around the shock wave. Additionally, the weight value of THINC ( $ \zeta i$ for MUSCLTHINC, and $ \zeta i \xi i$ in the case of MUSCLTHINCpρ) is displayed in the same figure. The density profiles generated by MUSCLTHINC and MUSCLTHINCpρ were nearly identical, with a slight but noteworthy difference. The shock wave (0.84 < x < 0.86) thickness for MUSCLTHINCpρ was slightly broader than that of MUSCLTHINC, 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 MUSCLTHINCpρ, as evidenced by the smaller value of $ \zeta i \xi i$ compared to $ \zeta 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), MUSCLTHINCpρ accurately captured it, similar to the performance of MUSCLTHINC. 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 steadystate flow simulations.
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)
Figure 11 displays the density profile, computed using MUSCL, MUSCLTHINC, and MUSCLTHINCpρ. 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 MUSCLTHINC and MUSCLTHINCpρ 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.
For M_{s} = 1.5 and 3.0, MUSCL showed a discontinuous profile. MUSCLTHINC and MUSCLTHINCpρ also showed a discontinuous profile. The only difference was the number of cells expressing the shock wave. MUSCLTHINC solved the shock wave with two cells for M_{s} = 1.5 and 3.0, whereas MUSCLTHINCpρ solved it with four cells (M_{s} = 1.5) or three cells (M_{s} = 3.0). As discussed in Sod's shocktube problem (Sec. III C), suppressing THINC at the shock wave using the nonlinearityweighted 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 onedimensional space of [0, 100], where 100 cells were arrayed in the xdirection (Δx = 1). A normal shock wave was installed at the interface between the 50th and 51st cells in the xdirection. 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 walllike boundary used in Ref. 24 was applied to the rightside 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, MUSCLTHINC, and MUSCLTHINCpρ 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 MUSCLTHINC and MUSCLTHINCpρ captured the shock wave sharply. Although a slight postshock oscillation was observed, the implementation of THINC significantly improved the resolution for a very weak shock wave.
MUSCLTHINC exhibited periodic postshock oscillations at M = 1.5 and 3.0. In contrast, MUSCLTHINCpρ yielded less oscillatory results due to the improved hybrid balance between MUSCL and THINC achieved through the nonlinearityweighted 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 highresolution and stable solution. The qualitative convergence performance is evaluated using the 2D bluntbody problem described in Sec. III F.
F. Blunt body problem (steady)
Figure 13 displays the computational grid [Fig. 13(a)] and density contours at 100 000 time steps for MUSCL [Fig. 13(b)], MUSCLTHINC [Fig. 13(c)], and MUSCLTHINCpρ [Fig. 13(d)]. MUSCL [Fig. 13(b)] and MUSCLTHINCpρ [Fig. 13(d)] successfully resolved the bow shock and the downstream flow field without apparent anomalies. However, MUSCLTHINC [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 ( $ \zeta i$ for MUSCLTHINC and $ \zeta i \xi i$ for MUSCLTHINCpρ) used in the density reconstruction in the idirection, along with density contour lines (in black), are illustrated in Fig. 14. The key difference between MUSCLTHINC and MUSCLTHINCpρ in these figures is the weight value of THINC around the shock wave. In MUSCLTHINC [Fig. 14(a)], the weight value of THINC $ \zeta i$ is close to 1 within and behind the shock wave. In contrast, the weight value of THINC in MUSCLTHINCpρ $ \zeta i \xi 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 MUSCLTHINC and MUSCLTHINCpρ.) As observed in previous numerical tests of onedimensional stationary shock waves (Sec. III E), using THINC leads to numerical oscillations behind strong shocks. Therefore, this test reveals that the weight value of MUSCLTHINC for strong shock waves was inappropriate for this steady, twodimensional problem.
Figure 15 shows the density residuals, which represent the maximum value of the density changes between every time step in all computational domains. In MUSCLTHINC, the residual reached 1 × 10^{−2} 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 MUSCLTHINCpρ, the density residual converged to 1 × 10^{−13} at a time step of approximately 15 000. Compared to MUSCLTHINC, 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 MUSCLTHINCpρ was suppressed near the strong shock wave compared to MUSCLTHINC. 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^{−15}. Although the residual of MUSCLTHINCpρ was approximately two orders larger than that of MUSCL, MUSCLTHINCpρ exhibited considerably improved convergence over MUSCLTHINC.
In summary, this test revealed that the proposed hybridization strategy successfully addressed steadystate problems involving strong shock waves. This was achieved owing to successful tuning of the THINC utilization at strong shock waves. The proposed approach, MUSCLTHINCpρ, 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 highresolution feature for weak shock waves.
G. Shu–Osher's shockentropy wave interaction (unsteady)
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^{−4}. MUSCL exhibited significant numerical dissipation in the entropy waves after their interaction with the shock. In contrast, both MUSCLTHINC and MUSCLTHINCpρ resolved the entropy waves after the interaction, comparable to the reference result with four times the number of cells. The difference between MUSCLTHINC and MUSCLTHINCpρ was minimal in this problem, confirming that the nonlinearityweighted parameter does not obviously compromise the resolution of entropy waves, as intended in Sec. II B 4.
H. 2D Riemann problem (unsteady)
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 lowdensity region where x < 0.5 and y < 0.5. MUSCL struggles to resolve the KH instability due to significant numerical dissipation. In contrast, both MUSCLTHINC [Fig. 17(b)] and MUSCLTHINCpρ [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 highdensity region with x < 0.5 and y < 0.5. Figure 17(d) presents the results obtained using WENO^{8} reconstruction for comparison with higherorder schemes. Reconstruction was carried out in a fifthorder componentwise manner. Notably, MUSCLTHINCpρ exhibits superior resolution for density discontinuity compared to the WENO reconstruction. This is because, for the discontinuity, the higherorder scheme is not necessarily optimized: the higherorder 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 secondorder framework, MUSCLTHINCpρ achieves resolution surpassing that of the approximately 26% more expensive higherorder 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.
I. DoubleMach 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.^{39} This problem, known as the doubleMach 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.^{20} previously tested the performance of the original hybrid MUSCLTHINC 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 MUSCLTHINC 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 MUSCLTHINC.
The density contours obtained from the simulations using MUSCL, MUSCLTHINC, and MUSCLTHINCpρ 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, MUSCLTHINC [Fig. 18(b)] exhibited the generation of vortices in the recirculation region, consistent with previous findings.^{20} 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 twodimensional 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 MUSCLTHINCpρ [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 realworld scenarios.^{21–24}
J. Viscous shocktube (unsteady, viscous)
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. Highresolution schemes are known to exhibit longitudinal elongation of the primary vortex at approximately (x, y)= (0.8, 0.05).^{40} The dissipative nature of the MUSCL scheme is evident from the horizontally lengthened vortex shape in Fig. 19(a). In contrast, both MUSCLTHINC [Fig. 19(b)] and MUSCLTHINCpρ [Fig. 19(c)], as well as the reference solution [Fig. 19(d)], demonstrated characteristics of a highorder scheme, particularly the longitudinal elongation of the primary vortex. Remarkably, the results obtained with MUSCLTHINCpρ closely resembled those of the reference solution. Furthermore, we observed that the MUSCLTHINCpρ solution exhibited fewer numerical oscillations compared to MUSCLTHINC. The density profiles along the bottom wall are displayed in Fig. 20. We confirmed that the MUSCLTHINCpρ 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.
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 waveentropy wave interactions (Sec. III G) and the twodimensional 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.
.  MUSCL .  MUSCLTHINC (original) .  MUSCLTHINCpρ (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 MUSCLTHINC and a 29.1% increase for MUSCLTHINCpρ 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 highresolution results with MUSCL, a significantly finer grid and a fourfold reduction in the time increment were required, resulting in an approximately 16fold increase in computational time.
In the context of the twodimensional Riemann problem, MUSCLTHINC exhibited a 20.7% increase in computational time, while MUSCLTHINCpρ 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 MUSCLTHINC and MUSCLTHINCpρ at the density discontinuity.
In summary, MUSCLTHINCpρ resulted in a computational time increase of 29.1%–31.5% compared to MUSCL and 8.87%–13.7% compared to the original MUSCLTHINC. 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 “MUSCLTHINCpρ” scheme, or briefly, TMUSCL, for both steady and unsteady singlephase 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 nonlinearityweighted parameter and a sloperatioweighted parameter. The proposed scheme offers the following distinctive features.

Enhanced accuracy for continuous flow simulations, providing secondorder 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 twodimensional blunt body problems. The density residual decreased significantly, surpassing ten orders of magnitude reduction when compared to the original MUSCLTHINC hybrid scheme.

Mitigation of shock anomalies in the presence of strong shock waves while maintaining highresolution computations for density discontinuities and entropy waves. Remarkably, the proposed scheme achieved results with resolution exceeding that of a more computationally intensive fifthorder reconstruction at slip lines in the twodimensional Riemann problem, despite being designed within a secondorder 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 nonlinearityweighted parameter based on pressure and density, it is multiplied by the slope ratioweighted parameter used in the original hybrid MUSCLTHINC scheme. Additionally, the stencil size for reconstruction aligns with that of the secondorder spatial accuracy MUSCL approach. The proposed scheme has some drawbacks: higher computational costs than that of MUSCL and the original hybrid approach, slightly insufficient smallscale resolution, and somewhat compromised convergence properties than MUSCL. Nevertheless, the proposed scheme consistently delivers robust, highresolution 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 realworld 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 weakshock 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 SHOCKTUBE PROBLEM
As discussed in the onedimensional weak shock formation problem (Sec. III A), β values of 2.4 and 2.8 were found to yield physically consistent compression trends for MUSCLTHINC and MUSCLTHINCpρ. In this Appendix, we further assess the impact of β on the viscous shocktube problem (Sec. III J). Figure 21 displays the density contours of the viscous shocktube problem, solved using MUSCLTHINC with β = 2.4 [Fig. 21(a)] and 2.8 [Fig. 21(b)], MUSCLTHINCpρ with β = 2.4 [Fig. 21(c)], and β = 2.8 [Fig. 21(d)]. Notably, when β was set to 2.8 for MUSCLTHINC, 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 MUSCLTHINCpρ 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 MUSCLTHINC computations.