A Chebyshev-Fourier discretization with shock fitting is used to solve the unsteady Euler equations. The method is applied to shock interactions with plane waves and with a simple model of homogeneous isotropic turbulence. The plane wave solutions are compared to linear theory

Hussaini, M. Y.

Kopriva, D. A.

Zang, T. A.

English

NASA Technical Reports Server

NASA Contractor Report 17213 NASA-CR-172133 
19830019681 
leASE 
PSEUDOSPECTRAL CALCULATION OF SHOCK TURBULENCE INTERACTIONS 
T. A. Zang 
D. A. Kopriva 
M.Y. Hussaini 
'Grant No. NAG1-109 . 
Contract Nos. NASl-l6394, NASl-l7070, NAS1-17130 
May 1983 
INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 
Operated by the UniverSities Space Research Association 
NI\S/\ 
National Aeronautics and 
Space Administration 
langley Research Center 
Hampton, Virginia 23665 
I IUIIIII IIU IIII UIII IIUI IIIII 11111 IIII lUI 
NF01585 ~13RARV COpy 
LANGLEY -
Ll6R.t.R 
HN.~PTON. VIRGINIA 
https://ntrs.nasa.gov/search.jsp?R=19830019681 2020-03-21T02:39:43+00:00Z
PSEUDOSPECTRAL CALCULATION OF SHOCK TURBULENCE INTERACTIONS 
* T. A. Zang 
College of William and Mary 
** D. A. Kopriva 
ICASE 
** M. Y. Hussaini 
ICASE 
ABSTRACT 
A Chebyshev-Fourier discretization with shock fitting is 
used to solve the unsteady Euler equations. The method is 
applied to shock interactions with plane waves and with a 
simple model of homogeneous isotropic turbulence. The plane 
wave solutions are compared to linear theory. 
* Research Associate at the College of William and Mary, 
Williamsburg, VA 23185, NASA Grant No. NAG1-109. 
** Staff Scientist and Senior Staff Scientist, respectively, at 
the Institute for Computer Applications in Science and 
Engineering at NASA Langley Research Center, Hampton, VA 
23665, NASA Contract Nos. NAS1-16394, NAS1-17070 and NAS1-
17130 • 
i 
1. INTRODUCTION 
The effects of a shock wave upon a turbulent boundary· 
layer have been studied extensively by both experimental and 
computational means. A persistent experimental result is that 
the turbulence levels downstream of the shock are higher, by 
at least a factor of 2, than the upstream levels {l]. The 
computational results, however, have consistently failed to 
exhibit this dramatic turbulence amplification [2]. This 
failure is almost surely due to inadequacies in the turbulence 
models used in the computations. Current turbulence models 
have achieved their greatest success in low-speed flows and 
attached high-speed flows which agree well with experiment 
even up to Mach 20 [3]. Evidently, some essential physics of 
separated compressible flow is missing from the models. Their 
improvement must await a better understanding of the 
interaction between shock waves and turbulence. 
Two theoretical tools are now available for the study of 
the inviscid aspects of this interaction. The basically 
linear effects may be examined analytically in a way developed 
by Ribner [4] and others, and used by Anyiwo and Bushnell [5] 
in their assessment of turbulence amplification in the shock-
boundary layer interaction. Numerical computations based on 
the nonlinear, two-dimensional Euler equations are more 
recent. Pao and Salas [6] used this approach to study the 
shock-vortex interaction. So did Zang, Hussaini and Bushnell 
[7] in their examination of the interaction of plane waves 
with shocks. Both sets of numerical computations used second-
order finite differences for the spatial discretization. In 
this paper we will describe a pseudospectral spatial 
discretization. 
2. NUMERICAL METHOD 
The model problem which is used to study the turbulence 
amplification and generation mechanisms is illustrated in 
Figure 1. At time t = 0 an infinite, normal shock at x = 0 
separates a rapidly moving, uniform fluid on the left from the 
fluid on the right which is in a quiescent state except for 
some specified fluctuation. The initial conditions are chosen 
so that in the absence of any fluctuation the shock moves 
uniformly in the positive x-direction with a Mach number 
2 
(relative 
presence 
ripples. 
to the fluid on the right) denoted' by Ms. In the 
of fluctuations the shock front will develop 
The structure of the shock is described by the 
IN FLOW 
\ 
PHYSICAL PLANE 
Xs 
o' 
o 
o 
o 
o 
. 
----------, ' .. 
----------" 
, , , 
----------" "'1' 
--------- .... , .... ~ I, 'I I I 
---------"-~. "'1' 
, , .' , 
, . , 
, , 
, I , , •• 
, , , , 
_________ ~,-.,.,., ,,/  ""'0, 
_________ , ..... - ........ I I' ' I I ,. 'I 
-------- ....................... 'III """ 
-------- ......... _ ............. I11I ,," •. , 
________ .......... _ .............. '" II ""1 
-------- .... ,-~~- I' 
--------- ...... - ...... ~ I' 
, , 1 , 
, , , •• , " I I 
• I I I 
_________ ..... - ...... - , , , I •• , ./ I ' 
.. _--------_ ...... ",,-. I • 
, , , 
, , , I 
I I I 
o 
o 
o 
---------_ ...... -... 
----------",-, 
----------", .... \. 
-----------'" 
-----------" ..... 
o 
o 
o 
, , I , 
, I I , •• , I , • 
"""111 
, , " , , I I 
, , , , 
, I I I 
. , 
.. , 
, , 
LEFT COMPUTATIONAL 
BOUNDARY 
SHOCK 
WAVE 
DOWNSTREAM 
FLOW 
Figure 1: Model problem in the physical domain. 
3 
function Xs (y, t). The numerical calculations are used to 
determine the state of the fluid in the region between the 
shock front and some suitable left boundary, xL{t), and also 
to determine the motion and shape of the shock front itself. 
In the applications given below the physical problem is 
periodic in y with period YL. The change of vari~b1es 
x = xs(y,t) - xL(t) 
Y = Y/YL 
T t 
produces the computational domain 
OC;XC;l 
o " Y " 1 
T ) O. 
(1) 
(2) 
In terms of the computational coordinates the two-dimensional 
Euler equations are 
(3 ) 
where Q = (P, u, v, S) , 
u yX yX 0 
c2 
x Y 
-x U 0 0 y x 
B = 
c2 
(4) 
-x 0 u 0 y y 
0 0 0 U 
and 
V yY yY 0 
c2 
x Y 
-Y V 0 0 y x 
C = • 
c2 
(5 ) 
-Y 0 V 0 
Y Y 
0 0 0 V 
The contravariant velocity components are given by 
and (6 ) 
V = Yt + uYx + vYy• 
A subscript denotes partial differentiation with respect to 
4 
the indicated variable. P, c and S are the natural 
logarithm of pressure, the sound speed and the entropy 
(divided by the specific heat at constant volume), 
respectively, all normalized by reference conditions at 
downstream infinity; u and v are velocity components in 
the x and y directions, both scaled by the characteristic 
velocity defined by the square root of the pressure-density 
ratio at downstream infinity. The ratio of specific heats is 
denoted by y; a value y = 1.4 has been used for all the 
calculations in this paper. 
Let n denote the time level and ~t the time 
increment. The time discretization of Eq. (3) is 
(7) 
(8) 
where L denotes the spatial discretization of Bd X + Cd y • 
The solution Q has the Chebyshev - Fourier series expansion 
Q(X,Y,T) = 
M N/2-1 I I Q (T)T (~)e2~iqy, 
p=o q=-N/2 pq P 
(9 ) 
where ~ = 2X-l and Tm is the Chebyshev polynomial of 
degree m. The derivatives QX and Qy are approximated by 
M N/2-1 
Q = 2 I I Q(I,O)(T)T (~)e2~lqy (10) 
X p=o q=-N/2 pq P 
M N/2-1 
Qy = 2~ I I Q(O,I)(T)T (~)e2~iqy (11) p=O q=-N/2 pq P 
where 
Q(1,O) 2 M 
=- I m ~q pq c 
m=p+1 p 
(12) 
p+m odd 
with 
{2 p = ° c = p 1 p ;> 1 (13) 
and 
Q(O,I) 
pq = iq ~q. (14) 
The most critical part of the calculation is the 
treatment of the shock front. The shock-fitting approach used 
here is desirable because it avoids the severe post-shock 
oscilla tions that plague shock-capturing methods. The time 
derivative of the Rankine-Hugoniot relations provides an 
equation for the shock acceleration. This equation is 
integrated to update the shock position. The jump conditions 
themselves are then used at xs(y,t) to determine the 
variables on the left (high perssure) side of the shock from 
the specified downstream flow (see [8] for details). For 
Ms > 2.08, the left boundary at XL(t) is a supersonic inflow 
boundary. Lower Mach numbers require more sophisticated 
boundary conditions as discussed in [8]. Spectral methods for 
compressible flows typically require some sort of filtering. 
3. INTERACTION OF PLANE WAVES WITH A SHOCK 
The non-linear interaction of plane waves with shocks was 
examined at length in [7]. The numerical method used there 
was similar to the one described above but employed second-
order finite differences in place of the present 
pseudospectral discretization. Detailed comparisons were made 
in [7] with the predictions of linear theory [41. The linear 
results turned out to be surprisingly robust, remaining valid 
at very low (but still supersonic) Mach numbers and at very 
high incident wave amplitudes. The only substantial 
disagreement occurred for incident waves whose wave fronts 
were nearly perpendicular to the shock front. This type of 
shock-turbulence interaction is a useful application of the 
pseudospectral technique. The method can be calibrated in the 
regions for which linear theory has been shown to be valid. 
The questionable regions can then be gainfully investigated by 
means of this alternative method •. 
We concentrate on Mach 3 normal shocks interacting with 
incident vorticity waves. The perturbation velocities ahead 
of the shock are taken to be 
(u' ,v') = A'( -k1,y,k1,x)(Cl!Plkl)cos(k1·.x), 
where the wavevector kl - (k1 x,k1 y) and A' is the 
amplitude. In terms of the incidJnce angle 81, kl-(k1cos 81,k1sin 81). The transmitted vorticity behind the 
shock can be expressed in the same manner with all 
subscripts 1 changed to 2. The. vorticity transmission 
coefficient is defined to be the ratio A2!Ai. At Mach 3 the 
ratio of the rotational part of the perturbation velocities 
is 0.158 times the transmission coefficient; if these 
perturbation velocities are scaled by the respective mean 
velocities (in the frame in which the shock is at rest), then 
this ratio is 0.611 times the transmission coefficient. 
A typical numerical simulation is shown in Figure 2. At 
time taO all the fluid in the computational domain, which 
then extended from xL = -0.5 to Xs = 0.0, was in uniform 
5 
6 
Figure 2: 
VORTICITY 
ENTROPY 
Responses at 
vorticity wave 
30° angle of 
!:I :: 
! 
: ; 
'1' II 
:: :: 
: : ;: 
... 
I I If 
: ::. 
! t!; 
I ", . . 
• I 
I • 
! { 
: 
: 
: 
. { 
• I: 
... 
: I 
I • 
· . 
· ,= 
! 
i 
: 
! 
: 
l 
: 
! 
• 
· .1 \: 
time t 
striking a 
= 0.80 
Mach 3 
incidence. The 
PRESSURE· 
to a 0.1% 
shock at a 
shock front is 
denoted by the solid line connecting the top 
and bottom boundaries. Dashed lines indicate 
negative contour levels. 
7 
flow. The flow ahead of the shock was given by Eq. (15) with 
an amplitude A"o = 0.001. In order to reduce the transients 
that would result if the fluid were to encounter this wave in 
a sudden fashion, the amplitude was turned on smoothly from 
A"o = 0.000 to A"o = 0.001 during the passage of the shock 
over the first half-wavelength of the incident wave. The 
interaction of the shock with the incident vorticity wave has 
produced a transmitted vorticity wave as well as generated 
acoustic and entropy waves. The acoustic wave clearly has a 
larger inclination angle and a faster propagation speed than 
the other two waves. 
e 
00 
0«> 
o o 0 
6 
<> 
2 o 
0~ ____________ ~ ___ 0~ __ 0 __ ~ ___________ L-__________ ~ 
1.0 1.5 2.0 2.5 3.0 
X 
Figure 3: Post-shock dependence at time t - 0.80 of the 
vorticity response to a vorticity wave incident 
at 30° to a Mach 3 shock. The solid line is 
the linear theory prediction. Computed 
responses are given for 0.1% (circles) and 
50% (diamonds) waves. 
Figure 3 indicates how we measure the transmitted 
vorticity. The computed velocities are differentiated 
pseudospectrally and combined to produce the response 
vorticity at each point on the computational grid. At each 
fixed value of X we perform a Fourier analysis in Y of the 
vorticity. . The wavenumber of interest is q = 1 (see Eq. 
(9». This amplitude is then scaled as indicated above to 
yield AZ and thence the ratio AzI Ai which is plotted in 
Figure 3 as a function of a typical value of the physical 
coordinate x corresponding to X. In the case shown in 
Figure 3 the incident wave does not reach full strength 
8 
until t = 0.56. A single value representing the transmission 
coefficient is obtained by averaging the x-dependent 
responses. The standard deviation of these values serves as 
an error estimate, indicated by the error bars in the 
following two figures. 
The angular dependence of the vorticity transmission 
coefficient for Mach 3 is displayed in Figure 4. Correspond-
ing results from the finite difference discretization used in 
[7J are shown for comparison. The amplitude dependence for 
30° incident waves is shown in Figure 5. 
18 
15 
IJ 
o 
IJ 
o 
3 
o~ __ ~ __ ~~ __ ~ __ ~ ____ ~ __ ~ ____ ~ __ ~ __ ~ 
o ~ ~ ~ 
INCIDENT ANGLE 
Figure 4: Incident angle dependence of the vorticity 
response to O.U: vorticity waves incident 
upon a Mach 3 shock. Pseudospectral (circles) 
and finite difference (squares) results are 
given. The solid line is the linear theory 
prediction for sub-critical angles. 
Linear theory predicts that the vorticity responses are 
peaked in the vicinity of 91 = 62°, the so-called critical 
angle. At larger angles the acoustic responses are 
evanescent. These predictions do not appear in the non-linear 
calculations. Tentative explanations for this failure of 
linear theory are provided in [7J. 
Below the critical angle, however, the linear theory is 
extremely robust. Figure 5 indicates that it remains valid 
for incident relative velocity fluctuations aproaching 50%. 
.. I 
9 
(In the frame in which the shock is at rest, the ratio of the 
fluctuating velocity to the mean velocity is one-third of the 
amplitude Ai.) 
>-
I--
t-I 
8 
7 
r-
~( 
~ 
~6 
--I--
a::: 
a 
:> 
5 -
q 
o 
~ () 
Figure 5 
I I I I 
-
- -
~ 
-
(~ 
() 
~ 
-
-
-
0 
I I I 
.5 1.0 1.5 2.0 2.5 
AMPLITUDE 
Amplitude dependence of the vorticity response 
to vorticity waves incident at 30° to a Mach 
3 shock. The solid line is the linear theory 
prediction. 
4. INTERACTION OF TURBULENCE WITH A SHOCK 
Low intensity turbulence can be viewed as the 
superposition of plane shear waves just discussed. Figure 6 
shows velocity vector plots of the interaction of a shock with 
a simple homogenous isotropic turbulence model in which the 
velocities are defined as sums of sinusoidal waves with random 
phases. The frame at t = 0.38 shows the undisturbed 
velocity field downstream of the shock and t = 1.6 shows the 
resul t after the shock has traveled approximately one 
wavelength of the longest wave component of the undisturbed 
turbulence. The development of the post-shock vortices seen 
in the figure coincides with the generation of rows of 
pressure and entropy spots behind the shock. 
5. CONCLUSION 
We have solved the two-dimensional unsteady non-linear 
Euler equations with a Chebyshev-Fourier pseudospectral 
10 
method. The linear theory of the interaction of a shock and 
plane shear waves has provided a test of the method. Details 
of the smoothing methods and questions of accuracy will appear 
in a later paper. But accuracy comparable to the finite 
difference calculations can be obtained with about half the 
number of Chebyshev modes in X and with only four modes in 
Y. The results indicate that it should be possible to compute 
more complicated shock-turbulence interactions. 
t t t . -~ 
~ : ; I .. t 
• 
, I , 
-
\ 
.. • 
, , 
-" 
. .. 
-
. 
· . ~ . 
t , , • 
" 
, 
t 1 ... ... ~ 
" t , 
--- ... 
• ..... 
· .......... 
· 
· 
",-~++;., 1" ..... t1,-' 
" .. _"" ..... 
-- .. ~ .. --
..---- .. ,. 
~ ........ _ ~ t ,,_ .. , " • -. -.. .... _ ~ , , .. -" 
• ~ .." t , , .. , • • , , " - I I I , • , 
• • ~ "'" t , , .. ~ • • , , , - " 1 I , • , 
..... , .... , t " ....... _- "., ........... '\t .. '-' 
........ · _ ........ - - -....-......- ....... -
_ # • " ..... , , , .... ____ .. _ ---. ....... -..~J' t ... _ , 
,1\"1' \',,,',., \o ... _.I1t\_~ ~ t , I • 1 , l ' .. , , • I I ~ - \ t f., 
t , ., ... \0 ~ , _ - .... " • ,.,.,._ .... '\ t ~ ..... 
. .". - ..... , , ..... - - ...... .....-....... .". - .. 
• • r- • ..... 
" . • • •• 
· _ ....... -
· , .... . 
- _ .......... 
. ,,----, I , • - .. 
. 
t 
. • .. • ..... .. t \ 
- • • • . .. ....... I' 
.... , •. -_.-.,r--
•••• -"""""li"-
.......... ""!-: 
- ••• , \ , I 
...... , .... , 
· ... ~ ... -........... ' 
•• ~ •• 4~.".'·, 
........ ,"'".,t, 
........... ,',/tl'-' 
--__ 'or- -
. . .. 
· .. .. 
-_ .... -..... ~ '.I{".,,~4 
•• t '1'\ 
.. .. "-". . "'-
· ... .................. .. ~~! 
· .. "" .".""--
... .... ,.,., 
-I-
-I-
., 
Figure 6: Velocity vector plots at t = O. 38 (top) and 
t = 1.6 (bottom) of the interaction of a Mach 
3 shock with low intensity homogeneous 
isotropic turbulence. 
... 
REFERENCES 
1. ROSE, W. C. - The Behavior of a Compressible Turbulent 
Boundary Layer in Shock-Wave-Induced Adverse Pressure 
Gradient. NASA TN D-7092, 1973. 
2. HORSTMAN, C. C., HUNG, C. M., SETTLES, G. S., VAS, I. E. 
and BOGDONOFF, S. M.- Reynolds Number Eff.ects on Shock-
Wave/Turbulent Boundary-Layer Interactions - A Comparison 
of Numerical and Experimental Results. AIAA Paper, 77-
42, 1977. 
3. BUSHNELL, D. M. CARY, A. M., JR., and HARRIS, J. E. -
calculation Methods for Compressible Turbulent Boundary 
Layers. NASA SP-422, 1976. 
4. RIBNER, H. S. - Shock-Turbulence Interaction and the 
Generation of Noise. NACA Report 1233, 1955. 
5. ANYIWO, J. 
Amplification 
Interactions. 
C. and BUSHNELL, D. M. Turbulence 
in Shock Wave Boundary Layer 
AIAA Journal, 20, pp. 893-899, 1982. 
6. PAO, S. P. and SALAS, M. D. - A Numerical Study of TWo-
Dimensional Shock Vortex Interaction. AIAA Paper, 81-
1205, 1981. 
7. ZANG, T. A., HUSSAINI, M. Y. and BUSHNELL, D. M. 
Numerical Computations of Turbulence Amplification in 
Shock Wave Interactions, AIAA Journal, to appear. 
8. HUSSAINI, M. Y., SALAS, M. D. and ZANG, T. A. - Spectral 
Methods for Inviscid Compressible Flows in Advances in 
Computational Transonics. Ed. Habashi, W. G., Pineridge 
Press, Swansea, United Kingdom, 1983. 
I 2. Government Accession No. 3. Recipient's Catalog No. 1. Report No. 
NASA CR-172133 
4. Title and Subtitle 5. Report Date 
May_ 1983 
PSEUDOSPECTRAL CALCULATION OF SHOCK TURBULENCE INTERACTION 6. Performing Organization Code 
7. Authorls) 
T. A. Zang, D. A. Kopriva, M. Y. Hussaini 
8. Performing Organization Report No.'. 
83-14 
I---------------------------------i 10. Work Unit No. 
9. Performing Organization Name and Address 
Institute for Computer Applications in Science 
and Engineering NAG~"1'a~ ,~A'!il~:&394 '. 
Mail Stop 132C, NASA Langley Research Center NASt-17070, NASl-17130 
~~H~a=~~·t~o~n~-V~A~~2~3~6~6~5~--------------------~11 Ty~~ReportmdP~~COver~ 
12. Sponsoring Agency Name and Address contractor report 
National Aeronautics and Space Administration 
Washington, DC 20546 
15. Supplementary Notes 
Technical Monitor: Robert H. Tolson 
Final Report 
16. -Abstract 
14. Sponsoring Agency Code 
A Chebyshev-Fourier discretization with shock fitting is used to solve the unsteady 
Euler equations. The method is applied to shock interactions with plane waves and 
with a simple model of homogeneous isotropic turbulence. The plane wave solutions are 
compared to linear theory. 
17. Key Words (Suggested by Author(s)) 18. Distribution Statement 
spectral methods, turbulence, shock waves 
19. Security Oassif. (of this report) 
Unclassified 
Unclassified-Unlimited 
Subject Category 01 
20. Security Classif. (of this page) 
Unclassified 
21. No. of Pages 
13 
22. Price 
A02 
N-305 For sa Ie by the Nationa I Technica I Information Service, Springfie Id. Virginia 22161 
.. 1 
-
End of Document 


Psuedospectral calculation of shock turbulence interactions

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server