- Email: [email protected]
boundary-layer interaction in a duct
Fluid Dynamics North-Holland
Research
217
5 (1989) 217-234
Computational study of the shock-wave/boundary-layer interaction in a duct Itaru Hataue Department Received
af Applied Physics, The University of Tokyo, Hongo, Bunkyo-ku,
20 September
Tokyo, 113, Japan
1988
Abstract. The compressible flow in the interior of a duct has been calculated by using the shock capturing method (Harten-type second-order TVD scheme) with LU factorization based on the two- and three-dimensional Navier-Stokes equations. The influences of some parameters contributing to the shock-wave/boundarylayer interaction are discussed. Some unsteady mechanisms of the interaction have been clarified and the computed results have made it possible to elucidate the details of the experimental observations.
1. In~~uction When a normal shock wave interacts strongly with a boundary layer, a short distance ahead of the point where the essentially perpendicular shock wave impinges on the boundary layer there appears a short leg forming a so-called h-shock. As the the interaction becomes stronger, the shape of the shock wave changes from A-type to X-type. Another interesting phenomenon caused by the shock-wave/boundary-layer interaction is the pseudo-shock wave. As a pseudoshock wave we consider a series of shock waves, followed by accelerated regions, that occurs when a strong normal shock wave interacts with a thick boundary layer. Back flow and separation of the boundary layer might occur along the walls. Finally, when the flow reaches subsonic speeds, it fills the duct. These phenomena are important not only as typical nonlinear problems but also as technicat problems in the design or operation of wind tunnel diffusers and in the design of high-speed centrifugal compressors. There have been many experimental investigations, using Schlieren photographs and pressure fluctuation measurements, of walls and interiors (Neumann and Lustwerk, 1949; Waltrup et al., 1973, 1974; Merkli, 1976; Om and Childs, 1985; Tamaki et al. 1969, 1970; Ikui et al. 1973; North and Stuart, 1962). Parameters contributing to the shock-wave/boundary-layer interaction, such as the upstream Mach number, the Reynolds number, the upstream boundary-layer thickness, and so on, have also been reported in detail in some papers. Theoretical studies have also been performed, for example that by use of the shockless model (Crocco, 1958). The rapid progress of the computer and the development of excellent schemes and algorithms for viscous gas dynamic equations have enabled us to obtain a great deal of detailed information (Viegas et al., 1978, 1979). In the present work, two-dimensional strong shock-wave/ boundary-layer interaction is been studied by solving the Navier-Stokes equation directly. Calculated examples are as follows: Case (1) normal shock-wave/bound~-layer interaction. Case (2) incident static shock-wave/bounda~-layer interaction. Case (3) reflected shock-wave/boundary-layer interaction. Case (4) propagating shock-wave/boundary-layer interaction. 0169-5983/89/$5.50
0 1989, The Japan
Society
of Fluid
Mechanics
218
I. Haiaue / Shock-wave/boundary-layer
interaction
Because the efficiency of computation in case (1) is not good, some parameter studies have been done in case (2). Furthermore, to discuss the three-dimensional effect, we studied two cases in which the cross-section of the duct is rectangular or circular and compared them with the two-dimensional case (result of case (2)). The purposes of this paper are to analyze the mechanism of unsteady bifurcation of some normal shock waves and the process of X-type shock wave change to X-type shock wave and/or pseudo-shock wave, and to discuss the influences of some parameters on the shock-wave/ boundary-layer interaction.
2. Method 2.1. The compressible
Navier-Stokes
Consider a two-dimensional conservation-law form under
equation
compressible Navier-Stokes the generalyzed independent
equation in nondimensional variable transformation:
~~+~~+~~=Re-‘(Rg+~~),
strong
(1)
where I
‘P\ fi=
PU
J-1
PU
i=
puu+
J-’
’
\
PV puv+
%P
pvV+
77yP
, (e+p)VI
s^Z
J-1
The variables
For viscous
kxy ,
sp
\
5 x%x + ~y~xy
hi=-’
+ E&
E,r+5,s
’ j
0 \ ~1*?xx+ 77yQy .
Il,r+17,s
U and V are called the contravariant
u=
0
I ,
71x%, + S,Tyy \
’
(e+p)U-&p,
I J-1
5,P
PVUf&P
e
k’=
1
PU
E, + 5,u + E,u,
v=
velocities
and defined
as
n)lt+ qxu + TJYV.
terms, we have 7,,=(X+2/+,+hvy, Y=U7xx+v7*Y+p(c*)x, p=(y-l)[e-&p(u2+v2)],
rx,=j_+y+v,),
7yy=(h+2~)vy+XU,,
s=U7,Y+v7YY+p(c~)Y,
p=
K Pr( y - 1) ’
c’=F,
with K and Pr denoting the thermal conductivity and the Prandtl number, respectively. Assuming an ideal gas of diatomic molecules, the specific heat ratio y is taken to be 1.4. The molecular viscosity is calculated by use of Sutherland’s formula. The turbulent eddy viscosity, pt, is also computed using a two-layer algebraic viscosity model developed by Baldwin and Lomax (1978) in some cases. J is the Jacobian and metrics are evaluated using second-order
I. Hataue / Shock-wave/boundary-layer
219
mm-action
central difference formulas for interior points and three-points one-sided formulas at boundary approximation has been used in eq. (1) points. In the three-dimensional case, the “thin-layer” only to the TI and the [ directions. The Navier-Stokes equation is written as follows:
(2) 2.2. Numerical
algorithm
(1) Explicit part. For convection terms, some TVD (total variation diminishing) schemes (Roe, 1981; Chakravarthy and Osher, 1985; Harten, 1984) have been proposed as shock-capturing techniques. In the present paper, a five-point TVD scheme proposed by Harten (1984) is adopted. Consider a two-dimensional system of hyperbolic conservation law: U, + E, + Fy = 0. Under
the generalized
(3)
coordinate
transformation
eq. (3) is rewritten
as follows:
U+@&~=o,
(4)
where l? = U/J, I? = (5, E + t,F)/J, and P = ( nXE + qF F)/J. Eq. (4) corresponds to the convection part in eq. (1). Using the Jacobian matrices A = ~E/S.J and B = aF/alJ, 2 and L? can be written as a=(&A+&B),
i=(q,A+q,B).
The Jacobian matrices of generalized coordinate system a and i can be diagonalized by the eigenvectors R, and R,. We can obtain eigenvalues a; (i = 1-4) of A and a; (i = l-4) of h. In view of the fact that any physical phenomenon propagates along the characteristics, the second-order TVD scheme can be written as
(5) The numerical flux functions gJ+ ,,Z,k and flux limiter @,+1,2 or @k+1,2 as ‘J+l,Z,k
=
@J.k
e.k+l,Z
=
t(4.k
+
+
‘J+l.k
e.k+l
+
+
4,k+ t,* can be expressed
R/+1,2.k@~+1,2.k)~
RJ,k+&jJ,ktl,z)>
by using
the concept
of
(6) (7)
where QJ+ 1/l k and QJ k+ 1,2 retain the TVD property. For the viscous terms, an accuracy up to second order in space is adopted. The MacCormack predictor-corrector, explicit scheme (MacCormack, 1969) retains an accuracy up to second order in time, and the implicit scheme (LU-factorized form described below) retains to first order. We computed some unsteady problems of the shock-wave/boundary-layer interaction by two different schemes, the MacCormack explicit scheme and the LU-AD1 implicit scheme, and compared the results. There is no large difference between them. Therefore, we discuss in the present paper the results calculated by use of the implicit scheme, because the implicit scheme is more efficient than the explicit scheme, especially in the three-dimensional calculation. (2) Implicit part. It is necessary to use an implicit scheme to improve the efficiency of the computation. The implicit scheme used in this paper is based on the Beam-Warming-Steger method (Beam and Warming, 1976; Steger, 1978) which is based on approximate factorization.
220
I. Hataue / Shock-wave/boundary-laker
This algorithm has first-order written as follows:
accuracy
in time and
i~teruct~o~
second-order
accuracy
in space
and is
where a = TsDAT;-‘, g = TvDBT;’ (both DA and D, are diagonal matrices (Pulliam et al., 1981, 1985). The central differencing can be decomposed into one-sided differencings by using the flux vector splitting technique (Obayashi and Kuwahara, 1984), i.e. (I+hs,A^)=T#+vcD;
+A,D,-)r,-‘,
(9a)
(I+~~~~)=T~(~+v~D~
+b,D,)Tq-‘,
t9b)
with D * = ih(D t_ 1D I), where h = At, 6 is the central finite-difference operator, and A and v are forward and backward difference operators, respectively. The right-hand side of eq. (9a) can be rewritten as T&A
+ MA f N,)T;‘,
where for three-point
upwinding
(10) L,,
L, = - %Q,;_, -I- +I&, The diagonally dominant here can be described as
factorization
MA and NA can be written MA=I+$(D;-D<), (Lombard
as follows: NA=:D;+,-;D&.
et al., 1983; Fujii and Obayashi,
1986) used
L, + MA + NA = (LA + MA) M,-‘( MA + N,) f O(h2). Here, the block pentadiagonal operator is factorized as the product of the lower and upper scalar tridiagonal ones. In order to treat viscous terms, the split Jacobian matrices a* are modified as
2.3. Grid systems The grid systems used in this paper are shown in fig. 1. The numbers of the grid points are (a) 350 X 77, (b) 350 X 153 in the case of two-dimensional calculations, and (c) 75 X 49 X 49 (rectangular cross section), (d) 175 x 60 X 39 (circular cross section) in three-dimensional computations, For the number of grid points (1 x J x K), I is the number of grid points of the 6 direction, J is that of the q direction and K is that of the { direction. Near the wall, the grid points are concentrated and the nondimensionalized least distance between adjacent grid points is set to be less than 0.1/a. The velocity is normalized by the velocity of sound, and the length is normalized by the width of the duct. 2.4. Initial conditions For case (l), the initial velocity is taken to be zero everywhere except at the inlet. For cases (2) and (4) a uniform flow and an incident shock wave are given. The slow-start technique is used during the first 200 time steps for the impulsive start in cases (1) (2), and (4). In case (3), there is a diaphragm that separates the duct in two rooms. There is a perfect gas in each room
I. Hataue / Shock-wave/boundary-layer
interaction
the
221
centrat‘axis((=l)
(4 Fig. 1. Grid systems in a duct for two-dimensional ((c) 75 x 49 x 49, (d) 175 x 60 x 39).
calculations
((a) 350 X 77, (b) 350 X 153) and three-dimensional
ones
with the same uniform temperature but different densities. The initial density ratio is taken to be 1: 30. The flow is started by the rupture of the diaphragm. The same initial conditions as in case (2) are assumed in the three-Dimensions calculations.
2.5.
Boundary
conditbons
The boundary conditions are as follows: on the wall of the duct, the no-slip condition is required for viscous flow and the density p is given by the zero gradient conditions (+/a?7 = 0). The pressure p is computed by solving the normal momentum equation. At the inlet (left side boundary), a fixed uniform flow condition is imposed in case (1) and the boundary conditions in cases (2)-(4) are given by the zero gradient conditions (at/at = 0). At the outlet (right side boundary) the boundary conditions are given in the same manner at the inlet in cases (2) and (4). As for the boundary condition of the right side boundary in case (3) where an incident shock-wave reftects, velocity is zero and other physical quantities are the same as those as inner
222
I. Hataue / Shock-wave/boundagAayer
w&r-action
adjacent points. The same boundary conditions are imposed also in the three-dimensional calculation in which the cross section of the duct in rectangular. The boundary conditions on the wall and at the inlet and outlet of the pipe with circular cross section are the same as the two-dimensional calculation. All values on the central axis (l = 1) are averaged from the extrapolated ones at adjacent points. The boundary condition of the r~ direction is given as follows. Two grid points are overlapped, for example p ,,,, k ,,,,,.+, = p ,,,, 1, p, ., ,k ml.,\= p ,,,, O, and so on. All the calculations were done on the Hitachi S820 of the Computer Center, the University of Tokyo.
Results 3.1. Two-dimensional
case
Case (1): The calculation for Mach number M = 2 and the Reynolds number Re = 10’ without turbulence model. The fixed uniform flow condition at the inlet is chosen in this case (fig. 2a). Fig. 2 shows density contours (figs. 2b-2f) and pressure contours (fig. 2g). Time proceeds from fig. 2b to 2f. Normalized time is from T = 2 to T = 12. A strong normal shock wave is formed near the inlet. Then, by the shock-wave/boundary-layer interaction, the shape of the shock wave changes from normal to X-type (fig. 2b). As the interaction becomes stronger, the shape changes to X-type. The boundary layer grows thick behind the front shock wave, while the rear one is not very clear (fig. 2~). Finally, after the first X-type shock wave, the second normal shock-wave is formed from both side walls to the center region (fig. 2d). The second shock wave is also of normal type at the beginning, but the shape eventually changes to h-type and X-type (fig. 2e). Furthermore, behind the second X-type shock wave, a third one is formed and grows (fig. 2f). The number of shock waves forming a pseudo-shock wave is reported to be about ten in the case of the strong shock-wave/ boundary-layer interaction on the basis of experimental data. Our calculation in this case ends here. However, if the calculation were continued further, we would have more then three shock waves. Fig. 3 shows the plots of velocity vectors superimposed on the density contours. The number of the density contours in fig. 3 is less than that in fig. 2. The front oblique shock wave changes the flow direction towards the center and several vortices appear near the separated boundary layer. When the rear oblique shock wave makes the flow direction change towards the side walls again, the flow is accelerated and then a second normal shock-wave is formed. Fig. 4 shows Schlieren photographs of the pseudo-shock wave by Tamaki et al. (1969). Though the calculation has not reached steady state yet, the calculated result agrees well with the experimental one, generally. Fig. 5 shows a comparison of the computed pressure fluctuation near the wall and in the interior with experimental data of Tamaki et al. (1970). The pressure values in all of the calculated region are normalized between maximum and minimum values. The global profile is that pressure increases gradually from first shock waves to next ones. The oscillation of the pressure becomes stronger in the interior, especially around the center line, that near the wall. This result is in good qualitative agreement with experimental data. The pressure oscillates a little near the wall. This is supposed to be due to the influence of strong shock-wave/ boundary-layer interaction. Case (2): The calculation for (a) Mach number M = 2 and M = 3, (b) Reynolds number Re = lo5 and Re = 107) and (c) with or without the turbulence model. In this case, incidentnormal-shock-wave/boundary-layer interaction is calculated (fig. 6a). We discuss the relationship between some parameters (the strength of the shock wave and the Reynolds number) and the shock-wave/boundary-layer interaction. The Reynolds number (Re = 105) and the upstream Mach number (M = 2) are the same as in case (1) but the upstream density and
I. Hataue / Shock-wave/boundary-layer
interaction
223
(a)
(b) T=2
(c) T=3
(d) T=4
(e) T=5
(f) T=12
(g) T=12
Fig. 2. Density nondimensional
contours ones.
(b-f)
and pressure
contours
(g) at M=
2. Re =lO’,
without
turbulence
model.
Times
are
pressure are different from those in case (1). Therefore, the strength of the shock wave is different in cases (1) and (2). Fig. 6b shows the density contours for M = 2, Re = lo5 (number of grid points = 350 X 77) without turbulence model. As time proceeds, the shape of the shock-wave changes from normal-type to X-type and X-type, and eventually a pseudo-shock wave consisting of many X-type shock waves appears. Fig. 6c shows the wall pressure fluctuation. A gradual pressure increase also appears in this case. Fig. 7 shows the density contours for M = 3, Re = 10’ (number of the grid points = 350 X 77) without turbulence model. In the case of M = 2 (fig. 6) each X-type shock wave is very clear, and the structure of the pseudo-shock wave is periodical. On the contrary, the structure is very complicated in the case of M = 3, so that the third shock wave already is not clear. Fig. 8 shows the experimental comparison of the structure of the pseudo-shock wave with different Mach numbers using Schlieren photographs (by Ikui et al. 1973). Those pictures show that a clear pseudo-shock wave cannot be seen when the upstream
I. Hataue / Shock-wave/boundary-layer
224
density
/
interaction
(a) T=5
\ velocity vectors
COD~OU~~.
/
[b) T=lZ
Fig. 3. Plots of velocity vectors superimposed
on the density
contours
at A4 = 2, Re = 105, without
turbulence
model.
Mach number is high. This agrees with the calculated results (figs. 6b and 7). Fig. 9 shows the calculated dependence on the number of grid points. Though it shows that sharper shock waves and finer structures of the shock-wave/boundary-layer interaction are captured in the finer grid system, there is little difference in the general structure. Therefore, in discussing the global structure of the shock-wave/ boundary-layer interaction, results calculated by use of a not very fine grid are regarded as sufficient in order to capture the shock-wave/boundary-layer interaction. In fig. 10 are compared the densities in the case of M = 3, Re = lo5 with those for M = 3, Re = 10’. There is only a minor difference owing to a difference in the grid systems. The interaction is still regarded as strong in the case of Re = 107. In fig. 11 the cases calculated with and without the turbulence model for A4 = 3, Re = 10’ are compared. There is no large difference between the two figures. Case (3): The calculation of the interaction of the reflected shock-wave/boundary-layer interaction, for Re = lo5 without the turbulence model. In this case, we discuss the relationship between the thickness fo the upstream boundary layer and the shock-wave/ boundary-layer interaction. The Mach number of the reflected shock wave relative to the upstream one is 2.2. Fig. 12 shows the density contours. The formed normal shock wave propagates (fig. 12a), reflects and propagates towards the opposite direction, then the interaction occurs (fig. 12b) and the shape of the reflected shock wave changes to X-type and X-type (fig. 12c-12d). Fig. 13 shows Schlieren photographs of the reflected shock-wave/boundary-layer interaction taken by North and Stuart (1962). Our calculated results are in very good agreement with their
Fig. 4. Schlieren
photograph
by Tamaki
et al
I. Hataue / Shock-wave/boundary-layer
interaction
225
h
Pmin Fig. 5. Comparison of computed by Tamaki et al. (b). (- - -)
wall and interior pressure fluctuation at T = 12 (a) with experimental determined y/h =l.O, (- - -- --) y/h = 0.5, (.-.-.) y/h = 0.25, ( -) y/h=O.
ones
experimental data. Fig. 13b corresponds to fig. 12~. In this case, the upstream boundary layer becomes thicker through the influence of the propagating incident shock wave. As a consequence, the amplitude of the reflected shock-wave/boundary-layer interaction becomes larger and the X-type shock wave is formed faster. Case (4): The calculation of the interaction of propagating shock-wave/ boundary-layer interaction for Re = 105. Comparing with the results of cases (2) and (3), we discuss here the relationship between the parameters (the strength of the shock wave and the thickness of the upstream boundary layer) and the shock-wave/boundary-layer interaction again. Fig. 14 shows the density contours. The Mach number of the shock wave (MS) in this case is about 2.6. The shock wave proceeds, firstly the bifurcation occurs, then the shape changes to X-type and the second shock wave (pseudo-shock wave) is formed. We computed the flows for several Mach numbers of the shock wave. When the Mach number of the shock wave is less than 2.4, bifurcation of the normal shock wave never occurs and the bifurcation begins very late for MS = 2.6. That is the reason why the boundary layer is very thin, and therefore the interaction is weak. 3.2. Three-dimensional
case
Two kinds of three-dimensional calculation are performed under the same condition case (2) of the two-dimensional case (M = 3, Re = 105, without the turbulence model).
as the
I. Hataue / Shock-wave/boundary-layer
226
incident
interaction shock -wave
~
Fig. 6. Density
contours
at M = 2, Re = 10’. without
Fig. 7. Density
contours
turbulence
at M = 3, Re = lo’, without
model and pressure
turbulence
fluctuation.
model.
(b) Fig. 8. Schlieren
photograph
by Ikui et al. at Re = IO’. (a) M = 2.15, (b) A4 = 2.68.
I. Hatme / Shock-waae/boundary-la~ver interaction
227
(4
(b) Fig. 9. Density contours at M = 3, Re = 105, without turbulence model. (a) 350 X 77 grid points (same calculations as fig. 7 but at different nondimensional time), (b) 350x 153 grid points. (b) 350 x 153 grid points.
Case (5): Duct with rectangular cross section. Fig. 15 shows the density contours in two different horizontal planes. Comparing the figure of the central plane (fig. 15b) with that of the two-Dimensions case (fig. 15c), we find that there is no large difference. In fig. 16 shows the velocity vectors are plotted, projected on each cross section superimposed on the density contours. The patterns of the secondary flow induced by the shock-wave/ boundary-layer interaction are different for the upstream side and the downstream side of the front oblique X-type shock wave (figs.16c, 16e). In the upstream region of the front oblique shock wave, the
(4
fbf Fig. 10. Density contours at (a) M = 3, Re = lo5 (b) M = 3, Re = lo’, without turbulence model.
228
1. Hatme / Shack-waoe/boundat-iuy~r
interartiort
(b) Fig. 11. Density
contours
(a) without
turbulence
model (b) with turbulence
model at M = 3, Re = 105.
(a) T-5.0
(b) Tz5.7
(c) T=6.4
(e) T=7.9
Fig. 12. Density
contours.
(a) Incident
shock wave, (b)-(e)
reflected
shock wave,
I. Hataue / Shock-waoe/boundary-layer
229
interaction
(4 Fig. 13. Schlieren
photograph
of reflected
shock-wave/boundary-layer
interaction.
secondary flow occurs from the corner towards the center along the side wall in the separated boundary layer. When the direction of the main flow is toward the wall behind the oblique shock wave, the flow in the cross section is in the radial direction. Fig. 17 shows the
I--
(a) T=4.8
incident
propagating
shock-wave
-
(b) T=5.5
Fig. 14. Density
contours
in the case of propagating
shock wave.
I. Hatme
230
/ Shock-wave/boundary-lacer
interactron
(b)
Fig. 15. Density plane.
contours
in three-dimensional
case (a, b) and the two-dimensional
case (c). (b) is the one of the central
instantaneous streamlines projected to each cross section. Many vortices appear in the separated boundary layer. Case (6): Pipe with circular cross section. Fig. 18 gives density contours in a plane of symmetry and velocity vectors projected on each cross section superimposed on the density contours. Fig. 19 shows the instantaneous streamlines in each cross section. The density contours are almost axially symmetrical. Between the density contours in the symmetrical plane and that in the central horizontal plane in the case of a rectangular duct (fig. 15b), we find no large difference. However, the streamlines in the cross section are almost axially symmetrical.
I. Hataue / Shock-wave/boundary-layer
interaction
231
fdl (40/75) Fig. 16. Plots of velocity vectors in each cross section mth grid plane numbered from the inlet of the total
superimposed on the density n grid planes.
(25175)
Fig. 17. Streamlines
in each cross section.
contours.
(e) (W75) Symbol
(m/n)
means the
I. Hataue / Shock-wave/boundary-layer
232
(20/175)
interaction
(60/175)
(80/175)
(140/175)
(120/175) Fig. 18. Density contours in the plane superimposed on the density contours.
of symmetry
(top figure),
(!60/175) and plots of velocity
vectors
in each cross section
4. Conclusion The shock-wave/boundary-layer interaction was successfully computed on the basis of twoand three-dimensional Navier-Stokes equations by the shock capturing method (Harten-type second-order TVD scheme) with LU factorization. The mechanisms of bifurcation and formation of X-type shock wave, X-type shock wave, and pseudo-shock wave were captured. By calculating several cases, the dependence of the flow pattern on the upstream Mach number and the upstream boundary-layer thickness was clarified. Furthermore, the secondary flow showing the three-dimensionally of the flow is clearly observed.
I. Hataue / Shock-waue/boundary-layer
(110/175)
(120/175) Fig. 19. Streamlines
interaction
(140/175)
233
(160/175)
in each cross section.
Acknowledgement
The author is grateful to Professor Hideo Takami of the University of Tokyo for his valuable discussion and to Professor Kunio Kuwahara of ISAS (Institute of Space and Astronautical Science) and Dr. Katsuya Ishii of ICFD (Institute of Computational Fluid Dynamics) for their kind suggestions.
References
Baldwin, B.S. and H. Lomax (1978) AIAA paper 78-0257. Beam, R. and R.F. Warming (1976) J. Comput. Phys. 22, 87-110. Chakravarthy, S.R. and S. Osher (1985) AIAA paper 85-0363. Crocco, L. (19581 Fundamentals of Gas Dy~amics~ Section B {P~nceton) p. 124. Fujii, K. and S. Obayashi (1986) AIAA paper 86-1831. Harten, A. (1984) SIAM J. Numer. Anal. 21, l-23. Ikui, T., K. Matsuo and M. Nagai (1973) Trans. Japan Sot. Mech. Eng. B 39-326. 3054-3064. Lombard, C.K., J. Bardina, E. Venkatapathy and J. Oliger (1983) AIAA paper 83 - 189SCP. MacCormack, R.W. (1969) AIAA paper 69-354. Merkli, P.E. (1976) AIAA J. 14, 168-172. Neumann, E.P. and F. Lustwerk (1949) J. Appl. Mech. 16, 195-202. North, R.J. and C.M. Stuart (1962) Nat. Phys. Lab. Aeron. Rept. 1029. Obayashi, S. and K. Kuwahara (1984) AlAA paper 48.1670. Om, D. and M.E. Cbilds (1985) AIAA J. 23, 1506-1511. Pulliam, T.H. and D.S. Chaussee (1981) J. Comput. Phys. 39, 347-363. Pulliam, T.H. and J.L. Steger (1985) AIAA paper 85-0360. Roe, P.L. (1981) J. Compuf. Phys. 43, 357-372. Steger, J.L. (1978) AIAA J. 17, 249-255. Tamaki, S., Y. Torn&a and R. Yamane (1969) Trans. Japan Sot. Mech. Eng. B 35-273, 1028-1035
234
I. Hataue / Shock-wave/boundary-layer
interaction
Tamaki, S., Y. Tomita and R. Yamane (1970) Trans. Japan Sot. Mech. Eng. B 36-292, 2056-2066 Viegas, J.R. and T.J. Coakley (1978) AIAA J. 16, 293-294. Viegas, J.R. and CC. Horstman (1979) AIAA J. 17, 811-820. Waltrup, P.J. and F.S. Billig (1973) AIAA J. Il, 1404-1408. Waltrup, P.J. and J.M. Cameron (1974) AIAA J. 12, 878-880.
Research
217
5 (1989) 217-234
Computational study of the shock-wave/boundary-layer interaction in a duct Itaru Hataue Department Received
af Applied Physics, The University of Tokyo, Hongo, Bunkyo-ku,
20 September
Tokyo, 113, Japan
1988
Abstract. The compressible flow in the interior of a duct has been calculated by using the shock capturing method (Harten-type second-order TVD scheme) with LU factorization based on the two- and three-dimensional Navier-Stokes equations. The influences of some parameters contributing to the shock-wave/boundarylayer interaction are discussed. Some unsteady mechanisms of the interaction have been clarified and the computed results have made it possible to elucidate the details of the experimental observations.
1. In~~uction When a normal shock wave interacts strongly with a boundary layer, a short distance ahead of the point where the essentially perpendicular shock wave impinges on the boundary layer there appears a short leg forming a so-called h-shock. As the the interaction becomes stronger, the shape of the shock wave changes from A-type to X-type. Another interesting phenomenon caused by the shock-wave/boundary-layer interaction is the pseudo-shock wave. As a pseudoshock wave we consider a series of shock waves, followed by accelerated regions, that occurs when a strong normal shock wave interacts with a thick boundary layer. Back flow and separation of the boundary layer might occur along the walls. Finally, when the flow reaches subsonic speeds, it fills the duct. These phenomena are important not only as typical nonlinear problems but also as technicat problems in the design or operation of wind tunnel diffusers and in the design of high-speed centrifugal compressors. There have been many experimental investigations, using Schlieren photographs and pressure fluctuation measurements, of walls and interiors (Neumann and Lustwerk, 1949; Waltrup et al., 1973, 1974; Merkli, 1976; Om and Childs, 1985; Tamaki et al. 1969, 1970; Ikui et al. 1973; North and Stuart, 1962). Parameters contributing to the shock-wave/boundary-layer interaction, such as the upstream Mach number, the Reynolds number, the upstream boundary-layer thickness, and so on, have also been reported in detail in some papers. Theoretical studies have also been performed, for example that by use of the shockless model (Crocco, 1958). The rapid progress of the computer and the development of excellent schemes and algorithms for viscous gas dynamic equations have enabled us to obtain a great deal of detailed information (Viegas et al., 1978, 1979). In the present work, two-dimensional strong shock-wave/ boundary-layer interaction is been studied by solving the Navier-Stokes equation directly. Calculated examples are as follows: Case (1) normal shock-wave/bound~-layer interaction. Case (2) incident static shock-wave/bounda~-layer interaction. Case (3) reflected shock-wave/boundary-layer interaction. Case (4) propagating shock-wave/boundary-layer interaction. 0169-5983/89/$5.50
0 1989, The Japan
Society
of Fluid
Mechanics
218
I. Haiaue / Shock-wave/boundary-layer
interaction
Because the efficiency of computation in case (1) is not good, some parameter studies have been done in case (2). Furthermore, to discuss the three-dimensional effect, we studied two cases in which the cross-section of the duct is rectangular or circular and compared them with the two-dimensional case (result of case (2)). The purposes of this paper are to analyze the mechanism of unsteady bifurcation of some normal shock waves and the process of X-type shock wave change to X-type shock wave and/or pseudo-shock wave, and to discuss the influences of some parameters on the shock-wave/ boundary-layer interaction.
2. Method 2.1. The compressible
Navier-Stokes
Consider a two-dimensional conservation-law form under
equation
compressible Navier-Stokes the generalyzed independent
equation in nondimensional variable transformation:
~~+~~+~~=Re-‘(Rg+~~),
strong
(1)
where I
‘P\ fi=
PU
J-1
PU
i=
puu+
J-’
’
\
PV puv+
%P
pvV+
77yP
, (e+p)VI
s^Z
J-1
The variables
For viscous
kxy ,
sp
\
5 x%x + ~y~xy
hi=-’
+ E&
E,r+5,s
’ j
0 \ ~1*?xx+ 77yQy .
Il,r+17,s
U and V are called the contravariant
u=
0
I ,
71x%, + S,Tyy \
’
(e+p)U-&p,
I J-1
5,P
PVUf&P
e
k’=
1
PU
E, + 5,u + E,u,
v=
velocities
and defined
as
n)lt+ qxu + TJYV.
terms, we have 7,,=(X+2/+,+hvy, Y=U7xx+v7*Y+p(c*)x, p=(y-l)[e-&p(u2+v2)],
rx,=j_+y+v,),
7yy=(h+2~)vy+XU,,
s=U7,Y+v7YY+p(c~)Y,
p=
K Pr( y - 1) ’
c’=F,
with K and Pr denoting the thermal conductivity and the Prandtl number, respectively. Assuming an ideal gas of diatomic molecules, the specific heat ratio y is taken to be 1.4. The molecular viscosity is calculated by use of Sutherland’s formula. The turbulent eddy viscosity, pt, is also computed using a two-layer algebraic viscosity model developed by Baldwin and Lomax (1978) in some cases. J is the Jacobian and metrics are evaluated using second-order
I. Hataue / Shock-wave/boundary-layer
219
mm-action
central difference formulas for interior points and three-points one-sided formulas at boundary approximation has been used in eq. (1) points. In the three-dimensional case, the “thin-layer” only to the TI and the [ directions. The Navier-Stokes equation is written as follows:
(2) 2.2. Numerical
algorithm
(1) Explicit part. For convection terms, some TVD (total variation diminishing) schemes (Roe, 1981; Chakravarthy and Osher, 1985; Harten, 1984) have been proposed as shock-capturing techniques. In the present paper, a five-point TVD scheme proposed by Harten (1984) is adopted. Consider a two-dimensional system of hyperbolic conservation law: U, + E, + Fy = 0. Under
the generalized
(3)
coordinate
transformation
eq. (3) is rewritten
as follows:
U+@&~=o,
(4)
where l? = U/J, I? = (5, E + t,F)/J, and P = ( nXE + qF F)/J. Eq. (4) corresponds to the convection part in eq. (1). Using the Jacobian matrices A = ~E/S.J and B = aF/alJ, 2 and L? can be written as a=(&A+&B),
i=(q,A+q,B).
The Jacobian matrices of generalized coordinate system a and i can be diagonalized by the eigenvectors R, and R,. We can obtain eigenvalues a; (i = 1-4) of A and a; (i = l-4) of h. In view of the fact that any physical phenomenon propagates along the characteristics, the second-order TVD scheme can be written as
(5) The numerical flux functions gJ+ ,,Z,k and flux limiter @,+1,2 or @k+1,2 as ‘J+l,Z,k
=
@J.k
e.k+l,Z
=
t(4.k
+
+
‘J+l.k
e.k+l
+
+
4,k+ t,* can be expressed
R/+1,2.k@~+1,2.k)~
RJ,k+&jJ,ktl,z)>
by using
the concept
of
(6) (7)
where QJ+ 1/l k and QJ k+ 1,2 retain the TVD property. For the viscous terms, an accuracy up to second order in space is adopted. The MacCormack predictor-corrector, explicit scheme (MacCormack, 1969) retains an accuracy up to second order in time, and the implicit scheme (LU-factorized form described below) retains to first order. We computed some unsteady problems of the shock-wave/boundary-layer interaction by two different schemes, the MacCormack explicit scheme and the LU-AD1 implicit scheme, and compared the results. There is no large difference between them. Therefore, we discuss in the present paper the results calculated by use of the implicit scheme, because the implicit scheme is more efficient than the explicit scheme, especially in the three-dimensional calculation. (2) Implicit part. It is necessary to use an implicit scheme to improve the efficiency of the computation. The implicit scheme used in this paper is based on the Beam-Warming-Steger method (Beam and Warming, 1976; Steger, 1978) which is based on approximate factorization.
220
I. Hataue / Shock-wave/boundary-laker
This algorithm has first-order written as follows:
accuracy
in time and
i~teruct~o~
second-order
accuracy
in space
and is
where a = TsDAT;-‘, g = TvDBT;’ (both DA and D, are diagonal matrices (Pulliam et al., 1981, 1985). The central differencing can be decomposed into one-sided differencings by using the flux vector splitting technique (Obayashi and Kuwahara, 1984), i.e. (I+hs,A^)=T#+vcD;
+A,D,-)r,-‘,
(9a)
(I+~~~~)=T~(~+v~D~
+b,D,)Tq-‘,
t9b)
with D * = ih(D t_ 1D I), where h = At, 6 is the central finite-difference operator, and A and v are forward and backward difference operators, respectively. The right-hand side of eq. (9a) can be rewritten as T&A
+ MA f N,)T;‘,
where for three-point
upwinding
(10) L,,
L, = - %Q,;_, -I- +I&, The diagonally dominant here can be described as
factorization
MA and NA can be written MA=I+$(D;-D<), (Lombard
as follows: NA=:D;+,-;D&.
et al., 1983; Fujii and Obayashi,
1986) used
L, + MA + NA = (LA + MA) M,-‘( MA + N,) f O(h2). Here, the block pentadiagonal operator is factorized as the product of the lower and upper scalar tridiagonal ones. In order to treat viscous terms, the split Jacobian matrices a* are modified as
2.3. Grid systems The grid systems used in this paper are shown in fig. 1. The numbers of the grid points are (a) 350 X 77, (b) 350 X 153 in the case of two-dimensional calculations, and (c) 75 X 49 X 49 (rectangular cross section), (d) 175 x 60 X 39 (circular cross section) in three-dimensional computations, For the number of grid points (1 x J x K), I is the number of grid points of the 6 direction, J is that of the q direction and K is that of the { direction. Near the wall, the grid points are concentrated and the nondimensionalized least distance between adjacent grid points is set to be less than 0.1/a. The velocity is normalized by the velocity of sound, and the length is normalized by the width of the duct. 2.4. Initial conditions For case (l), the initial velocity is taken to be zero everywhere except at the inlet. For cases (2) and (4) a uniform flow and an incident shock wave are given. The slow-start technique is used during the first 200 time steps for the impulsive start in cases (1) (2), and (4). In case (3), there is a diaphragm that separates the duct in two rooms. There is a perfect gas in each room
I. Hataue / Shock-wave/boundary-layer
interaction
the
221
centrat‘axis((=l)
(4 Fig. 1. Grid systems in a duct for two-dimensional ((c) 75 x 49 x 49, (d) 175 x 60 x 39).
calculations
((a) 350 X 77, (b) 350 X 153) and three-dimensional
ones
with the same uniform temperature but different densities. The initial density ratio is taken to be 1: 30. The flow is started by the rupture of the diaphragm. The same initial conditions as in case (2) are assumed in the three-Dimensions calculations.
2.5.
Boundary
conditbons
The boundary conditions are as follows: on the wall of the duct, the no-slip condition is required for viscous flow and the density p is given by the zero gradient conditions (+/a?7 = 0). The pressure p is computed by solving the normal momentum equation. At the inlet (left side boundary), a fixed uniform flow condition is imposed in case (1) and the boundary conditions in cases (2)-(4) are given by the zero gradient conditions (at/at = 0). At the outlet (right side boundary) the boundary conditions are given in the same manner at the inlet in cases (2) and (4). As for the boundary condition of the right side boundary in case (3) where an incident shock-wave reftects, velocity is zero and other physical quantities are the same as those as inner
222
I. Hataue / Shock-wave/boundagAayer
w&r-action
adjacent points. The same boundary conditions are imposed also in the three-dimensional calculation in which the cross section of the duct in rectangular. The boundary conditions on the wall and at the inlet and outlet of the pipe with circular cross section are the same as the two-dimensional calculation. All values on the central axis (l = 1) are averaged from the extrapolated ones at adjacent points. The boundary condition of the r~ direction is given as follows. Two grid points are overlapped, for example p ,,,, k ,,,,,.+, = p ,,,, 1, p, ., ,k ml.,\= p ,,,, O, and so on. All the calculations were done on the Hitachi S820 of the Computer Center, the University of Tokyo.
Results 3.1. Two-dimensional
case
Case (1): The calculation for Mach number M = 2 and the Reynolds number Re = 10’ without turbulence model. The fixed uniform flow condition at the inlet is chosen in this case (fig. 2a). Fig. 2 shows density contours (figs. 2b-2f) and pressure contours (fig. 2g). Time proceeds from fig. 2b to 2f. Normalized time is from T = 2 to T = 12. A strong normal shock wave is formed near the inlet. Then, by the shock-wave/boundary-layer interaction, the shape of the shock wave changes from normal to X-type (fig. 2b). As the interaction becomes stronger, the shape changes to X-type. The boundary layer grows thick behind the front shock wave, while the rear one is not very clear (fig. 2~). Finally, after the first X-type shock wave, the second normal shock-wave is formed from both side walls to the center region (fig. 2d). The second shock wave is also of normal type at the beginning, but the shape eventually changes to h-type and X-type (fig. 2e). Furthermore, behind the second X-type shock wave, a third one is formed and grows (fig. 2f). The number of shock waves forming a pseudo-shock wave is reported to be about ten in the case of the strong shock-wave/ boundary-layer interaction on the basis of experimental data. Our calculation in this case ends here. However, if the calculation were continued further, we would have more then three shock waves. Fig. 3 shows the plots of velocity vectors superimposed on the density contours. The number of the density contours in fig. 3 is less than that in fig. 2. The front oblique shock wave changes the flow direction towards the center and several vortices appear near the separated boundary layer. When the rear oblique shock wave makes the flow direction change towards the side walls again, the flow is accelerated and then a second normal shock-wave is formed. Fig. 4 shows Schlieren photographs of the pseudo-shock wave by Tamaki et al. (1969). Though the calculation has not reached steady state yet, the calculated result agrees well with the experimental one, generally. Fig. 5 shows a comparison of the computed pressure fluctuation near the wall and in the interior with experimental data of Tamaki et al. (1970). The pressure values in all of the calculated region are normalized between maximum and minimum values. The global profile is that pressure increases gradually from first shock waves to next ones. The oscillation of the pressure becomes stronger in the interior, especially around the center line, that near the wall. This result is in good qualitative agreement with experimental data. The pressure oscillates a little near the wall. This is supposed to be due to the influence of strong shock-wave/ boundary-layer interaction. Case (2): The calculation for (a) Mach number M = 2 and M = 3, (b) Reynolds number Re = lo5 and Re = 107) and (c) with or without the turbulence model. In this case, incidentnormal-shock-wave/boundary-layer interaction is calculated (fig. 6a). We discuss the relationship between some parameters (the strength of the shock wave and the Reynolds number) and the shock-wave/boundary-layer interaction. The Reynolds number (Re = 105) and the upstream Mach number (M = 2) are the same as in case (1) but the upstream density and
I. Hataue / Shock-wave/boundary-layer
interaction
223
(a)
(b) T=2
(c) T=3
(d) T=4
(e) T=5
(f) T=12
(g) T=12
Fig. 2. Density nondimensional
contours ones.
(b-f)
and pressure
contours
(g) at M=
2. Re =lO’,
without
turbulence
model.
Times
are
pressure are different from those in case (1). Therefore, the strength of the shock wave is different in cases (1) and (2). Fig. 6b shows the density contours for M = 2, Re = lo5 (number of grid points = 350 X 77) without turbulence model. As time proceeds, the shape of the shock-wave changes from normal-type to X-type and X-type, and eventually a pseudo-shock wave consisting of many X-type shock waves appears. Fig. 6c shows the wall pressure fluctuation. A gradual pressure increase also appears in this case. Fig. 7 shows the density contours for M = 3, Re = 10’ (number of the grid points = 350 X 77) without turbulence model. In the case of M = 2 (fig. 6) each X-type shock wave is very clear, and the structure of the pseudo-shock wave is periodical. On the contrary, the structure is very complicated in the case of M = 3, so that the third shock wave already is not clear. Fig. 8 shows the experimental comparison of the structure of the pseudo-shock wave with different Mach numbers using Schlieren photographs (by Ikui et al. 1973). Those pictures show that a clear pseudo-shock wave cannot be seen when the upstream
I. Hataue / Shock-wave/boundary-layer
224
density
/
interaction
(a) T=5
\ velocity vectors
COD~OU~~.
/
[b) T=lZ
Fig. 3. Plots of velocity vectors superimposed
on the density
contours
at A4 = 2, Re = 105, without
turbulence
model.
Mach number is high. This agrees with the calculated results (figs. 6b and 7). Fig. 9 shows the calculated dependence on the number of grid points. Though it shows that sharper shock waves and finer structures of the shock-wave/boundary-layer interaction are captured in the finer grid system, there is little difference in the general structure. Therefore, in discussing the global structure of the shock-wave/ boundary-layer interaction, results calculated by use of a not very fine grid are regarded as sufficient in order to capture the shock-wave/boundary-layer interaction. In fig. 10 are compared the densities in the case of M = 3, Re = lo5 with those for M = 3, Re = 10’. There is only a minor difference owing to a difference in the grid systems. The interaction is still regarded as strong in the case of Re = 107. In fig. 11 the cases calculated with and without the turbulence model for A4 = 3, Re = 10’ are compared. There is no large difference between the two figures. Case (3): The calculation of the interaction of the reflected shock-wave/boundary-layer interaction, for Re = lo5 without the turbulence model. In this case, we discuss the relationship between the thickness fo the upstream boundary layer and the shock-wave/ boundary-layer interaction. The Mach number of the reflected shock wave relative to the upstream one is 2.2. Fig. 12 shows the density contours. The formed normal shock wave propagates (fig. 12a), reflects and propagates towards the opposite direction, then the interaction occurs (fig. 12b) and the shape of the reflected shock wave changes to X-type and X-type (fig. 12c-12d). Fig. 13 shows Schlieren photographs of the reflected shock-wave/boundary-layer interaction taken by North and Stuart (1962). Our calculated results are in very good agreement with their
Fig. 4. Schlieren
photograph
by Tamaki
et al
I. Hataue / Shock-wave/boundary-layer
interaction
225
h
Pmin Fig. 5. Comparison of computed by Tamaki et al. (b). (- - -)
wall and interior pressure fluctuation at T = 12 (a) with experimental determined y/h =l.O, (- - -- --) y/h = 0.5, (.-.-.) y/h = 0.25, ( -) y/h=O.
ones
experimental data. Fig. 13b corresponds to fig. 12~. In this case, the upstream boundary layer becomes thicker through the influence of the propagating incident shock wave. As a consequence, the amplitude of the reflected shock-wave/boundary-layer interaction becomes larger and the X-type shock wave is formed faster. Case (4): The calculation of the interaction of propagating shock-wave/ boundary-layer interaction for Re = 105. Comparing with the results of cases (2) and (3), we discuss here the relationship between the parameters (the strength of the shock wave and the thickness of the upstream boundary layer) and the shock-wave/boundary-layer interaction again. Fig. 14 shows the density contours. The Mach number of the shock wave (MS) in this case is about 2.6. The shock wave proceeds, firstly the bifurcation occurs, then the shape changes to X-type and the second shock wave (pseudo-shock wave) is formed. We computed the flows for several Mach numbers of the shock wave. When the Mach number of the shock wave is less than 2.4, bifurcation of the normal shock wave never occurs and the bifurcation begins very late for MS = 2.6. That is the reason why the boundary layer is very thin, and therefore the interaction is weak. 3.2. Three-dimensional
case
Two kinds of three-dimensional calculation are performed under the same condition case (2) of the two-dimensional case (M = 3, Re = 105, without the turbulence model).
as the
I. Hataue / Shock-wave/boundary-layer
226
incident
interaction shock -wave
~
Fig. 6. Density
contours
at M = 2, Re = 10’. without
Fig. 7. Density
contours
turbulence
at M = 3, Re = lo’, without
model and pressure
turbulence
fluctuation.
model.
(b) Fig. 8. Schlieren
photograph
by Ikui et al. at Re = IO’. (a) M = 2.15, (b) A4 = 2.68.
I. Hatme / Shock-waae/boundary-la~ver interaction
227
(4
(b) Fig. 9. Density contours at M = 3, Re = 105, without turbulence model. (a) 350 X 77 grid points (same calculations as fig. 7 but at different nondimensional time), (b) 350x 153 grid points. (b) 350 x 153 grid points.
Case (5): Duct with rectangular cross section. Fig. 15 shows the density contours in two different horizontal planes. Comparing the figure of the central plane (fig. 15b) with that of the two-Dimensions case (fig. 15c), we find that there is no large difference. In fig. 16 shows the velocity vectors are plotted, projected on each cross section superimposed on the density contours. The patterns of the secondary flow induced by the shock-wave/ boundary-layer interaction are different for the upstream side and the downstream side of the front oblique X-type shock wave (figs.16c, 16e). In the upstream region of the front oblique shock wave, the
(4
fbf Fig. 10. Density contours at (a) M = 3, Re = lo5 (b) M = 3, Re = lo’, without turbulence model.
228
1. Hatme / Shack-waoe/boundat-iuy~r
interartiort
(b) Fig. 11. Density
contours
(a) without
turbulence
model (b) with turbulence
model at M = 3, Re = 105.
(a) T-5.0
(b) Tz5.7
(c) T=6.4
(e) T=7.9
Fig. 12. Density
contours.
(a) Incident
shock wave, (b)-(e)
reflected
shock wave,
I. Hataue / Shock-waoe/boundary-layer
229
interaction
(4 Fig. 13. Schlieren
photograph
of reflected
shock-wave/boundary-layer
interaction.
secondary flow occurs from the corner towards the center along the side wall in the separated boundary layer. When the direction of the main flow is toward the wall behind the oblique shock wave, the flow in the cross section is in the radial direction. Fig. 17 shows the
I--
(a) T=4.8
incident
propagating
shock-wave
-
(b) T=5.5
Fig. 14. Density
contours
in the case of propagating
shock wave.
I. Hatme
230
/ Shock-wave/boundary-lacer
interactron
(b)
Fig. 15. Density plane.
contours
in three-dimensional
case (a, b) and the two-dimensional
case (c). (b) is the one of the central
instantaneous streamlines projected to each cross section. Many vortices appear in the separated boundary layer. Case (6): Pipe with circular cross section. Fig. 18 gives density contours in a plane of symmetry and velocity vectors projected on each cross section superimposed on the density contours. Fig. 19 shows the instantaneous streamlines in each cross section. The density contours are almost axially symmetrical. Between the density contours in the symmetrical plane and that in the central horizontal plane in the case of a rectangular duct (fig. 15b), we find no large difference. However, the streamlines in the cross section are almost axially symmetrical.
I. Hataue / Shock-wave/boundary-layer
interaction
231
fdl (40/75) Fig. 16. Plots of velocity vectors in each cross section mth grid plane numbered from the inlet of the total
superimposed on the density n grid planes.
(25175)
Fig. 17. Streamlines
in each cross section.
contours.
(e) (W75) Symbol
(m/n)
means the
I. Hataue / Shock-wave/boundary-layer
232
(20/175)
interaction
(60/175)
(80/175)
(140/175)
(120/175) Fig. 18. Density contours in the plane superimposed on the density contours.
of symmetry
(top figure),
(!60/175) and plots of velocity
vectors
in each cross section
4. Conclusion The shock-wave/boundary-layer interaction was successfully computed on the basis of twoand three-dimensional Navier-Stokes equations by the shock capturing method (Harten-type second-order TVD scheme) with LU factorization. The mechanisms of bifurcation and formation of X-type shock wave, X-type shock wave, and pseudo-shock wave were captured. By calculating several cases, the dependence of the flow pattern on the upstream Mach number and the upstream boundary-layer thickness was clarified. Furthermore, the secondary flow showing the three-dimensionally of the flow is clearly observed.
I. Hataue / Shock-waue/boundary-layer
(110/175)
(120/175) Fig. 19. Streamlines
interaction
(140/175)
233
(160/175)
in each cross section.
Acknowledgement
The author is grateful to Professor Hideo Takami of the University of Tokyo for his valuable discussion and to Professor Kunio Kuwahara of ISAS (Institute of Space and Astronautical Science) and Dr. Katsuya Ishii of ICFD (Institute of Computational Fluid Dynamics) for their kind suggestions.
References
Baldwin, B.S. and H. Lomax (1978) AIAA paper 78-0257. Beam, R. and R.F. Warming (1976) J. Comput. Phys. 22, 87-110. Chakravarthy, S.R. and S. Osher (1985) AIAA paper 85-0363. Crocco, L. (19581 Fundamentals of Gas Dy~amics~ Section B {P~nceton) p. 124. Fujii, K. and S. Obayashi (1986) AIAA paper 86-1831. Harten, A. (1984) SIAM J. Numer. Anal. 21, l-23. Ikui, T., K. Matsuo and M. Nagai (1973) Trans. Japan Sot. Mech. Eng. B 39-326. 3054-3064. Lombard, C.K., J. Bardina, E. Venkatapathy and J. Oliger (1983) AIAA paper 83 - 189SCP. MacCormack, R.W. (1969) AIAA paper 69-354. Merkli, P.E. (1976) AIAA J. 14, 168-172. Neumann, E.P. and F. Lustwerk (1949) J. Appl. Mech. 16, 195-202. North, R.J. and C.M. Stuart (1962) Nat. Phys. Lab. Aeron. Rept. 1029. Obayashi, S. and K. Kuwahara (1984) AlAA paper 48.1670. Om, D. and M.E. Cbilds (1985) AIAA J. 23, 1506-1511. Pulliam, T.H. and D.S. Chaussee (1981) J. Comput. Phys. 39, 347-363. Pulliam, T.H. and J.L. Steger (1985) AIAA paper 85-0360. Roe, P.L. (1981) J. Compuf. Phys. 43, 357-372. Steger, J.L. (1978) AIAA J. 17, 249-255. Tamaki, S., Y. Torn&a and R. Yamane (1969) Trans. Japan Sot. Mech. Eng. B 35-273, 1028-1035
234
I. Hataue / Shock-wave/boundary-layer
interaction
Tamaki, S., Y. Tomita and R. Yamane (1970) Trans. Japan Sot. Mech. Eng. B 36-292, 2056-2066 Viegas, J.R. and T.J. Coakley (1978) AIAA J. 16, 293-294. Viegas, J.R. and CC. Horstman (1979) AIAA J. 17, 811-820. Waltrup, P.J. and F.S. Billig (1973) AIAA J. Il, 1404-1408. Waltrup, P.J. and J.M. Cameron (1974) AIAA J. 12, 878-880.