Large Eddy Simulation of Pseudo-Shock Waves Using Wall Model

Sangwook Lee^*, Yousuke Ogino, Keisuke Sawada

Department of Aerospace Engineering, Tohoku University, Sendai, Japan.

DOI: 10.4236/ojfd.2014.43024                

Abstract

A supersonic turbulent flowfield involving the pseudo-shock waves in an isolator of a supersonic combustion ramjet is computed using two different LES codes which are a high-order upwind finite volume scheme, and a sixth order compact differencing scheme utilizing the localized artificial diffusivity method for stabilizing shock waves and employing a wall model to enable the use of coarse mesh. In the validation study where a supersonic turbulent boundary layer flow over a flat plate is examined, both LES codes are well validated using velocity profile in the boundary layer given by the hot-wire anemometry and normal stress given by the laser Doppler anemometry. In particular, the sixth order compact differencing scheme gives closer agreements with these experimental data. Then, the validated LES codes are applied to solve the Mach number 2.5 supersonic turbulent flowfield involving the pseudo-shock waves. It is shown that typical features of unsteady flowfield of the pseudo-shock waves are well obtained by both schemes. Again, it is indicated that the sixth order compact differencing scheme gives closer agreements with the existing velocity data obtained by particle image velocimetry and pressure fluctuation data on the wall surface. Besides, the computational cost of the compact differencing scheme is found to be 1/7 of that for the upwind finite volume scheme, even though a wall model is solved at each grid point on the wall surface. Therefore, the obtained results in the present study allow recommending the sixth order compact differencing scheme with a wall model for simulating supersonic turbulent flowfield in an isolator involving the pseudo-shock waves.

Keywords

Large Eddy Simulation, Wall Model, Pseudo-Shock Wave

1. Introduction

For future high-speed propulsion system of aerospace vehicles, dual-mode ramjet engines have been investigated extensively in various institutions. However, according to the recent research reviews [1] [2] , there yet remain a lot of technical issues that must be resolved before putting them into practical use. One of such remaining issues is the proper control of pseudo-shock waves, known as pre-combustion shock waves that appear in the isolator of the engine. The pseudo-shock waves are a noteworthy interaction phenomenon between shock waves and boundary layer, which can alter the combustion mode from supersonic combustion to subsonic one. In most cases, the pseudo-shock waves composed of a series of oblique shock waves, X or lambda shaped shock waves, depend on the Mach number of the freestream and the boundary layer thickness. They include several shock waves, expansion fans, separation of boundary layer and slip lines. The pseudo-shock waves are known to help mix the injected fuel with oxidizer air flow, and to improve combustion efficiency. However, if the pseudo-shock waves go into the inlet due to higher back pressure caused by too much heat release in combustor, most of air cannot enter into the inlet, and the inlet then falls into an unstart condition. The aim of setting an isolator between the inlet and combustor is to contain and stabilize the pseudo-shock waves. A longer isolator can certainly contain the pseudo-shock waves stably, but it causes to increasing engine mass, frictional loss, and cooling requirement. In order to properly determine the isolator length, it is extremely important to have a capability of accurately predicting the possible length of the pseudo-shock waves at various inlet conditions.

Experiments [3] -[5] are very important to obtain the structure of the pseudo-shock waves, but it is sometimes difficult to collect detailed information about the flowfield of the pseudo-shcok waves because they emerge in supersonic flow in a duct. This shock wave is intrinsically an unsteady phenomenon and needs precise adjustment of the back pressure in the isolator. Furthermore, experimental studies are expensive because large scale facilities capable of operating over a wide range of Mach number are needed. An attempt to measure the pseudo-shock waves in flight experiment will be extremely high-cost. On the other hand, Computational Fluid Dynamics (CFD) can provide all the information of flowfield involving the pseudo-shock waves at a reasonable cost. The flowfield of the pseudo-shock waves was numerically studied using the Reynolds Averaged Navier-Stokes (RANS) equations [6] [7] . Because the RANS simulation cannot reproduce temporal evolution of turbulent structures, numerical studies have been rapidly shifted to utilize Large Eddy Simulation (LES). Although LES requires much higher computational cost than that for RANS simulation, LES is expected to provide both qualitative and quantitative information about complicated flowfield inside of an engine, particularly in the isolator where significant unsteady interaction of turbulent boundary layer and shock wave occurs.

