Research Article  Open Access
Hybrid Flux Method in MonotonicityPreserving Scheme for Accurate and Robust Simulation in Supersonic Flow
Abstract
The fifthorder monotonicitypreserving (MP5) scheme is an accurate and low dissipative numerical method. As a finitevolume method, MP5 adopts the Roeflux scheme for solving the numerical flux in the compressible Euler equation. However, due to the deficiency of the MP limiter and Roeflux in maintaining positive density and pressure, the calculation could fail in cases of extreme flow involving small values of density and pressure. In this study, to overcome such a limitation but still to achieve a highaccuracy of MP5, we propose a hybrid flux method: the Roeflux is used in the global computational domain, but the firstorder LaxFriedrich (LF)flux is adopted only for trouble grids. The numerical results of shocktube and complicated interaction problems indicate that the present scheme is more accurate at discontinuities and local extrema compared to the previous scheme, maintaining positive density and pressure values. For twodimensional applications, a supersonic jet is explored with different Mach numbers and temperature conditions. As a result, small vortices induced by the shear layer can be clearly captured by the proposed scheme. Furthermore, a simulation was successfully conducted without blowup of calculation even in the extreme jet flow condition.
1. Introduction
The usage of an efficient highorder shockcapturing scheme is important in the large eddy simulation (LES)/direct numerical simulation (DNS) particularly in supersonic flow conditions because of the requirements to capture a sharp discontinuity and simultaneously to resolve small structures of turbulence. In addition, supersonic flow sometimes causes very low density and pressure values that have the potential to induce failure of calculation. Therefore, the numerical robustness is an additional requirement for a stable simulation.
The wellknown essentially nonoscillatory (ENO) [1] and weighted ENO (WENO) schemes [2] are very robust highorder method to capture the shockwave without introducing any nonphysical oscillation; however, they have been proven as a large numericaldissipative scheme. The weighting terms of WENO scheme degenerate the accuracy in nondiscontinuous solution fields by allowing evident dissipation. The large dissipation yields the lack of performance for resolving small flow structures.
On the other hand, following several classical limiting strategies [3, 4] to prevent nonphysical oscillation, fifthorder monotonicitypreserving (MP5) scheme [5] was proposed by Suresh and Huynh. In this method, the Roeflux scheme [6] was adopted as the basis method for solving numerical flux of compressible Euler equation. The welldesigned limiting part of MP5 scheme enables preserving the monotonicity at discontinuities and the accuracy at local extrema. Consequently, accurate results can be obtained in complicated flow physics such as shockwave and smallscale turbulent flows. Hence, the MP5 scheme has been applied for the DNS and LES of compressible turbulence interacting with and without shock waves [7, 8]. In addition, the MP5 scheme has even also been used for nonlinear acoustic wave simulation [9]. Several attempts have been made toward improving the accuracy of the MP5 scheme. Li et al. [10] increased the stencil size of interpolation in order to enhance the order of accuracy and found the proper coefficients of interpolation to reduce numerical dissipation. Fang et al. [11, 12] designed the interpolation in a similar way, but they considered a shocksensor for an efficient adaptation of an MP limiter. In addition, although not related to accuracy, Daru and Tenaud [13] proposed a onestep scheme by reinterpreting the MP limiter as a TVDlike one.
However, it was found that the calculation using MP5 scheme could blow up in specific situations when the flow involves very low density or pressure. On the region of local minimum value, nonphysical negative density or pressure can be generated during simulation. He et al. [14] firstly dealt with this problem in MP schemes and discussed the cause, which is a small overshoot/undershoot oscillating solution behavior due to the weakness of the MP limiter. Therefore, to overcome such limitation, they considered the following modifications for MP5 scheme: (i) total variation diminishing (TVD) condition [15] was implemented to the limiter part of MP5 scheme to prevent overshoot/undershoot oscillating solution behavior, and (ii) LaxFriedrichs (LF)flux [16] was used for the basis flux method because it has a positivitypreserving capability (Roeflux used in the MP5 scheme has no positivitypreserving capability [17]). They successfully showed the robustness of their method (MP5R) in several benchmark problems which make the calculations using the original MP5 scheme fail. However, their method has disadvantage related to the overall accuracy. Since the LFflux method is a highly dissipative [18, 19], the overall accuracy is degenerated from that of the original MP5. For reference, a few limiters have been proposed to solve the negative problem in the improved versions of WENO and discontinuous Galerkin (DG) schemes [20–22]. However, those methods also rely on the dissipative LFflux method.
In this study, we try to develop a new efficient MP scheme to solve the negative density and pressure problem and simultaneously to achieve the highaccuracy. The key strategy of present scheme is the hybridization of Roe and the firstorder LFfluxes depending on local density and pressure state. Initially, the solution in the whole flowfield is obtained using the original MP5 scheme; however, only for the local trouble region detected by criteria, the firstorder LFflux is applied to modify the numerical flux. The firstorder never introduces the nonphysical oscillation; therefore, the additional TVD condition used in MP5R is not necessary.
The order of accuracy test is conducted to assess a global level of error and convergence in a different number of grids. After that, the shocktube and complicated interaction problems with normal/extreme conditions are carried out to evaluate the performances of schemes in various situations. In normal conditions, we compare the accuracy between MP5R and MP5Hyb scheme. However, in the extreme conditions involving low density and pressure, the robustness of scheme can be investigated as well. For the application, supersonic jet flow simulations are considered, which require the same numerical properties of highaccuracy and robustness. Initially, the position of initial vortex induced by a cold jet with Mach 1.56 was investigated to verify numerical solutions and to compare the resolving performance for small flow structures. In addition, a Mach 3.0 jet with an extremely high temperature of 2,100K was considered to ensure the robustness of scheme in such tough condition. Note that the solutions of original MP5 just before blowup of calculation are presented as well in tests involving low density and pressure.
2. Background Methodology of MP Scheme
2.1. Discretization
Consider the methods for onedimensional linear advection equation with speed of wave
For the moment, assume that positive wave speed . Let be grid size and be the index of cellcentered grid. After integrating (1) over the grid , it yieldsor
For accurate calculations, the highorder interpolation for the numerical flux is performed. In this study, we consider the monotonicitypreserving (MP) scheme framework as introduced in the next section.
2.2. MonotonicityPreserving (MP) Scheme Framework
The procedure of MP scheme is separated in two parts: (i) highorder interpolation and (ii) application of MP limiter (or constraints) for shockcapturing. The former procedure enhances the order of accuracy by obtaining the cellinterface flux using neighboring grids. Let us consider the cellaveraged variable on an arbitrary grid as , then the fifthorder upwind interpolation adopted in the MP scheme is given by [5],where , and . Equation (4) is expressed for the rightgoing wave. For the opposite direction (), the formulation has the symmetrical form. However, with only (4), nonphysical oscillation appears near a discontinuous solution. Therefore, the limiting procedure is performed for a clear shockcapturing capability and stable solution. Some constraints generate the local specific interval for every cell interfaces. Near the discontinuous region, is changed to one of the bound values or to satisfy the monotonicitypreserving condition. On the other hand, the constraints do not alter the original one so that the accuracy of (4) is preserved. The local interval is constructed by the following terms:whereand is defined as Here, indicates the local measure of curvature as follows:
The minmod function is nothing but selecting the minimum if signs of terms are same; otherwise, it gives zero. The definition of others and details can be found in [5].
The error characteristics of the scheme can be investigated by Fourier analysis (or modified wavenumber analysis). The relation determining the spatialaccuracy can be obtained by combining (4) with the spatialderivative term (without constant c) of governing equation as
Then Fourier transform is adapted to estimate spectral properties for wide range of wavenumbers. Since we can simply obtain from analytical differential aswe can deriveFrom this equation, the (scaled) modified wavenumber can be derived as
The wavenumber (15) can be decomposed into real and imaginary parts which provide dispersive and dissipative behaviors of solution, respectively. Figure 1 summarizes dispersion and dissipation characteristics of MP scheme, and those of and order schemes are also presented for references. In particular, the values in the dissipation part are negative for each corresponding wavenumber , which means that the amplitude of solution always decreases; thus, the MP scheme is designed to be spatially stable.
(a)
(b)
2.3. Extension to Euler Equation
Consider the onedimensional compressible Euler equationwhere are density, velocity, and pressure, respectively. The total energy can be expressed with the internal energy as E =. In the same manner with (1)(3), it yieldswhere is the cell average. Note that upper hat is missing on the numerical flux term for simplicity. Since the system of numerical flux is nonlinear, we need a specific method such as wellknown Roe [6] or LaxFriedrich (LF)flux [16] to obtain a solution. In general, the procedure of highorder method is usually carried out in the characteristic field. Let us take as the Jacobian matrix. Then, the system of equation can be decoupled by multiplying the left eigenvector () of Jacobian matrix as follows:where is the righteigenvector of Jacobian matrix. The details of eigenvectors can be found in the Appendix. Then, (18) can be simplified with the characteristic variable aswhere and . is the diagonal matrix composed of eigenvalues of Jacobian matrix as . Meanwhile, the flux of (16) can be split according to the sign of flux aswhere and . The matrices and are nothing but the positive and negative components of the Jacobian matrix, respectively, i.e., and . Here, and . Therefore, (20) also can be decoupled by multiplying the left eigenvector aswhere and .
3. Different Flux Methods in Previous MP Schemes
3.1. RoeFlux in Original MP5 Scheme for HighAccuracy
The original fifthorder monotonicitypreserving (MP5) scheme proposed by Suresh and Huynh [5] was developed based on the Roeflux [6] which is the approximated version of Godunov flux scheme. The merit of Roeflux is the highaccuracy and good resolution of shock waves because it considers not only the direction of characteristic wave propagation, but also the waves themselves at the cellinterface [18, 19].
The procedure of MP5 scheme can be summarized as follows: (i) calculate the characteristic variable shown in (19) on stencils for j+1/2 flux as (ii) Perform the procedure of (4)~(9) on each vector of in left and right direction at the cellinterface to obtain and . (iii) Return to the physical space by multiplying the righteigenvector: and . (iv) Left and right directional conservative variables, and , are used in the Roeflux asThe differential of the characteristic vector is denoted by as follows:
The notation stands for Roeaveraged values between left and right values aswhere is the total enthalpy as . However, the calculation using the MP5 scheme could fail on the region involving very low density and pressure owing to the following reasons:
Problem I. Density and pressure values at the cellinterface ( and ) obtained by MP scheme could be negative due to undershoot oscillations [14]. This could occur when the averaged values on stencils ( and , ) are slightly positive.
Problem II. Roeflux scheme has no positivitypreserving capability [17]. It has a potential to give negative cellaveraged density and pressure ( and ) even though slightly positive cellinterface values ( and ) are used in the Roeflux.
3.2. LaxFriedrich Flux in MP5R Scheme for Robustness
To overcome the aforementioned limitation of original scheme, MP5R [14] was suggested and successfully verified its robustness in severe numerical test problems. In this method, to solve Problem I, the total variation diminishing (TVD) [15] condition is imposed in the limiter part of MP5 scheme to strictly prevent nonphysical oscillations. In addition, to solve Problem II, the Roeflux was substituted with the LFflux scheme which has positivitypreserving capability [19]. As a flux splitting method, the LFflux is monotone itself, which means never introducing nonphysical oscillation so that density and pressure values are always positive. However, the LFflux is quite dissipative than Roeflux [18, 19]. Therefore, the overall accuracy of MP5R scheme is degenerated from that of the original MP5 scheme.
The MP5R scheme based on the LFflux can be summarized as follows: (i) Calculate the characteristic variable shown in (21) on stencils for j+1/2 flux asHere, is obtained using LFflux as (ii) Perform the procedure of MP scheme shown in (4)~(9) on each vector of in left and right direction at the cellinterface to obtain and . (iii) Return to the physical space by multiplying the righteigenvector: and . (iv) Then, the final numerical flux can be obtained as
4. Hybrid Roe and FirstOrder LFFlux for HighAccuracy and Robustness
As introduced in the previous section, MP5 and MP5R schemes have different numerical characteristics by adopting a different flux scheme. The former is quite accurate due to using Roeflux; however, it is not robust in low density or pressure conditions. In contrast, the latter is robust; however, the accuracy is degenerated from the original MP5 scheme due to dissipative LFflux.
However, in practice, the numerical problem owing to low density or pressure does not appear in the entire computational domain. Therefore, on the same MP scheme framework, we suggest the efficient method as follows: use the LFflux only for trouble grids to modify the negative density and pressure; otherwise, apply the Roeflux to obtain accurate solution. By implementing such a hybridconcept, we try to achieve both of the respective advantages of the MP5 and MP5R schemes.
Let us explain the detailed methodology of the hybrid flux method. If we know the solution at th time step, then the residue shown in the second term of (17) can be expressed as . By giving information of the fluxin/out on the jth cell, we can move to the next time step. At first, to obtain the accurate solution, residues on the entire grids are calculated by using MP5 scheme with Roeflux.where is the total number of grids. After that, regardless of local flow state, just perform a timemarching to obtain numerical solutions in the next time step as . In this step, if the entire flowfield has no problem, then the calculation is straightforward. However, the calculation is failed when negative density or pressure is observed in at least one local grid. Therefore, the modification of flux for the trouble grid is necessary. We use the following condition to find the trouble grids at the n+1th time step:Note that we used the asterisk mark () as superscript to specify the index of the trouble grid. The density and pressure are directly related to several terms such as speed of sound , three eigenvalues , eigenvectors of Jacobian matrix, and so on. Therefore, negative values of density and pressure make the system equation of (16) to be nonhyperbolic and illposed. Furthermore, in practice, negative density or pressure at the cellinterface ( and ) can give a NaN (Not a Number) during the Roeflux calculation (see, Roeaveraged values in (25)); therefore, the following condition is additionally considered to find trouble grids.
Indeed, we do not need to recalculate numerical fluxes for all grids because solutions are normal except the trouble grids. Therefore, we call the solutions at the nth time step and calculate residue only for the trouble grids. As mentioned before, we try to use the LFflux to obtain new flux values in this step. However, unlike the MP5R scheme, we adopt the firstorder of LFflux which does not require the highorder procedures as described in Section 2.2. The cellaveraged conservative and flux values at j and j+1 are taken for the left and right directional values of (28), respectively, as
Here, the candidates for the conservative and flux terms are changed so that can be obtained as
Let us explain the positivity of density and pressure [20]. As mentioned before, the solution in the next time step (n+1) can be obtained through the residue.Therefore, using the firstorder LFflux, the solution can be expressed as
We can see that is determined by three set of vectors and . Under the condition of CourantFriedrichsLewy (CFL) number <1.0, the first term on the righthandside of (35) is positive. In addition, the density and pressure for the rest of the terms are unconditionally positive owing to the following reasons (for simplicity, we drop the notations): for the density in the first component of vectorand for the pressureThus, it should satisfy
Equation (38) is always valid. Therefore, has positive density and pressure values.
For all these reasons, on a local trouble grid, Problem I and Problem II which were discussed in Section 3.1 are solved automatically. Through this approach, the additional numerical treatment such as TVD condition adopted in MP5R is not required anymore. While the accuracy for local trouble grids may seem low, the simulations we are targeting are not dominated by the active behaviors of vacuum flows. Rather, we are interested in the flows which are dominated by shockwave/turbulence interaction, sometimes involving low density and pressure on few local regions, e.g., supersonic jet flow (which will be discussed in Section 6). To this end, we pursue the numerical scheme that represents highaccuracy in overall domain while maintaining a stable calculation giving reasonable positive density and pressure values.
Therefore, we can calculate a new local residue asand finally we obtain after the timeintegration. Henceforward, MP5Hyb (Hyb is the abbreviation for Hybrid) refers to the present scheme, and its performance will be proved in numerical validations in Section 5.
Figure 2 shows the schematic for the hybrid flux method. The notation indicates the hybridfactor which allows the modification for neighbors of the trouble grid. It is adjustable by users; however, in this study we have used the simplest case .
In general, the multistage RungeKutta method is used for the timeintegration. In this case, the modification procedures for negative density or pressure should be applied in every internal stages of RungeKutta method. Figure 3 summarizes the procedures of the hybrid flux method with kth stages RungeKutta method as a flow chart.
Implementation in a multidimensional case is quite straightforward. The procedures for flux modification are adopted in direction x, y, and z. One thing to note is that the numerical fluxes at all grid interfaces (four interfaces in 2D and six interfaces in 3D) should be obtained for one trouble grid. Therefore, for example, in a twodimensional problem, the residue is changed as follows:where i and j indexes stand for x and y coordinates, respectively.
5. Numerical Validations
In this section, several problems related to the normal and extreme flow conditions have been investigated for accuracy and robustness of the present method. Classical test problems in onedimensional cases and supersonic jet flow simulations in asymmetrical twodimensional domain are performed, which is the main application of the present scheme. For all the following computations, a thirdorder TVD RungeKutta timeintegration method [2] is adopted. Regardless of the flux scheme, a CFL number limitation of MP scheme for stable calculation is determined by the coefficient of (7). Theoretically, should be satisfied with RungeKutta timeintegration; however, still yields a stable solution [5, 12]. In this study, the CFL number is set to 0.1. If not mentioned otherwise, the specific heat ratio of standard air is used as . Note that the solutions of original MP5 and MP5Hyb are the same unless negative density or pressure appears. Therefore, for a normal condition, the results of original MP5 are excluded for brevity.
5.1. Order of Accuracy and Convergence Test
A typical way to verify the order of accuracy and convergence of numerical scheme is to simulate the advection of wave and compare the standard errors of solutions in different number of grid situations. For the compressible Euler equation, on the domain , a density wave with a constant freestream can be given by [23]A periodic boundary condition was enforced on both boundaries. Therefore, the location of the density wave is the same as that of the initial condition in every time cycle. At the first cycle t=2, we compared the norm value of the error which is given byWe progressively increased the number of grids from N=15 to N=240 and plotted the corresponding errors as can be seen in Figure 4. Note that each axis was set by a logarithm scale. We can observe that MP5R and MP5Hyb schemes have the almost the same slope as five (slope indicates the order of accuracy); however, the absolute error values of MP5Hyb is less than MP5R scheme in all different number of grids. Such results are acceptable because MP5Hyb scheme includes Roeflux whereas MP5R takes only dissipative LFflux scheme. Details of error and its slope are listed in Table 1.

