1 Introduction

The prediction of laminar-turbulent transition is of great importance for the design of long distance flying vehicles ^[1], and it has been and still remains a subject of extensive investigation ^[2-4].

In designing flying vehicles, usually, experiment plays a very important role. Therefore, people expect that experiment can also be very helpful in the prediction of laminar-turbulent transition. For this purpose, quiet wind tunnels were designed and are being used. Generally speaking, to simulate the real flight, certain similarity rules have to be satisfied, as shown in Ref. [5]. However, to simulate the laminar-turbulent transition of the boundary layer of the super/hypersonic flying vehicles, more stringent requirements are necessary. For example, first, transition location is sensitive to viscosity, which depends on temperature. Thus, the temperature distribution of the flow fields must be the same for both the real flying vehicle and the model. Second, transition location is also sensitive to the free-stream disturbances. Thus, the disturbances in the free stream in their non-dimensional form must be the same for both the experiment and the real flying vehicle. Obviously, these two requirements are impossible to be realized in practice.

The other way to model the transition process is to do direct numerical simulation (DNS), if computer resource allows. In this case, proper numerical scheme has to be used. One criterion for its choice is how well it can treat the interaction of the disturbances with the shock, because, for super/hypersonic flight, disturbances in the oncoming flow would certainly first interact with the shock, and then disturbances behind the shock resulting from the interaction would interact with the boundary layer to trigger instability waves in the boundary layer. Therefore, whether the numerical scheme can faithfully reproduce the interaction is crucial for its success. Usually, the shock-fitting method can produce reliable results ^[6-10], but it can hardly be used for three-dimensional problems ^[11]. There have been people using the shock-capturing scheme for this purpose ^[12-14]. However, they assumed implicitly that the shock-capturing scheme can faithfully reproduce the disturbances behind the shock, but with no proof or convincing discussion to confirm its reliability. Therefore, how well the shock-capturing scheme can reproduce the interaction of small amplitude disturbances with the shock is a problem worth being investigated, as pointed out in the latest review article "Two problems in the transition and turbulence for the near space hypersonic flying vehicles" by Zhou and Zhang ^[15].

Therefore, in this paper, our main concern is the accuracy of shock-capturing method in its application to the interaction of small disturbance and the shock.

In a linear regime, a general small amplitude free-stream disturbance in a uniform flow can be decomposed of three elementary types of components ^[16], namely, acoustic wave, vorticity wave, and entropy wave. The acoustic wave propagates at the speed of sound relative to the moving fluids, while the vorticity and entropy waves are convected passively by the mean flow. Furthermore, the acoustic wave can also be distinguished as fast or slow acoustic wave according to whether its group velocity is larger or smaller than the moving velocity of the flow. All types of disturbances may be generated behind the shock no matter what type of disturbance wave interacts with the shock ^[17].

However, in the high altitude atmosphere, the velocity and temperature are bound to be non-uniform, which are the cause of entropy wave and vorticity wave, but there can hardly be external object generating acoustic waves, while the generation of sound wave by the motion of the fluid itself, is far less effective. Therefore, in this paper, we only pay attention to the incidence of entropy and vorticity waves. Besides, for super/hypersonic flying vehicles, the free stream disturbances can always be regarded as small amplitude disturbances.

Thus, the problem we are dealing with is, when small amplitude free-stream disturbances meet a stationary oblique shock, what kind of disturbances would be induced downstream of the shock, and whether the numerical simulation using the shock-capturing scheme can obtain the results with sufficient accuracy. For this purpose, we need reliable theoretical or experimental results to compare with numerical results. Reference ^[17] can serve for this purpose. The paper is organized as follows. In Section 2, we use the method of McKenzie and Westphal ^[17] and the errata given by Anyiwo and Bushnell ^[18] to analyze the interaction of small amplitude disturbances with an oblique shock. In Section 3, we perform numerical simulations using the shock-capturing scheme with moderate order of accuracy for the same problem. In Section 4, numerical results are compared with the theoretical results to check if the shock-capturing method can produce sufficiently accurate results. In Section 5, we conclude the paper with a summary discussion.