LES calculation of the pseudo-shock waves in an isolator is not entirely new. For example, John et al. [8] carried out one such simulation. In their study, simulations of pseudo-shock waves were accomplished by utilizing a hybrid LES and RANS code. However, their computed results were not thoroughly validated through experimental data. Before utilizing LES for practical sizing of an isolator, it is indeed necessary to validate LES code using experimental data. In this study, we will first focus on the comparison of the computed mean streamwise velocity profiles and normal stress profiles with existing experimental data for a supersonic turbulent boundary layer on a flat plate. In addition, we will also attempt to justify the use of wall model in LES simulation combined with a sixth order compact differencing scheme with the localized artificial diffusivity model [9] in order to reduce computing time substantially. The computed results given by the higher-order upwind finite volume scheme employed in the previous study [10] are also compared with the experimental data. It is interesting to see the possible influence of numerical viscosity inherently involved in the computed results using an upwind scheme. Then, the pseudo-shock waves in an isolator are computed using these LES codes. Detailed comparisons are made with the existing experimental data to show that one of the present LES codes can reproduce turbulent flowfield involving the pseudo-shock waves fairly well.

In the following, we will first describe the numerical methods with the computed results of validation study in section 2, and then show the computational setup in section 3. The detailed comparisons of the computed results and experimental data are shown, and the related discussions are given in section 4. Finally, in section 5, we will give the conclusions obtained in this study.

2. Numerical Methods

2.1 Governing Equation and Subgrid-Scale Model

The governing equations are the spatially filtered mass, momentum, and energy conservation equations and the perfect gas equation of state [10] given below

(1)

(2)

(3)

(4)

where the tilde ^~ denotes the Favre filtering defined as, in which the overbar ⁻ represents the spatial filtering. The Favre filtered form of the total energy per unit mass is defined by

(5)

In above equations, is the density, (i = 1, 2, and 3) are the velocity vector components in Cartesian coordinates, is the static pressure, T is the static temperature, and is the Kronecker’s delta, R (= 287 J/Kg×K) is the gas constant of air, (= 1.4) is the ratio of the constant specific heats. These, , and are given by

(6)

(7)

(8)

where is the filtered rate of strain tensor, μ is the molecular viscosity obtained using Sutherland’s law, κ is the thermal conductivity. The thermal conductivity is obtained using the constant Prandtl number.

The SGS (Subgrid-Scale) terms, denoted with the superscript sgs in the governing equations, need to be modeled in order to close the system of equations. The SGS stress tensor can be closed using the eddy viscosity model given by

(9)

The selective mixed-scale model [11] is used to evaluate the SGS eddy viscosity. This model is categorized as a dynamic eddy viscosity model and was validated in supersonic channel flow [11] and subsonic crossflow jet [12] . The eddy viscosity is given by

(10)

where is the second invariant of the shear stress tensor, is selection function, the characteristic length scale, and the kinetic energy of the smallest-scale resolved by grid. In the present simulation, , and are used. In above equations, the subgrid-scale kinetic energy is denoted by.

The SGS total enthalpy flux vector is calculated from the gradient diffusion model given by

(11)

in which the filtered total enthalpy is defined by

(12)

and the turbulent Prandtl number of.

2.2. Numerical Schemes

In this study, two different schemes described below are assessed for conducting LES of the pseudo-shock waves in an isolator.

1) Scheme A

It is a high order MUSCL TVD upwind finite volume scheme in which a third order interpolation is employed to evaluate the dependent variables at the interface of computational cells [13] . Because of the third order interpolation, the formal spatial accuracy of the scheme becomes fourth order for 1D problems but remains second order for multidimensional problems due to lower spatial accuracy in flux integration of finite volume formulation. This scheme has been successfully applied to yield highly resolved solutions for various problems including the previous study of LES for compressible flows [10] . The numerical convective flux at the cell interface is evaluated by the Simple Low-dissipation AUSM [14] (SLAU) scheme, which is a member of AUSM type upwind scheme, while the viscous terms are evaluated by the second order central differencing scheme. The time advancement is performed by using the classical third order Runge-Kutta scheme. The non-dimensional time step, defined as, where H is the duct height, is a user specified constant time interval which is determined as 5 μs in this study, and U is the inflow velocity, is. With this definition of time step, the maximum inviscid Courant-Freiedrichs-Lewy (CFL) number becomes about 0.4.

2) Scheme B