5.2. ShockTube Problems
Shocktube problem demonstrates the dynamics of two different state of gas generating a shockwave, contactdiscontinuity, and expansion fan. Through this test, we can compare the accuracy of discontinuous solution. We consider two different initial conditions of normal and extreme cases, and all numerical solutions were compared with the analytical solution.
5.2.1. Sod ShockTube Problem (Normal Condition)
This is a normal condition case of shocktube problem. On the domain , the initial condition was given as [24]The number of grids was 80 and a zerogradient boundary condition was imposed on both sides to allow the flow naturally leaving out from the domain. The numerical solution was obtained at t=0.4. Numerical solutions of density obtained by MP5R and MP5Hyb are shown in Figure 5. As a result, both schemes gave a good trend of solutions; however, at the contactdiscontinuity (in the enlarged figure) the result of MP5Hyb is slightly close to the analytical solution than the MP5R. Aforementioned, Roeflux part in the hybrid flux method makes such differences. Figure 6 represents the solutions of density with coarse and fine grids. The trend of accuracy according to the number of grids is coherent; thus the gap between MP5R and MP5Hyb is similar in coarse and fine grids. The gradient of solution can be steepened with a fine grid. Since the Sod shocktube problem does not involve low density and pressure values, the local modification by the firstorder LFflux was not activated in the MP5Hyb scheme.
(a)
(b)
5.2.2. Leblanc ShockTube Problem (Extreme Condition)
The LeBlanc problem [20] is the extreme case of the shocktube problem. This test involves an extremely low density and pressure state at the initial stage; therefore, the calculation of a nonrobust scheme will fail. The initial condition was given as follows:On the domain filled with a perfect gas with the specific heat ratio . Total 400 grids were used and a zerogradient boundary condition was imposed on both sides. Numerical solutions were obtained at t=6. Figure 7 depicts the numerical solutions on density and pressure. In particular, the values presented in this problem are very small; thus, yaxis is expressed in the logarithmic scale. The calculations using MP5 scheme do not hold up due to the tough initial condition. However, it is observed that MP5R and MP5Hyb schemes can solve this problem without blowsup of calculation. On the other hand, the accuracy between MP5R and MP5Hyb is quite different. In the density plot, the solution in the flat region between x ~ 6.0 and 6.8 is more accurately predicted with MP5Hyb scheme. In addition, the location of shockwave presented at x ~ 8 in density and pressure plots is more precise in the case of MP5Hyb. Figure 8 summarizes the changes of density solutions with different number of grids. In both cases of coarse and fine grids, MP5Hyb shows the better location of shockwave than MP5R. However, even with the finer grid, the solution of MP5R is not much improved. These results confirm that MP5Hyb scheme can be an efficient numerical method to solve shockwave even in the extreme flow case.
(a)
(b)
(a)
(b)
5.3. Complicated Interaction Problems
In this subsection, we consider two different interaction problems. One is shock and wave interaction, and the other is two shock waves interaction. In these problems, there is no analytical solution because the interaction behaviors are complicated and highly nonlinear. Therefore, we produce a reference solution using the WENO scheme [2].
5.3.1. ShockEntropy Wave Interaction Problem
The simulation for shockentropy wave interaction problem [1] was carried out in this subsection. From this test we can assess the capability of shockcapturing and resolving accuracy on a highfrequency local extrema region. On the domain , the initial condition was given as follows:
A shockwave of Mach=3.0 interacts with the linear entropy waves. A zerogradient boundary condition was imposed on both sides, and total 200 grids were used. The results were compared at the simulation time of t=1.8. For the comparison, the solution obtained by WENO scheme with 2,000 grids was used as a reference solution. As can be seen in Figure 9, both results follow the overall trend of the reference solution. In particular, we present enlarged figures on two parts of the solution: (a) Nshape discontinuities and (b) highfrequency region. On the region (a), the solution of MP5R scheme is overestimated at the peaks of Nshape whereas MP5Hyb scheme presents a clear solution fit on the reference solution. In particular, on the region (b), we can confirm that the accuracy of MP5Hyb is obviously improved compared to MP5R by producing the better smooth peak solutions. Figure 10 demonstrates the solutions of density with coarse and fine grids. Compared to the grid studies of shocktube tests, the accuracy differences are dramatic because the solution in a highfrequency region is smooth. N=150 is too small (only five points given per wave); thus both schemes are hard to resolve the local extrema. The fine grid case, N=250, improves the peak to peak solutions on highfrequency region, and we can still observe the accuracy differences between MP5R and MP5Hyb.
5.3.2. TwoBlastWave Interaction Problem
We considered the twoblastwave interaction problem [25]. In particular, this test includes very low density and pressure state in the shockshock interaction moment, and they could lead to the failure of calculation. The initial condition was given by
On the domain , total 400 grids were used and a reflective boundary condition was imposed on both sides. Therefore, three different state of gas initially generate two strong shock waves, and they move face to face after reflection from the wall boundaries. A simulation time to compare numerical solutions was t=0.038. The solution obtained by WENO scheme with 4000 grids was used as a reference solution.
Figure 11 shows the density and pressure distributions. We can observe that MP5R and MP5Hyb schemes give reasonable results without blowsup of calculation while the MP5 scheme does not maintain its calculation at t=0.027 (the results of blowsup moment are presented in plots). In particular, the solution of MP5Hyb is much accurate than MP5R on the local extrema of density distribution.
(a)
(b)
Figure 12 shows the changes of density solutions with different numbers of grids. With a coarse grid, the MP5Hyb presents slightly oscillating solution near x=0.71, but the peaks of local extrema can be much more accurately predicted than MP5R. With a fine grid, the maximum values of extrema can be predicted by MP5Hyb whereas grids seem still insufficient for MP5R. These facts ensure that the hybrid flux method of present scheme achieve the robustness in low density and pressure state and simultaneously improve the accuracy of solution compared to the MP5R.
(a)
(b)
6. Application to Supersonic Jet Simulation
The main purpose of developing the present scheme is to solve the supersonic jet flow that requires the numerical capabilities we have investigated because (1) the supersonic jet flow contains the shockwave and shear layer to be solved using the method of highaccuracy for a good quality of solution, and (2) low density and pressure appear at the core of vortices induced by highspeed jet flows.
First, we assess that MP schemes successfully simulate the supersonic jet flow by comparing with experiment data, and we compare the detailed flow structures. After that, we evaluated the robustness of present scheme in a highspeed supersonic jet case which generates the vortices with very low density and pressure. The governing equation is the twodimensional axisymmetric compressible Euler equation,where
6.1. Mach 1.56 Supersonic Cold Jet Simulation
We tried to validate and evaluate the performance of MP schemes for the supersonic jet flow. Therefore, the pulse jet problem [26] was chosen as the validation case, and we investigate the flow structures and position of the initial vortex to be compared with the experimental data [22]. For this test, on the twodimensional domain , a grid size was used. Figure 13 summarizes the boundary conditions and domain size for computation. In particular, ideal gas isentropic relations were used to obtain the flow properties at the jet exit boundary [27]. The jet flow condition was cold jet with Mach number = 1.56. Figure 14 shows a density contour at a certain initial stage (t=100ms). This figure introduces the initial structures of supersonic jet such as the shock, reflected shock, and initial vortex. Figure 15 represents the position of the initial vortex along simulation time, and the other simulation results are also referred for comparison [26, 27]. As a result, a good agreement with the experimental data is confirmed by the MP5R and MP5Hyb schemes.
(a)
(b)
Figure 16 shows the comparison of instantaneous density contours obtained by MP5R and MP5Hyb at t=500 ms. Although the location of initial vortex is similar, the whole flow structures are quite different. In particular, clear vortex structures induced by the shear layer can be captured using the MP5Hyb scheme while they are smeared due to large numerical dissipation in the case of MP5R scheme. These facts confirm that the Roeflux scheme in the hybrid flux method increases the resolving performance for smallscale flow structures.
(a)
(b)
6.2. Mach 3.0 Supersonic Hot Jet Simulation
The present scheme was applied for a highMach number and hightemperature supersonic jet simulation to investigate the robustness of scheme. When jet flow speed is high, low density and pressure can be observed at the core of initial vortex and small vortices induced by shear layer. This can lead to negative density or pressure value, and thereby, the calculation may fail unless an adequate numerical treatment is adopted. The computational domain has been changed a little from the previous problem. The domain size in xdirection was doubled owing to higher jet flow speed, and a stretching grid was generated to make the actual computing environment ( and used at the nozzle lip). The total number of grids was about 35,000. The nozzle exit Mach number () and fully expanded Mach number () were set by and , respectively. The total temperature ratio of nozzle and ambient state was 7.
Figure 17 shows the results of MP5 and MP5Hyb schemes at very initial time. At ( is a nondimensional time), the calculation of MP5 scheme is blowup at the core of initial vortex due to negative pressure caused by the fast rotating flow. On the other hand, MP5Hyb can overcome this situation, and maintains the calculation. To investigate another blowup location at a different time, we kept the calculations of both schemes. However, as the calculation of the MP5 scheme was fail at the initial stage, the flowfield of MP5Hyb scheme at was used as the initial condition of MP5 simulation. At , we observed that the MP5 scheme again fails as can be seen in Figure 18; however, the blowup location is not at the core of initial vortex. Hence, not only the strong initial vortex but also the flows induced by shear layer can be the cause of negative density or pressure. On the other hand, the calculation of MPHyb is still stable. Figure 19 presents density and pressure contours of fully developed jet flow obtained by the MP5Hyb scheme (at ). Such flowfield can be obtained by overcoming lots of blowup situations; therefore, we confirm that MP5Hyb scheme can be a useful tool to simulate the supersonic jet flow with the extreme condition.
(a)
(b)
(a)
(b)
7. Conclusions
In this study, MP5Hyb scheme was proposed for robust and accurate simulation in supersonic flow conditions particularly involving low density/pressure values. The original MP5 scheme based on the Roeflux is an accurate shockcapturing scheme; however, the calculation could fail when very low density/pressure appears. MP5R scheme was proposed to overcome this limitation by adopting the LFflux and additional TVD condition. However, the LFflux scheme is more dissipative than the Roeflux scheme; therefore, the overall accuracy is degenerated from the original MP5.
To overcome such limitations but still achieve a highaccuracy of MP5, we have proposed a hybrid flux method. The strategy is that the Roeflux is used in the global computational domain, but the firstorder LFflux is adopted only for trouble grids resulting in negative density or pressure. The firstorder LF never introduces numerical oscillation, and it has a positivitypreserving capability itself; therefore, the additional treatment such as TVD condition is not necessary.
The results of shocktube and complicated interaction problems proved that the MP5Hyb scheme has the robustness without resulting in calculation failure. In addition, it is beneficial to capture the shock waves and local extrema. The MP5Hyb scheme was applied to simulate supersonic jet flows in different conditions. As a result, the MP5Hyb scheme has better accuracy than the MP5R scheme for smallscale structures of vortices. In addition, the MP5Hyb scheme can maintain a stable calculation even with the extreme jet flow condition whereas the MP5 scheme failed. Therefore, we conclude that MP5Hyb scheme can be an efficient tool to simulate the supersonic flow involving low density or pressure.
Appendix
Eigen Vectors of Jacobian Matrix
Left eigenvector of Jacobian matrix can be expressed aswhere and . is a specific heat ratio of gas. The righteigenvector is nothing but the inverse of as
Data Availability
All data can be accessed in the Numerical Validations section of this article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by Space Core Technology Development Program through the National Research Foundation of Korea (NRF), which is funded by the Ministry of Science, ICT & Future Planning NRF2014M1A3A3A020348. This work was also conducted at HighSpeed Vehicle Research Center of KAIST with the support of the Defense Acquisition Program Administration and the Agency for Defense Development under Contract UD170018CD.
References
 C.W. Shu and S. Osher, “Efficient implementation of essentially nonoscillatory schemes,” Journal of Computational Physics, vol. 77, pp. 439–471, 1988. View at: Publisher Site  Google Scholar  MathSciNet
 G. S. Jiang and C. W. Shu, “Efficient implementation of weighted ENO schemes,” Journal of Computational Physics, vol. 126, no. 1, pp. 202–228, 1996. View at: Publisher Site  Google Scholar  MathSciNet
 B. Van Leer, “Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov's method,” Journal of Computational Physics, vol. 32, pp. 101–136, 1979. View at: Google Scholar
 P. Colella and P. R. Woodward, “The piecewise parabolic method (PPM) for gasdynamical simulations,” Journal of Computational Physics, vol. 54, no. 1, pp. 174–201, 1984. View at: Publisher Site  Google Scholar
 A. Suresh and H. T. Huynh, “Accurate monotonicitypreserving schemes with RungeKutta time stepping,” Journal of Computational Physics, vol. 136, no. 1, pp. 83–99, 1997. View at: Publisher Site  Google Scholar  MathSciNet
 P. L. Roe, “Approximate Riemann solvers, parameter vectors and difference schemes,” Journal of Computational Physics, vol. 43, pp. 357–372, 1981. View at: Publisher Site  Google Scholar  MathSciNet
 A. Jammalamadaka, Z. Li, and F. A. Jaberi, “Subgridscale models for largeeddy simulations of shockboundarylayer interactions,” AIAA Journal, vol. 51, no. 5, pp. 1174–1188, 2013. View at: Publisher Site  Google Scholar
 Z. Li and F. A. Jaberi, “A highorder finite difference method for numerical simulations of supersonic turbulent flows,” International Journal for Numerical Methods in Fluids, vol. 68, no. 6, pp. 740–766, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 C. ChiHoon, L. DuckJoo, J. P. Breard, and J. H. Lan, “Timedomain simulation of nonlinear acoustic propagation in a lined duct,” AIAA Paper, pp. 2008–2830, 2008. View at: Google Scholar
 X.l. Li, Y. Leng, and Z. He, “Optimized sixthorder monotonicitypreserving scheme by nonlinear spectral analysis,” International Journal for Numerical Methods in Fluids, vol. 73, no. 6, pp. 560–577, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 J. Fang, Z. Li, and L. Lu, “An optimized lowdissipation monotonicitypreserving scheme for numerical simulations of highspeed turbulent flows,” Journal of Scientific Computing, vol. 56, no. 1, pp. 67–95, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 J. Fang, Y. Yao, Z. Li, and L. Lu, “Investigation of lowdissipation monotonicitypreserving scheme for direct numerical simulation of compressible turbulent flows,” Computers & Fluids, vol. 104, pp. 55–72, 2014. View at: Publisher Site  Google Scholar
 V. Daru and C. Tenaud, “High order onestep monotonicitypreserving schemes for unsteady compressible flow calculations,” Journal of Computational Physics, vol. 193, no. 2, pp. 563–594, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 Z. He, Y. Zhang, F. Gao, X. Li, and B. Tian, “An improved accurate monotonicitypreserving scheme for the Euler equations,” Computers & Fluids, vol. 140, pp. 1–10, 2016. View at: Publisher Site  Google Scholar  MathSciNet
 A. Harten, “High resolution schemes for hyperbolic conservation laws,” Journal of Computational Physics, vol. 49, no. 3, pp. 357–393, 1983. View at: Publisher Site  Google Scholar  MathSciNet
 P. D. Lax, “Weak solutions of nonlinear hyperbolic equations and their numerical computation,” Communications on Pure and Applied Mathematics, vol. 7, no. 1, pp. 159–193, 1954. View at: Publisher Site  Google Scholar  MathSciNet
 B. Einfeldt, C.D. Munz, P. L. Roe, and B. Sjogreen, “On Godunovtype methods near low densities,” Journal of Computational Physics, vol. 92, no. 2, pp. 273–295, 1991. View at: Publisher Site  Google Scholar  MathSciNet
 J. Blazek, Computational Fluid Dynamics: Principles and Applications, Elsevier, 2001.
 E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, 2009. View at: Publisher Site  MathSciNet
 X. Zhang and C.W. Shu, “On positivity preserving high order discontinuous galerkin schemes for compressible euler equation on rectangular meshes,” Journal of Computational Physics, vol. 229, pp. 8918–8934, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 X. Y. Hu, N. A. Adams, and C.W. Shu, “Positivitypreserving method for highorder conservative schemes solving compressible Euler equations,” Journal of Computational Physics, vol. 242, pp. 169–180, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 Y. Guo, T. Xiong, and Y. Shi, “A maximumprinciplesatisfying highorder finite volume compact WENO scheme for scalar conservation laws with applications in incompressible flows,” Journal of Scientific Computing, vol. 65, pp. 83–109, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 D. Ghosh and J. D. Baeder, “Compact reconstruction schemes with weighted ENO limiting for hyperbolic conservation laws,” SIAM Journal on Scientific Computing, vol. 34, pp. A1678–A1706, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 G. A. Sod, “A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws,” Journal of Computational Physics, vol. 27, no. 1, pp. 1–31, 1978. View at: Publisher Site  Google Scholar  MathSciNet
 P. Woodward and P. Colella, “The numerical simulation of twodimensional fluid flow with strong shocks,” Journal of Computational Physics, vol. 54, no. 1, pp. 115–173, 1984. View at: Publisher Site  Google Scholar  MathSciNet
 R. Ishii, H. Fujimoto, N. Hatta, and Y. Umeda, “Experimental and numerical analysis of circular pulse jets,” Journal of Fluid Mechanics, vol. 392, pp. 129–153, 1999. View at: Publisher Site  Google Scholar
 C. Lee and D. J. Lee, “An analysis of flow and screech tone from supersonic axisymmetric jet at the initial stage,” Journal of Mechanical Science and Technology, vol. 22, no. 4, pp. 819–826, 2008. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 MyeongHwan Ahn and DuckJoo Lee. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.