2 Theoretical analysis

The method of McKenzie and Westphal ^[17] is adopted to formulate the problem. For the readers' convenience, we outline the main procedures. In addition to their work, in Subsection 2.2, we complement another physical condition to remove the nonphysical solution that may arise for the solutions of wave vectors. Furthermore, we deduce the condition to determine whether fast or slow acoustic wave can be generated behind the shock, which is also in Subsection 2.2.

2.1 Problem formulation

The coordinate system (x, y) is set up, as shown in Fig. 1, where x and y are normal to and along with the shock, respectively. u₁ and u₂ denote the velocities of the uniform flows at upstream and downstream of the shock, respectively. Their angles with respect to the shock are denoted by β₁ and β₂, respectively. We assume that the free-stream disturbance, which can be an acoustic, entropy, or vorticity wave, is of the form δAe^{i(k_xx+k_yy-ω t)}, where k_x and k_y are the components of wave vector k in (x, y), ω is the frequency, and δA represents the amplitude of any perturbed quantity. k₁, k₂, and k₃ represent the wave vectors of the incident wave, the acoustic wave leaving the shock, and the diverging entropy-vorticity wave, respectively. Their angles of wave vector with respect to the x-axis are θ₁, θ₂, and θ₃, respectively.

The total flow field can be written as the summation of steady base flow and a small perturbed quantity as follows:

(1)

Substitute (1) into the linearized Euler equations and solve the eigenvalue systems. We can get the dispersion relations and the corresponding eigenvectors for three elementary types of free-stream disturbances.

For acoustic waves, the dispersion relation is

(2)

and the corresponding eigenvector is

(3)

where c is the sound speed, , and θ is the angle between the wave vector and the x-axis. Plus (minus) sign corresponds to fast (slow) acoustic waves.

For entropy and vorticity waves, the dispersion relations are both

(4)

and the eigenvectors are, respectively,

(5)

(6)

It is worth noting that the direction of k indicates the phase direction of the disturbance, which is perpendicular to the wavefronts. It is different from the direction of group velocity v_g, i.e., ray. We depict the phase and the ray for the incident entropy/vorticity wave in Fig. 2. For the entropy/vorticity wave, v_g is the velocity of the moving fluid, i.e., u₁, while for the acoustic wave, the group velocity is given by

(7)

where plus (minus) corresponds to fast (slow) acoustic wave.

2.2 Wave vectors

Three conditions need to be satisfied across the shock:

(ⅰ) The frequency is continuous on both sides of the shock.

(ⅱ) The component of wave vector tangential to the shock, say, k_y is continuous.

(ⅲ) The component of group velocity normal to the shock should be larger than zero.

Condition (ⅲ) is deduced from the assumption that all the disturbances behind the shock should propagate downstream, which is true for the interaction taking place in the free space.

Since the waves are of small amplitude, these conditions can be applied to the unperturbed shock in the linear approximation. For acoustic waves generated by the entropy/vorticity wave, considering the dispersion relations given by (2) and (4), the application of Condition (ⅰ) gives

(8)

where , and subscripts 1 and 2 of the unperturbed quantities denote the quantities at upstream and downstream of the shock, respectively. From the first-order of continuity of unperturbed velocity tangential to the shock, we have v₁=v₂. From Condition (ⅱ), we have k_1y=k_2y=k_y, and then (8) reduces to

(9)

where u₂ and c₂ can be computed by the Rankine-Hugoniot (R-H) conditions across the unperturbed shock. k_2x can be solved from (9), and then the angle of the wave vector is given by tanθ₂=k_y/k_2x.

Similarly, for the entropy/vorticity wave behind the shock, the wave vector can be solved by

According to the existent condition of (9), we can deduce that the incident angle θ₁ should satisfy (details are referred to Appendix A)

(10)

where the subscript n represents components normal to the shock. Considering that the physical frequency should be larger than zero, by rewriting the dispersion relation of (4), the incident angle θ₁ for the entropy/vorticity wave should also satisfy