The Navier-Stokes equations are solved using a sixth order accurate compact differencing scheme coupled with an eighth order lowpass spatial filter [15] . In order to capture shock waves stably, the localized artificial diffusivity method (LAD) [9] is employed. The fourth order Runge-Kutta method is adopted for time integration. The non-dimensional time step of is chosen to assure the CFL number to be less than 1.0. The scheme B assumes to employ an equilibrium wall model proposed by Kawai and Larsson [16] in order to substantially reduce the computational cost. The set of equations for the present wall model is given by

(13)

(14)

where the eddy viscosity is evaluated based on the mixing length model given by

(15)

in which is the distance from the wall in wall units, the model parameter and. The set of equations for the wall model are solved separately on a different grid having of 43 grid points in y direction that covers from the wall to the inner portion of the boundary layer typically up to about. At the outer boundary of the wall grid referred to as the matching location, the instantaneous temperature and velocity components are provided from the LES results as the boundary conditions, while an adiabatic wall and nonslip conditions are assumed on the wall surface. A second order central differencing scheme is used to discretize the velocity and temperature derivatives. The discretized equations for the wall model are implicitly integrated using a tri-diagonal matrix solver. Using the obtained solution for the wall model, the wall friction and heat flux on the wall surface are provided as the boundary condition for LES computation. The matching location is assumed at the 5th grid points from the wall where in the present study. According to Kawai and Larsson [16] , this treatment enables us to avoid so called log-layer mismatching.

2.3. Validation of LES Schemes

A supersonic turbulent boundary layer on a flat plate is one of the fundamental compressible flowfields suitable for validation of LES codes. The simulated flow conditions are summarized in Table 1. The turbulent inflow is given by the rescaling method which will be mentioned in section 3. The computed results given by scheme A and scheme B are compared with the existing experimental data such as the velocity profile given by the hot- wire anemometry (HWA) and the normal stress given by the laser Doppler anemometry (LDA) appeared in [17] . The existing LES result of Sagaut et al. [17] is also compared with the present LES results.

Three different computational grids referred to as coarse (113 × 91 × 45), medium (159 × 121 × 61), and fine (225 × 170 × 89) are employed for scheme A in the present validation study, while another mesh summarized in Table 2 is employed for scheme B. The minimum mesh interval at the wall surface for scheme A computation is less than for all three meshes. On the other hand, the minimum mesh interval of the computational mesh for scheme B calculation is as large as. Because the wall model is solved simultaneously in scheme B, a use of such coarse mesh is possible without deteriorating the computed results as shown below. The wall grid for scheme B employs 43 grid points in y direction assigned at each grid point on the wall surface. The minimum grid spacing of each wall grid is less than.

Figure 1 shows the computed profiles in wall units, where is Van Driest-transformed mean streamwise velocity, with the experimental data obtained by HWA and also LES result of Sagaut et al. [17] . One can immediately find that the computed profiles given by scheme A, even for the case of using the fine mesh, become larger than those experimental data and that given by LES of Sagaut et al. On the other hand, a fair agreement is obtained for the profile given by scheme B with those experimental data and that by LES, although the number of grid points of the computational mesh for scheme B is 159 × 132 × 63, which is similar to that for the medium grid for scheme A.

In Figure 2, the computed normal stress profiles are compared with LDA experimental data and the LES result of Sagaut et al. [17] , where is the instantaneous streamwise velocity fluctuation and is the wall shear stress. The computed normal stress profiles given by scheme A have larger peak values at than

Table 1. Simulated flow conditions.

Table 2. Grid parameters for scheme B.

Figure 1. Comparison of profiles plotted in wall units.

Figure 2. Comparison of normal stress profiles.

those given by LDA data and also LES result of Sagaut et al. In addition, the normal stress profiles show a plateau in between 0.25δ and 0.75δ. On the other hand, the normal stress between 0.1δ and 0.25δ given by scheme B agrees well with LDA data, though a plateau also appears in between 0.25δ and 0.75δ. Similar tendency was also observed in the LES result of Sagaut et al. in between 0.1δ and 0.5δ.

It should be noted that scheme A does not attain a grid converged solution using the fine mesh. Among several possible reasons, this can be attributed most likely to the use of TVD limiter and use of upwind numerical flux in scheme A. The LES code attempts to resolve vortices of grid size, but the TVD limiter will inevitably limit the velocity peak of small vortices to yield velocity jump of order of grid size at the cell interface. This jump then results in to introduce numerical viscosity through the use of an upwind scheme. The larger value given by scheme A is probably caused by the above reason. On the other hand, scheme B gives fair agreements with the experimental data of both HWA and LDA, and the LES result of Sagaut et al. Based on these results, although we do not check the grid convergence yet, scheme B is judged validated.

