Abstract
The fifth-order monotonicity-preserving (MP5) scheme is an accurate and low dissipative numerical method. As a finite-volume method, MP5 adopts the Roe-flux scheme for solving the numerical flux in the compressible Euler equation. However, due to the deficiency of the MP limiter and Roe-flux 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 high-accuracy of MP5, we propose a hybrid flux method: the Roe-flux is used in the global computational domain, but the first-order Lax-Friedrich (LF)-flux is adopted only for trouble grids. The numerical results of shock-tube 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 two-dimensional 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 blow-up of calculation even in the extreme jet flow condition.
1. Introduction
The usage of an efficient high-order shock-capturing 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 well-known essentially nonoscillatory (ENO) [1] and weighted ENO (WENO) schemes [2] are very robust high-order method to capture the shock-wave without introducing any nonphysical oscillation; however, they have been proven as a large numerical-dissipative 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, fifth-order monotonicity-preserving (MP5) scheme [5] was proposed by Suresh and Huynh. In this method, the Roe-flux scheme [6] was adopted as the basis method for solving numerical flux of compressible Euler equation. The well-designed 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 shock-wave and small-scale 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 shock-sensor for an efficient adaptation of an MP limiter. In addition, although not related to accuracy, Daru and Tenaud [13] proposed a one-step scheme by reinterpreting the MP limiter as a TVD-like 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) Lax-Friedrichs (LF)-flux [16] was used for the basis flux method because it has a positivity-preserving capability (Roe-flux used in the MP5 scheme has no positivity-preserving capability [17]). They successfully showed the robustness of their method (MP5-R) 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 LF-flux 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 LF-flux 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 high-accuracy. The key strategy of present scheme is the hybridization of Roe- and the first-order LF-fluxes depending on local density and pressure state. Initially, the solution in the whole flow-field is obtained using the original MP5 scheme; however, only for the local trouble region detected by criteria, the first-order LF-flux is applied to modify the numerical flux. The first-order never introduces the nonphysical oscillation; therefore, the additional TVD condition used in MP5-R 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 shock-tube 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 MP5-R and MP5-Hyb 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 high-accuracy 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 blow-up 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 one-dimensional linear advection equation with speed of wave
For the moment, assume that positive wave speed . Let be grid size and be the index of cell-centered grid. After integrating (1) over the grid , it yieldsor
For accurate calculations, the high-order interpolation for the numerical flux is performed. In this study, we consider the monotonicity-preserving (MP) scheme framework as introduced in the next section.
2.2. Monotonicity-Preserving (MP) Scheme Framework
The procedure of MP scheme is separated in two parts: (i) high-order interpolation and (ii) application of MP limiter (or constraints) for shock-capturing. The former procedure enhances the order of accuracy by obtaining the cell-interface flux using neighboring grids. Let us consider the cell-averaged variable on an arbitrary grid as , then the fifth-order upwind interpolation adopted in the MP scheme is given by [5],where , and . Equation (4) is expressed for the right-going 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 shock-capturing 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 monotonicity-preserving 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 spatial-accuracy can be obtained by combining (4) with the spatial-derivative 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 one-dimensional 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 well-known Roe- [6] or Lax-Friedrich (LF)-flux [16] to obtain a solution. In general, the procedure of high-order 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 right-eigenvector 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. Roe-Flux in Original MP5 Scheme for High-Accuracy
The original fifth-order monotonicity-preserving (MP5) scheme proposed by Suresh and Huynh [5] was developed based on the Roe-flux [6] which is the approximated version of Godunov flux scheme. The merit of Roe-flux is the high-accuracy and good resolution of shock waves because it considers not only the direction of characteristic wave propagation, but also the waves themselves at the cell-interface [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 cell-interface to obtain and . (iii) Return to the physical space by multiplying the right-eigenvector: and . (iv) Left and right directional conservative variables, and , are used in the Roe-flux asThe differential of the characteristic vector is denoted by as follows:
The notation stands for Roe-averaged 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 cell-interface ( 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. Roe-flux scheme has no positivity-preserving capability [17]. It has a potential to give negative cell-averaged density and pressure ( and ) even though slightly positive cell-interface values ( and ) are used in the Roe-flux.
3.2. Lax-Friedrich Flux in MP5-R Scheme for Robustness
To overcome the aforementioned limitation of original scheme, MP5-R [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 Roe-flux was substituted with the LF-flux scheme which has positivity-preserving capability [19]. As a flux splitting method, the LF-flux is monotone itself, which means never introducing nonphysical oscillation so that density and pressure values are always positive. However, the LF-flux is quite dissipative than Roe-flux [18, 19]. Therefore, the overall accuracy of MP5-R scheme is degenerated from that of the original MP5 scheme.
The MP5-R scheme based on the LF-flux can be summarized as follows: (i) Calculate the characteristic variable shown in (21) on stencils for j+1/2 flux asHere, is obtained using LF-flux as (ii) Perform the procedure of MP scheme shown in (4)~(9) on each vector of in left and right direction at the cell-interface to obtain and . (iii) Return to the physical space by multiplying the right-eigenvector: and . (iv) Then, the final numerical flux can be obtained as
4. Hybrid Roe- and First-Order LF-Flux for High-Accuracy and Robustness
As introduced in the previous section, MP5 and MP5-R schemes have different numerical characteristics by adopting a different flux scheme. The former is quite accurate due to using Roe-flux; 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 LF-flux.
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 LF-flux only for trouble grids to modify the negative density and pressure; otherwise, apply the Roe-flux to obtain accurate solution. By implementing such a hybrid-concept, we try to achieve both of the respective advantages of the MP5 and MP5-R 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 flux-in/out on the j-th 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 Roe-flux.where is the total number of grids. After that, regardless of local flow state, just perform a time-marching to obtain numerical solutions in the next time step as . In this step, if the entire flow-field 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+1-th 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 ill-posed. Furthermore, in practice, negative density or pressure at the cell-interface ( and ) can give a NaN (Not a Number) during the Roe-flux calculation (see, Roe-averaged 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 n-th time step and calculate residue only for the trouble grids. As mentioned before, we try to use the LF-flux to obtain new flux values in this step. However, unlike the MP5-R scheme, we adopt the first-order of LF-flux which does not require the high-order procedures as described in Section 2.2. The cell-averaged 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 first-order LF-flux, the solution can be expressed as
We can see that is determined by three set of vectors and . Under the condition of Courant-Friedrichs-Lewy (CFL) number <1.0, the first term on the right-hand-side 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 MP5-R 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 shock-wave/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 high-accuracy 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 time-integration. Henceforward, MP5-Hyb (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 hybrid-factor 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 Runge-Kutta method is used for the time-integration. In this case, the modification procedures for negative density or pressure should be applied in every internal stages of Runge-Kutta method. Figure 3 summarizes the procedures of the hybrid flux method with k-th stages Runge-Kutta 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 two-dimensional 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 one-dimensional cases and supersonic jet flow simulations in asymmetrical two-dimensional domain are performed, which is the main application of the present scheme. For all the following computations, a third-order TVD Runge-Kutta time-integration 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 Runge-Kutta time-integration; 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 MP5-Hyb 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 MP5-R and MP5-Hyb schemes have the almost the same slope as five (slope indicates the order of accuracy); however, the absolute error values of MP5-Hyb is less than MP5-R scheme in all different number of grids. Such results are acceptable because MP5-Hyb scheme includes Roe-flux whereas MP5-R takes only dissipative LF-flux scheme. Details of error and its slope are listed in Table 1.
5.2. Shock-Tube Problems
Shock-tube problem demonstrates the dynamics of two different state of gas generating a shock-wave, contact-discontinuity, 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 Shock-Tube Problem (Normal Condition)
This is a normal condition case of shock-tube problem. On the domain , the initial condition was given as [24]The number of grids was 80 and a zero-gradient 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 MP5-R and MP5-Hyb are shown in Figure 5. As a result, both schemes gave a good trend of solutions; however, at the contact-discontinuity (in the enlarged figure) the result of MP5-Hyb is slightly close to the analytical solution than the MP5-R. Aforementioned, Roe-flux 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 MP5-R and MP5-Hyb is similar in coarse and fine grids. The gradient of solution can be steepened with a fine grid. Since the Sod shock-tube problem does not involve low density and pressure values, the local modification by the first-order LF-flux was not activated in the MP5-Hyb scheme.
(a)
(b)
5.2.2. Leblanc Shock-Tube Problem (Extreme Condition)
The LeBlanc problem [20] is the extreme case of the shock-tube 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 zero-gradient 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, y-axis 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 MP5-R and MP5-Hyb schemes can solve this problem without blows-up of calculation. On the other hand, the accuracy between MP5-R and MP5-Hyb 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 MP5-Hyb scheme. In addition, the location of shock-wave presented at x ~ 8 in density and pressure plots is more precise in the case of MP5-Hyb. Figure 8 summarizes the changes of density solutions with different number of grids. In both cases of coarse and fine grids, MP5-Hyb shows the better location of shock-wave than MP5-R. However, even with the finer grid, the solution of MP5-R is not much improved. These results confirm that MP5-Hyb scheme can be an efficient numerical method to solve shock-wave 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. Shock-Entropy Wave Interaction Problem
The simulation for shock-entropy wave interaction problem [1] was carried out in this subsection. From this test we can assess the capability of shock-capturing and resolving accuracy on a high-frequency local extrema region. On the domain , the initial condition was given as follows:
A shock-wave of Mach=3.0 interacts with the linear entropy waves. A zero-gradient 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) N-shape discontinuities and (b) high-frequency region. On the region (a), the solution of MP5-R scheme is overestimated at the peaks of N-shape whereas MP5-Hyb scheme presents a clear solution fit on the reference solution. In particular, on the region (b), we can confirm that the accuracy of MP5-Hyb is obviously improved compared to MP5-R 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 shock-tube tests, the accuracy differences are dramatic because the solution in a high-frequency 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 high-frequency region, and we can still observe the accuracy differences between MP5-R and MP5-Hyb.
5.3.2. Two-Blast-Wave Interaction Problem
We considered the two-blast-wave interaction problem [25]. In particular, this test includes very low density and pressure state in the shock-shock 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 MP5-R and MP5-Hyb schemes give reasonable results without blows-up of calculation while the MP5 scheme does not maintain its calculation at t=0.027 (the results of blows-up moment are presented in plots). In particular, the solution of MP5-Hyb is much accurate than MP5-R 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 MP5-Hyb presents slightly oscillating solution near x=0.71, but the peaks of local extrema can be much more accurately predicted than MP5-R. With a fine grid, the maximum values of extrema can be predicted by MP5-Hyb whereas grids seem still insufficient for MP5-R. 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 MP5-R.
(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 shock-wave and shear layer to be solved using the method of high-accuracy for a good quality of solution, and (2) low density and pressure appear at the core of vortices induced by high-speed 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 high-speed supersonic jet case which generates the vortices with very low density and pressure. The governing equation is the two-dimensional 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 two-dimensional 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 MP5-R and MP5-Hyb schemes.
(a)
(b)
Figure 16 shows the comparison of instantaneous density contours obtained by MP5-R and MP5-Hyb 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 MP5-Hyb scheme while they are smeared due to large numerical dissipation in the case of MP5-R scheme. These facts confirm that the Roe-flux scheme in the hybrid flux method increases the resolving performance for small-scale flow structures.
(a)
(b)
6.2. Mach 3.0 Supersonic Hot Jet Simulation
The present scheme was applied for a high-Mach number and high-temperature 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 x-direction 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 MP5-Hyb schemes at very initial time. At ( is a nondimensional time), the calculation of MP5 scheme is blow-up at the core of initial vortex due to negative pressure caused by the fast rotating flow. On the other hand, MP5-Hyb can overcome this situation, and maintains the calculation. To investigate another blow-up 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 flow-field of MP5-Hyb 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 blow-up 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 MP-Hyb is still stable. Figure 19 presents density and pressure contours of fully developed jet flow obtained by the MP5-Hyb scheme (at ). Such flow-field can be obtained by overcoming lots of blow-up situations; therefore, we confirm that MP5-Hyb 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, MP5-Hyb 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 Roe-flux is an accurate shock-capturing scheme; however, the calculation could fail when very low density/pressure appears. MP5-R scheme was proposed to overcome this limitation by adopting the LF-flux and additional TVD condition. However, the LF-flux scheme is more dissipative than the Roe-flux scheme; therefore, the overall accuracy is degenerated from the original MP5.
To overcome such limitations but still achieve a high-accuracy of MP5, we have proposed a hybrid flux method. The strategy is that the Roe-flux is used in the global computational domain, but the first-order LF-flux is adopted only for trouble grids resulting in negative density or pressure. The first-order LF never introduces numerical oscillation, and it has a positivity-preserving capability itself; therefore, the additional treatment such as TVD condition is not necessary.
The results of shock-tube and complicated interaction problems proved that the MP5-Hyb scheme has the robustness without resulting in calculation failure. In addition, it is beneficial to capture the shock waves and local extrema. The MP5-Hyb scheme was applied to simulate supersonic jet flows in different conditions. As a result, the MP5-Hyb scheme has better accuracy than the MP5-R scheme for small-scale structures of vortices. In addition, the MP5-Hyb scheme can maintain a stable calculation even with the extreme jet flow condition whereas the MP5 scheme failed. Therefore, we conclude that MP5-Hyb 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 right-eigenvector 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 NRF-2014M1A3A3A020348. This work was also conducted at High-Speed Vehicle Research Center of KAIST with the support of the Defense Acquisition Program Administration and the Agency for Defense Development under Contract UD170018CD.