and thus

(11)

Combining (10) and (11), we depict the range of incident angle for entropy/vorticity wave to generate acoustic wave in Fig. 3. In Region Ⅰ, only fast acoustic wave can be generated behind the shock. In Region Ⅱ, only a slow acoustic wave can be generated. Outside Regions Ⅰ and Ⅱ, no acoustic wave can be generated.

For Ma₂ < 1,

For Ma₂>1,

2.3 Wave amplitude

Small perturbations are introduced in the R-H jump conditions on either side of the shock, and the linearized jump relations for small perturbations across the shock obtained are

(12a)

(12b)

(12c)

(12d)

where subscripts n and T denote the components normal and tangential to the shock, respectively. h is enthalpy.

Assume that the state equation is in the form of ρ=ρ(p, s), where s is entropy. Small fluctuations in density are related to pressure and entropy fluctuations by

(13)

where for a perfect gas, and c_p is the heat capacity at the constant pressure.

Combine the shock adiabatic, which is in the form of ρ₂=ρ₂(p₂, ρ₁, p₁), and the relation of small-perturbation of pressure, entropy, and density can be obtained to replace (12d),

(14)

where

The deformation of the shock can be assumed to be in the same form of the incident disturbance,

(15)

where η is the amplitude of the shock distortion and obviously η≪1. The velocity of the shock u_s is

(16)

Let , and thus . The velocity of the fluid relative to the shock is u+δ u -u_s. Substitute its normal and tangential components into (12a) - (12c), and the term can be eliminated. Combining with (13) and (14), a system of four linear equations with four unknowns, i.e., δd, δ p₂, δs₂, and δu₂, is obtained, which can be solved.

3 Numerical simulation

Numerical simulation is performed to get the same result for the interaction of the small disturbance with the shock. The governing equations are the two-dimensional unsteady compressible Euler equations in the conservation form,

(17)

where U is the flow quantities, and E and F are the nonlinear terms, including the pressure gradient. The Cartesian coordinate (x, y) is set as shown in Fig. 1. The unperturbed flow quantities before the shock, i.e., u₁, ρ₁, and ρ₁ u₁², are used to non-dimensionalize the velocity, density, and pressure, respectively. We choose a common and moderate accuracy shock-capturing scheme, i.e., a third-order weighted essentially nonoscillatory (WENO) scheme, to test its accuracy to simulate the fluctuations behind the shock. A third-order total-variation-diminishing (TVD) Runge-Kutta scheme is used for time advancement. As for the boundary conditions, the periodic condition is used for the upper and lower boundaries, and a sponge region is used as outflow boundary condition.

For our problem, since the wave number keeps unchanging in the y-direction, the computation can actually be reduced to a one-dimensional problem by performing the Fourier transform in y. However, in this paper, we still perform a two-dimensional computation with only considering one wave length in the y-direction so that the total computational work is well acceptable. Figure 4 shows the distribution of the unperturbed shock which locates at x=0.3. The Mach number of the oncoming flow is taken as 8, and the angle between the flow and the shock β₁ is 35°.

4 Results and discussion

The computation is performed first without imposing any disturbance. When the flow filed is completely steady, the small amplitude disturbance is introduced at the inlet. The flow quantities at the inlet are given by

(18)

(19)

where δA=10^-5, k_y=16π, θ₁=45°, and ω can be computed from the dispersion relation, i.e., ω=(u₁/tanθ₁+v₁) k_y=66.4. Note that there is no length scale in our problem. Therefore, the value of k_y is taken arbitrarily. The flow parameters and quantities are given in Table 1. Case 1 and Case 2 correspond to the cases where the flow behind the shock is subsonic and supersonic, respectively.

4.1 Case 1 where Ma₂ < 1