The use of wall model in scheme B certainly enables us to employ a larger time step due to use of coarse mesh in near wall region. For the present supersonic turbulent flow over a flat plate, the computational cost to carry out LES using scheme A with the medium mesh is found seven times higher than that using scheme B.

3. Computational Setup

In the works of Choi et al. [18] and Hoshino et al. [19] , the pseudo-shock waves were measured in a duct having a square cross section whose height H was 30 mm and the length of the duct was 290 mm (9.67 H). The Mach number M_a was 2.5, the total pressure p₀ was 100 kPa, and the total temperature T₀ was 290 K at the inlet of the duct. The Reynolds number based on the duct height was 3 × 10⁵. As shown in Figure 3, the computational meshes (Mesh A and Mesh B) discretize the entire duct region. Table 3 shows the parameters of these computational meshes. Mesh A is employed for scheme A, and mesh B for scheme B. The minimum mesh interval of mesh A satisfies both on the wall surface, while that of mesh B is. The wall grid for Mesh B has 43 grid points at each grid point on the wall surface where the adiabatic and no slip conditions are assumed. According to experimental data, the back pressure at the outlet boundary is chosen as 0.27 p₀ in order to generate the pseudo-shock waves in the computational domain.

The turbulent inflow condition at the inlet is given by the multiwall rescaling method proposed by Boles et al. [20] which is an extension of the rescaling method introduced by Sagaut et al. [17] . The multiwall rescaling method is conducted in the domain enclosed by the dotted lines in Figure 3. As shown in the next section, the inflow exhibits large scale turbulent bulges and low speed streaks which are typical features of turbulent boundary layer flows. The validity of the given turbulent inflow is examined in terms of the boundary layer thickness and velocity profiles at the rescaling point chosen at 12δ form the inlet boundary.

4. Results and Discussions

In Figure 4, the instantaneous distribution of the streamwise velocity component computed by scheme B are plotted on the center plane, the exit plane and the horizontal plane which is slightly lifted from the bottom wall. One can find turbulent bulges and low speed streaks in between the inlet and the first shock wave. In addition, the x-shaped first shock wave on the center plane can be clearly seen. In the downstream of the first shock wave, there appear the second and third shock waves. Those shock waves are the characteristic feature of the pseudo-shock waves. The interaction of these shock waves with the unsteady turbulent boundary layer forms complicated flow structures involving separated flow regions that appear on the wall surface as well as at the corner region.

Figure 5 shows a comparison of the schlieren image taken at an arbitrary instant in the experiment [19] (top) with the computed images by scheme A (middle) using the medium mesh and that by scheme B (bottom). One can say that these schlieren images resemble each other. In particular, the intervals of the shock waves look quite similar. It seems, however, that the turbulent boundary layers in front of the foot of the first shock wave computed by both scheme A and scheme B become slightly thicker than that in the experiment. A detailed comparison of the computed schlieren images indicates that the first shock wave in near wall region given by scheme B is more clearly visible than that given by scheme A. Therefore, it can be said that the present compact scheme utilizing the local artificial diffusivity model (scheme B) can capture shock waves quite well in the complicated turbulent flowfield.

Figure 6 compares the mean streamwise velocity component plotted on several cross sections. The top figure shows those obtained by particle image velocimetry (PIV) [18] , while the middle figure shows those given by scheme A and the bottom figure by scheme B. The horizontal axis is given by, where the origin

Table 3. Grid parameters.

Figure 3. Schematic illustration of computational domain.

Figure 4. Instantaneous distribution of the streamwise velocity component.

Figure 5. Comparison of schlieren images in the centerplane. Experiment (Top), scheme A (middle), and scheme B (bottom).

Figure 6. Comparison of streamwise mean velocity distributions plotted in cross sections at = 0.0, 0.33, 0.67, 1.0, 1.33, 1.67, 2.0, 2.33, and 2.67.

of is determined at the foot of the first shock wave, and is non-dimensionalized by the duct height H. because the position of the first shock wave does not necessarily agree each other. The mean velocity distribution of PIV data is formed from 100 image pairs sampled at a rate of 15 Hz, while those given by scheme A and scheme B are formed from 2000 instantaneous field data. From these figures, one can clearly see the characteristic feature of the pseudo-shock waves that repeated deceleration and acceleration appears. At and 2.33, the streamwise velocity in the core region becomes larger in the PIV data while the computed results, in particular those given by scheme B, also show the same tendency. This is of course a natural consequence of the reasonable agreement in shock intervals indicated in Figure 5. On the other hand, behind the first shock wave at and 1.0, low speed regions appear at the middle of the upper and lower wall boundaries, while in further downstream, the low speed region rather appears at the corner. This tendency also agrees with the experimental data.

Figure 7 shows the velocity fluctuation intensities in the cross sections at and. In this figure, those shown in the first row are the experimental data [18] , the second row the computed results of

Figure 7. The distributions of velocity fluctuation intensity for x, y, and z components in cross sections at and 2.5.

scheme A and the third row the scheme B, while those in the first column are the intensity of velocity fluctuation of x component, the second column the y component and the third column the z component. One can find that the computed velocity fluctuation intensities in the core region become weaker than the surrounding region. This tendency agrees well with the experimental data. The high intensities in x component of velocity fluctuation at four corners of the square cross section at appear only in the experimental data. Such high intensities of velocity fluctuation, however, are absent in y and z components. What flow structure can give such high intensities only in x component is yet unknown, Although one can say that the turbulence is obviously anisotropic at these locations. It should be noted that the PIV data is sampled over a few seconds while the computed results are sampled over only 10ms which corresponds to one period of shock excursion. Therefore, in the experiment, the shock wave can have sufficient chances to move across which could induce the unknown flow structure, while the computed period was too short for the shock wave to come across that position. For y and z components, the patterns of higher intensity region in the computed results show qualitative agreement with the experimental data, though the computed intensities in the surrounding region become higher than those obtained in the experiment.

The wall pressure distributions obtained in the experiment [19] , that computed by scheme A and by scheme B are shown in Figure 8. The computed pressures are averaged over 10 ms. The onset location of pressure rise on the wall surface was determined in the experiment. The computed onset locations on the wall surface given by scheme A and scheme B agree well with the experiment. In particular, scheme B gives the onset location and successive pressure rise fairly close to those obtained in the experiment. One can find that significant pressure oscillation occurs along the centerline of the duct in the computed result of scheme B. This of course corresponds to the repeated flow deceleration and acceleration which is typical feature of the pseudo- shock waves.

In Figure 9, the wall pressure fluctuations at, 6.3, and 8.0 obtained in the experiment are compared with those computed fluctuations given by scheme A and scheme B at the same positions. At the upstream side of the first shock wave, the amplitude of pressure fluctuations are all less than 0.01, though the amplitude given by scheme A seems slightly larger than others. At the foot of the first shock wave, the low frequency components begin to appear due to unsteady shock motions. At, the low frequency components seem to diminish at this position and high frequency fluctuations become dominant again both in the experimental data and computed results. The amplitude of fluctuation given by scheme A is twice larger than others at this position. The corresponding power spectrum densities of pressure fluctuations computed by the fast Fourier transform at the chosen positions are shown in Figure 10, where G is the power spectrum and f is the frequency. At the upstream side of the first shock wave, only the high frequency components (>10 kHz) due to turbulence motions are dominant both in the experimental data and those computed results. On the other hand, at the foot of the shock wave, lower frequency components

Figure 8. Comparison of wall pressure profiles.

Figure 9. Comparison of wall pressure fluctuations in the center plane at, 6.3, and 8.0.

Figure 10. Comparison of power spectrum density of pressure fluctuations at, 6.3, and 8.0.

(≈300 Hz) do appear as observed in Figure 9. In addition, high frequency components remain only in the computed pressure fluctuations. At, those frequency components below 3 kHz remain in the experiment, while they are almost diminished in the computed pressure fluctuations. This implies that some unsteady shock motions are yet significant in the experiment. Unfortunately, the cause for absence of low frequency components in the computed results at this position is not identified yet.

Finally, the root mean square (RMS) of pressure fluctuations along the centerline of the bottom wall is shown in Figure 11. One can find that these RMS rapidly increase at the foot of the first shock wave due to significant unsteady shock motions. Although the RMS of the computed pressure fluctuations upstream of the first shock wave becomes twice larger than that of the experimental data due to thicker turbulent boundary layer in