Figures 5(a) - 5(c) show the instantaneous perturbations along x at a specific y location. At this time, the generated fast acoustic wave has already passed through the outlet boundary of the computational domain, while the entropy and vorticity waves still remain in the computational domain. As shown in Figs. 5(a) and 5(b), the density and streamwise velocity fluctuations show a pattern of a wave packet in the range of x=0.30 -0.70, where the transmitted entropy wave, the generated vorticity wave, and the fast acoustic wave are mixed. The pressure perturbation shown in Fig. 5(c) has a sine-wave pattern because only the generated acoustic wave has the component of pressure fluctuation.

Since there is no nonlinear interaction between disturbances, we can separate the different types of disturbance from the flow perturbation field. Figures 6(a) and 6(b) show δρ for the entropy wave and δu for the vorticity wave, respectively.

From the numerical results above, we can get the wave length of generated/transmitted wave along the x-and y-axes, respectively, i.e., λ_x, and λ_y. Then, the wave angle θ can be computed from

(20)

The wave length λ of transmitted/generated wave is given by

(21)

Since we are especially interested in the acoustic wave generated behind the shock, to evaluate the efficiency of generation of the acoustic wave, a coefficient of generation is defined as, for the entropy incidence,

and for the vortical incidence,

where subscripts E and V denote the fluctuations of entropy and vorticity, respectively.

According to Section 2, the wave length of the disturbances behind the shock given by the theoretical analysis is

(22)

(23)

where sign plus (minus) corresponds to fast (slow) acoustic wave.

We compare the numerical results with the theoretical results in Table 2. δp_A denotes the pressure fluctuation of acoustic wave. The negative value of δA_V means that its phase has a lag of π compared with the incident wave. We can see very good agreement between the results from the two methods. Furthermore, note that the pressure fluctuations, indicated as well as in the generation coefficient C, which represent the acoustic wave of more importance in the receptivity problem, are in perfect agreement.

We also use the vorticity wave as the inlet disturbance to perform the same computation. The comparison of results given by two methods is shown in Table 3. Also, very good agreement can be found.

4.2 Case 2 where Ma₂>1

For the case that the flow is supersonic behind the shock, both slow and fast acoustic waves can be generated behind the shock. Two angles are chosen for the incident entropy wave, say, 45° and 135°. Table 4 shows the parameters used in two cases. Note that θ₁=45°, which lies in Regions I of Fig. 3(b), while θ₁=135° lies in Region Ⅱ. Therefore, Cases 2a and 2b correspond to the cases that fast and slow acoustic waves are generated behind the shock, respectively.

Figures 7(a) - 7(d) show the instantaneous contour of flow fluctuations for Cases 2a and 2b at a specific time when the wave with higher group velocity has passed through the computational domain, while the other waves still remain in the domain. We compare the results with those given by theoretical analysis in Table 5. As before, very good agreement can be found. Especially, for the acoustic waves, the shock-capturing method gives very accurate results.

It should be noticed that for Case 2b, the generated slow acoustic wave still propagates faster than the entropy/vorticity wave, whose group velocity is u₂, so that it goes out of the computational domain first. This can be explained by recalling the group velocity of the acoustic wave in (7). For Case 2b, since cosθ₂ < 0, the x-component of group velocity is v_gx=u₂-c₂ cosθ₂>u₂. Therefore, we see that the slow acoustic wave goes out of the computational domain first.

Taking the vorticity wave as the inlet disturbance, no essential difference can be found. Therefore, we omit the results here.

5 Summary and conclusions

In this paper, the interaction of small amplitude free-stream disturbances with an oblique shock is investigated. Considering the real situation in high altitude environment that flying vehicles may encounter, only the incidence of entropy/vorticity wave is considered. The numerical simulation is performed, using a common and moderate order accuracy shock-capturing method, say, the third-order WENO. The results are compared with those given by theoretical analysis. Satisfactory agreement is achieved, confirming that the shock-capturing method can faithfully reproduce the interaction of weak external disturbances with the shock.

Acknowledgements The first author would like to thank Professor Heng ZHOU in Tianjin University for his continuous encouragement and Professor Jisheng LUO in Tianjin University for his inspiring discussion and helpful suggestions.