Figure 11. RMS profiles of pressure fluctuations.

the computed result as indicated in Figure 5, the obtained RMS by scheme B reasonably agrees with the experimental data behind the first shock wave, while that given by scheme A obviously overestimates the fluctuation level.

The computed results indicate that both scheme A and scheme B can obtain typical features of the pseudo- shock waves quite well. In particular, scheme B gives the computed results that are close to experimental data of velocity and pressure fluctuations, in addition to the overall unsteady flowfield where significant interaction of shock waves and turbulent boundary layer occurs. Because the computational cost of scheme B is substantially lower than that of scheme A, LES simulation of the pseudo-shock waves in an isolator now seems feasible. Future works will aim to enhance mixing of fuel and oxidizer air flow in combustion chamber utilizing the pseudo- shock waves in an isolator, and these will be presented elsewhere.

5. Conclusions

In the present study, a supersonic turbulent flowfield in an isolator involving the pseudo-shock waves is computed using two different LES codes. Scheme A is a high-order upwind finite volume scheme. Scheme B is a sixth order compact differencing scheme utilizing localized artificial diffusivity method to capture shock waves stably and employing a wall model to enable the use of coarse mesh at the wall surface. Followings are obtained in the course of the present study:

1) A supersonic turbulent boundary layer flow over a flat plate is computed by two LES codes (scheme A and scheme B). By comparing with the existing experimental data and the LES results, these two codes are successfully validated. In particular, closer agreements with the experimental data and also the LES results are shown when scheme B is employed in LES.

2) The validated LES codes are applied to solve the pseudo-shock waves in an isolator. Typical features of unsteady flowfield of the pseudo-shock waves are well obtained by both scheme A and scheme B. Again, it is shown that scheme B gives closer agreements with the existing experimental data of the pseudo-shock waves.

3) The computational cost of scheme B is found to be 1/7 of that of scheme A, even though the wall model is solved at each grid point on the wall surface.

4) Through comparisons, the obtained results in the present study allow recommending the compact differencing scheme with a wall model (scheme B) for simulating supersonic turbulent flowfield in an isolator involving the pseudo-shock waves.

Acknowledgements

The first author gratefully thank to Dr. Byuong-il Choi and Prof. Goro Masuya for providing experimental data and also for stimulating discussions. All the computations are conducted by SX-9 supercomputer at the Cyberscience Center, Tohoku University.

Nomenclature

E : Total energy

f : Frequency

G : Power spectrum

H : Duct height

H' : Total enthalpy

M_a : Mach number

p : Pressure

q : Heat flux

Re : Reynolds number

T : Temperature

t : Time

u : x component of velocity vector

v : y component of velocity vector

x : Streamwise coordinate

y : Vertical coordinate

z : Spanwise coordinate

μ : Molecular viscosity

δ : Boundary layer thickness

ρ : Density

κ : Thermal conductivity

τ : Shear stress

σ : Viscous stress

Subscripts

ave : Average value

match : Matching location

rms : Root mean square

w : Wall value

∞ : Free stream

Superscripts

sgs : Sub grid scale

+ : Wall unit quantity

’ : Fluctuation

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Ronald, S.F. (2004) A Century of Ramjet Propulsion Technology Evolution. Journal of Propulsion and Power, 20, 27-58. http://dx.doi.org/10.2514/1.9178
[2]	Bouchez, M., Perillat, V., Avrashkov, V. and Kopchenov, V. (2011) Numerical and Experimental Scientific Investigation of Combustion in a Translating Cowl Dual-Mode Ramjet. 17th AIAA international Space Planes and Hypersonic Systems and Technologies Conference, San Francisco, 11-14 April 2011, 2011-2313
[3]	Tamaki, T., Tomita, Y. and Yamane, R. (1970) A Study of Pseudo-Shock: 1st Report, λ-Type Pseudo-Shock. Japan Society of Mechanical Engineers, 13, 51-58. http://dx.doi.org/10.1299/jsme1958.13.51
[4]	Waltrup, P.J. and Billig, F.S. (1973) Structure of Shock Wave in Cylindrical Duct. AIAA Journal, 11, 1404-1408. http://dx.doi.org/10.2514/3.50600
[5]	Sullins. G. and McLafferty, G. (1992) Experimental Results of Shock Trains in Rectangular Ducts. AIAA 4th International Aerospace Planes Conference, Orlando, 1-4 December 1992, 92-5103.
[6]	Sun, L.Q., Sugiyama, H., Mizobata, K. and Fukuda, K. (2003) Numerical and Experimental Investigations on the Mach 2 Pseudo-Shock Wave in a Square Duct. Journal of Visualization, 6, 363-370. http://dx.doi.org/10.1007/BF03181743
[7]	Kawatsu, K., Koike, S., Kumasaka, T., Masuya, G. and Takita, K. (2005) Pseudo-Shock Wave Produced by Backpressure in Straight and Diverging Rectangular Ducts. AIAA/CIRA 13th International Space Planes and Hypersonic Systems and Technologies, Capua, 16-20 May 2005, 2005-3285.
[8]	John, A.B., Hagenmaier, M.A. and Hsuet, K.Y. (2011) Analysis of Hybrid LES/RANS Simulations of a Back Pressured Supersonic Isolator. 47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibition, San Diego, 31 July-3 August 2011, 2011-5825.
[9]	Kawai, S. and Lele, S.K. (2008) Localized Artificial Diffusivity Scheme for Discontinuity Capturing on Curvilinear meshes. Journal of Computational Physics, 227, 9498-9526. http://dx.doi.org/10.1016/j.jcp.2008.06.034
[10]	Watanabe, J., Kouchi, T., Takita, K. and Masuya, G. (2011) Numerical Study on Turbulent Structure of Transverse Jet into Supersonic Flow. AIAA Journal, 49, 2057-2067. http://dx.doi.org/10.2514/1.J051067
[11]	Lenormand, E., Sagaut, P. and Ta Phuoc, L. (2000) Large Eddy Simulation of Subsonic and Supersonic Channel Flow at Moderate Reynolds Number. International Journal for Numerical Methods in Fluids, 32, 369-406. http://dx.doi.org/10.1002/(SICI)1097-0363(20000229)32:4<369::AID-FLD943>3.0.CO;2-6
[12]	Dandois, J., Garnier, E. and Sagaut, P. (2006) Unsteady Simulation of a Synthetic Jet in a Crossflow. AIAA Journal, 44, 225-238. http://dx.doi.org/10.2514/1.13462
[13]	Yamamoto, S. and Daiguji, H. (1993) Higher-Order—Accurate Upwind Schemes for Solving the Compressible Euler and Navier-Stokes Equations. Computers & Fluids, 22, 259-270. http://dx.doi.org/10.1016/0045-7930(93)90058-H
[14]	Shima, E. and Kitamura, K. (2009) On New Simple Low-Dissipation Scheme of AUSM-Family for All Speeds. 47th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Orlando, 5-8 January 2009, 2009-136.
[15]	Lele, S.K. (1992) Compact Finite Difference Schemes with Spectral-Like Resolution. Journal of Computational Physics, 103, 16-42. http://dx.doi.org/10.1016/0021-9991(92)90324-R
[16]	Kawai, S. and Larsson, J. (2012) Wall-Modeling in Large Eddy Simulation: Length scales, Grid Resolution, and Accuracy. Physics of Fluids, 24, 015105. http://dx.doi.org/10.1063/1.3678331
[17]	Sagaut, P., Garnier, E., Tromeur, E., Larcheveque, L. and Labourasse, E. (2004) Turbulent Inflow Conditions for Large-Eddy Simulation of Compressible Wall-Bounded Flows. AIAA Journal, 42, 469-477. http://dx.doi.org/10.2514/1.3461
[18]	Choi, B., Takae, T., Kouchi, T. and Masuya, G. (2012) Mean and Turbullent Velocity Field for Jet Injected into Supersonic Flow with Pseudo-Shock Wave. Journal of Propulsion and Power, 28, 971-981. http://dx.doi.org/10.2514/1.B34371
[19]	Hoshino, T., Kouchi, T. and Masuya, G. (2010) Stability of Pseudo-Shock Waves in a Straight Rectangular Duct installing Corner-Bumps. Proceedings of the 54th Space Sciences and Technology Conference, Shizuoka, 17-19 November 2010, JSASS Paper, 2010-4146. (in Japanese)
[20]	Boles, J.A., Choi, J.I., Edwards, J.R. and Baurle, R.A. (2010) Multi-Wall Recycling/Rescaling Method for Inflow Turbulence Generation. 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Orlando, 4-7 January 2010, 2010-1699.