Supersonic turbulent boundary layers over two-dimensional protuberances are investigated, using the numerical finite difference alternating direction implicit (ADI) method. The turbulence is modeled mathematically. The turbulence is represented here by the eddy viscosity approach. The turbulent boundary layer structure as well as an interest in thick boundary layers and much larger protuberance heights than in the laminar case lead to new difficulties. The problems encountered and the means to remove them are discussed

Polak, A.

Werle, M. J.

English

NASA Technical Reports Server

General Disclaimer 
One or more of the Following Statements may affect this Document 
 
 This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 
 
 This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 
 
 This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 
 
 This document is paginated as submitted by the original source. 
 
 Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 
 
 
 
 
 
 
 
Produced by the NASA Center for Aerospace Information (CASI) 
https://ntrs.nasa.gov/search.jsp?R=19760006955 2020-03-22T17:33:18+00:00Z
University of Cincinnati, Cincinnati,.Ohio 45221
NUMERICAL STUDY OF SUPERSONIC TURBULENT
FLOW OVER SMALL PROTUBERANCES
Semiannual Progress Report,
July - December 1975
V	 S;
^p1 R
NASA Grant No. NSG 1208;
(NASA-CR-145924) NUMERICAL STUDY OF 	 N76-14043
SUPERSONIC TURBULENT FLOW OVER SMALL
PROTUBERANCES Semiannual Progress Report,
Jul. - Dec. 1975 . (Cincinnati Univ.) 16 p HC 	 Unclas
$3.50	 CSCL 20D G3/02 0.6795 i
Principal Investigators are A. Polak and M. J. Werle,
Department of Aerospace Engineering.
The NASA Technical Officer for this grant is Mr. James
C. Dunavant, NASA Langley Research Center, Hampton,
Virginia 23665.
NUMERICAL-STUDY OF SUPERSONIC TURBULENT
FLOW OVER SMALL PROTUBERANCES
NASA Grant No. NSG 1208
Semi-Annual Status Report
•	 July 1, 1975 - December 31,-1975
by
A. Polak and M.J. Werle
University of Cincinnati, Cincinnati, Ohio
This report summarizes the effort expended during the first six
months on the study of supersonic turbulent boundary layers over two-
dimensional protuberances. The approach taken is to extend an earlier
work (Ref. 1) to the turbulent regime, using the numerical finite-
difference alternating direction implicit (ADI) method. An obvious
departure from the previous case (Ref. 1) is forced by the need to
model mathematically the turbulence. The turbulence is represented,
here by the eddy viscosity approach as used in Refs. (2) and (3).
The turbulent boundary layer structure as well as an interest in thick
boundary layers and much larger protuberance heights than in the
laminar case lead to new difficulties. The problems encountered in
the course of the present investigation and the means to remove them
will be discussed later in this report.
The Siverning equations pertinent to the problem are
Continuity
V  + F + 2EF C = 0
-2-
Momentum
. ae
(k7F n ) n — VFn + ( l+a/2) (S+ atT (g—F2 ) - 29FF^ = 0
Energy
(
keg
 P n) — Vg 2Rr	 [R(e—e/Pr) FF n ]^ — 2&Fg^ = 0
n
with
F (4,0) = V( ,0) = 0 . g (4,0) = Hw/He
F (C, -) = g (E,-) = 1	 .
Ths momentum equation is cast into the present form to facilitate the
use of the ADI algorithm. The notation is the same as in Refs. (1)
and (3). Hence, F = normalized streamwise component of velocity,
g = normalized total enthalpy, V = transformed normal velocity function,
AT = total displacement body, k = molecular viscosity parameter, a = An
A	 .
inviscid parameter, and a and a are eddy viscosity parameters defined
for fully turbulent flow as e = 1 + e /u, e = 1 + pr 	 where e =
T
eddy viscosity.	 The Prandtl numbers, Pr and Pr  are taken to be
constant.
The study commenced with programming the above equations using
the Cebeci-Smith two -layer eddy viscosity model, applied successfully
to attached turbulent boundary layers earlier (Refs. (2), (3)). A
* Initially, this choice is adequate since the immediate aim is to
develop an efficient numerical algorithm; comparison with experi-
mental'data will be delayed for later.
-3
-3-
base case was identified with flow conditions corresponding to the
NASA test data (Ref. (4)). The free stream Mach number and static
temperature are Mm = 2.5 and T= = 2520R respectively. The Reynolds
number based on free stream conditions and a reference length L =
15.25 cm is Re
.,
 = 1.647x10 6 . The sine-wave protuberance profile is
given by y = h {1 - cos (wn (x - xA)l} , where h is the height and
w is the width of the protuberance. In all the base case studies
w = 0.24 (i.e. w = 0.24 x 15.25 cm = 3.66 cm; star quantities are
dimensional). For a single wave placed on a-flat plate far downstream
from the leading edge, at s = 25.85, and h = 0.02 a calculation was
performed with Tw/To = 0.81, Pr = 0.72 and Pr  = 0.9. An attached
solution was obtained (Fig. 1). In this case the boundary layer
displacement thickness 6 is about three times the protuberance height
(but only about 1/3 as thick as the NASA Langley tunnel wall boundary
layer). The calculation commenced at s  = 1.0 with a laminar boundary
layer, passing through transition and becoming fully turbulent at
s = 1.4. After sweeping only once in C direction up to station so =
25.5 the ADI algorithm is employed between stations s o and s  = 26.5.
The protuberance is placed between x = 25.85 and 26.09. The stream-
wise step size, As was taken 0.02 in this calculation so that there
are at least 12 grid points over the protuberance surface.- In 6.7
minutes of computer time 22 time sweeps were performed and the skin
friction is varying (between the last two time steps at a given station)'
in the third or higher decimal place. This was considered a converged
solution. It is worthwhile noticing that in the laminar case (Ref. 1)
the same number of time steps required about 15-20 minutes computer
time. The present efficiency was accomplished by linearization of the
-4-
coefficients of the corresponding difference equation around the
previous station value, eliminating thus the need to update these
coefficients in subsequent iterations.
Increasing the protuberance height at otherwise identical flow
conditions as given above, a separation region appeared ahead of the
protuberance. But the values of skin friction (the parameter we
monitor for testing the quality of the solution) did not converge in
time. On the contrary, wide oscillations occured in time indicating
some numerical instability, and eventually the calculation terminated
abruptly. Examination of the boundary layer velocity profiles showed
slight oscillations (of order 10-6 around the value F = •1) in the
boundary layer edge region. At this point it was not clear whether
this was due to a stronger interaction (higher protuberance), the
reverse flow region, an improper eddy model or some other effects.
To isolate the source of trouble it was decided to simplify the problem
by taking Pr = Pr  = 1, thus eliminating the need to solve the energy
equation. (This corresponds to adiabatic wall condition with g =
constant = 1); the rationale being that since the artificial time de-
pendent term appears only in the momentum equation, the source of the
time-like oscillations maybe due to finite-difference representation
of this equation. Because the difficulty appeared for the separated
cases, the eddy viscosity model became suspect as well. 'Freezing ,
 the
eddy viscosity distribution (this idea was borrowed from Ref. 5) near
the forward plate-protuberance junction, and the use of Alber's model
(Ref. 6) for the separated region did not resolve the difficulty.
It was brought to our attention (Ref. 7) that upwind differencing of
the F n term in the momentum equation maybe required to satisfy the
^I
-5-
convergence criteria of the numerical scheme (see also Ref. 8). In
the present case the momentum equation, after linearization, may be
written as
Fnn+ alFn + a 2 + a3 + a4FE + a5 V = 0
where al ,	 a5 are the linearized coefficients. In the boundary
f'	 layer near the wall the diffusion term 
Fnn 
is predominant over the
convection-like term a1F n and a central difference scheme for Fn
is appropriate. From a so-called diagonal dominance test of the model-
equation F 	 aF n
 = 0 it is found that with central differencing the
criteria IaAnJ < 2 must be met. Because of additional terms in.the
momentum equation we have found (by numerical experimentation) that
the use of the criteria JaAnJ < 1 eliminated the oscillations in the
profile near the boundary layer edge. Hence, F n
 is central differenced
when (aAn) < 1 and upwind differencing is used when (aAnJ > 1. Davis
(Ref. 7) .
 also pointed out that the use of forward differencing of F
term in the continuity equation is more appropriate because the upstream
propagation becomes important for the strongly interacting boundary
layers. Having incorporated these ideas into our algorithm difficulties
mentioned above have been removed; the boundary layer profile approaches
smoothly the edge values (F = 1) and the time-like oscillations also
J
disappeared. For the base flow conditions, with the forward junction
of the plate-protuberance placed at s  = 3.35 and protuberance height
of 0.024 results are shown in Figures 2 and 3.' The dashed line in
Figure 2 represents the variation of wall shear with time steps at a
fixed station (s = 3.33) located just downstream of the plate-protu-
berance junction. This typically unstable solution, obtained before
•	 _6_
the above mentioned changes were introduced, terminates abruptly after
ten time steps. After the modifications in the numerical algorithm
were made an apparently convergent solution was obtained (full line,
Fig. 2). The variation in C  appears only in the third or higher
decimal place after the first fifteen time steps. The spatial dis-
tribution of C  at the last time step (33rd sweep) is plotted in Figure
3. A small separation bubble appears near the forward plate-protuberance
junction and the maximum shear is about double that of the flat plate
value. The boundary layer displacement thickness ahead of the protu-
berance is relatively thin, d ti 0.010 (or d * = SL*	0.1525 cm) and
d/h = 0.4.
A major interest of this research project is in comparisons of
the numerical predictions with experimental data of very thick separated
turbulent boundary layers produced on the tunnel wall (Ref. 4). The
separation characteristics over the test plate placed in the wind tunnel
wall are strongly dependent on the approaching turbulent boundary layer
profile. For this reason calculations were performed for the ••_)undary
layer as it develops along the wall of the UPWT Langley Tunnel. The
following test section free-stream conditions were taken: M
.,
 2.535,
Re
,,
/cm = 1.08x10 5 , To = 5670R (3150K), T. = 252oR, Tw/To = 0.81,
Pr = 0.72, PrT	0.9. A reference length L = 15.25 cm (* 05 ft).
Correspondingly, the reference Reynolds number is Re = p= .=
u .,
1.08x10 5x15.25 = 1.647x10 6 ; . and the nondimensional protuberance width
w of 0.24 corresponds to w* of 3.66 cm. The boundary layer calculation
was carried out as a non-interacting 2-D laminar-transitional-turbulent
boundary layer developing from ahead of the nozzle throat under a
favorable pressure gradient. In the supersonic region downstream of
the throat the pressure distribution was calculated from the sidewall
f-7-
Mach number distribution, using isentropic relations. The Mach
number distribution was obtained from a characteristic net (Ref. 9).
From these data a cubic representation for M was assumed of the form
2
M=MT -K1 (1-s-) ( 1 -K2 s )	 ,
2	 T
where MT - test section Mach number at s = s T. The above polynomial
representation satisfies the condition dM/ds = 0 at s = s T . At the
location s of the first characteristic near the throat the Mach number
a
was estimated to be Ma = 1.11. The measured distance along the wall
from s  to s T
 is 39.6 , (-. 20 ft). At a location s s 8.86 units downstream
of s  the Mach number is M S = 1.84. Letting s a = 1, s o = 9.86 and s 
40.6. Using these values, the constants K1 and K2 can be determined.
For K1 , K2 > 0 and K2 < 1 (which is the case here) the polynomial rep-
resentation given above yields monotomically increasing Mach number.
The subsonic-transonic section of the tunnel was also assumed to be
represented by a cubic polynomial
2
M = Mi +C1
 (i - s ) (1-C2 s )
with the following properties:
At s = si = 0.3, M Mi and dM/ds = 0.
At s = s a , M = Ma = 1.11 and ^ =
s
a+
Three different values for Mi were chosen: 0.01, 0.03 and 0.05.
It turned out that the boundary layer development downstream of the
throat vas not sensitive to these initial values. Taking M i = 0.03
i
liiitG a:vtulrsc%.a LOVMiNi3LY LCLY=A. Yivlllaa L+utat;.tcu vu 1.tXt vQLUD tJCLWtlCtt
locations 40.6 and 72.6. (At s = 40.6 the Mach number becomes con-
stant on the side wall; at s' ti 50.6 the straight section begins).
These profiles are presently being used as initial data for separation
studies of thick turbulent boundary layers. The table belcw gives the
calculated boundary layer displacement thickness distribution along
the constant Mach number section of the UPWT Langley tunnel wall at
ten stations
s	 4.0..6 .44.6. 48..6 . 52.6...5.6..6 .6.0.6 	 64.6	 68.6 72.6	 76.6
d (cm)	 .1.69 1.84 1.98. .2.12 ..2..26 ..2.40. 2.54 .2.67 . 2.80
	 2.93
Using the profile produced at station s = 72.60 and placing a
single sine-wave protuberance on the flat witti its junction at s  =
73.25 calculations were initiated to obtain separation characteristics
of thick turbulent boundary layers. The calculation starts at s = 72.60
by sweeping once to station 72.96 and then the ADI algorithm was employed
between this station and s = 73.96 using 51 grid points in the s direct-
ion ' (Os = 0.02) with 105 point variable mesh in n direction. Figure
4 shows results obtained for the base case flow conditions with h =
0.08, w = 0.24, Pr = 0.72, Pr  = 0.90 and T w/To = 0.81. Note that here
the energy equation was solved simultaneously with the continuity and
momentum equations. Figure 4 shows the variation of the surface skin
friction and pressure levels predicted by the present method. It is
seen that ahead of the protuberance the skin friction drops and a
small separation region develops on the front face of the wave. The
-9-
skin friction veaks just aft of the wave top with a value three times
the undisturbed flat plate value. The pressure distribution is some-
what different and shows a monct:omic increase up to a peak near the
top of the wave followed by a mild recovery back to the flat plate
value.
While the above results demonstrv,te the capability of the present
method to handle flow fields typical of present interest, there re-
main several refinements to by made before the technique can be con-
sidered iperational.
The first of these involves the siza of the finite difference
fresh applied to the problem. Using the present mesh from s = 73 to
74 it is noticed that the surface pressure distribution for all points
in the mesh ahead of an aft of the protuberance is influenced by the
protuberance. This effect is manifest by the depressed pressure levels
(p/pm < 1) near s = 73 and the monotonically increasing levels of C 
ahead of protuberance. A test was conducted with the ends of the finite
difference mesh moved away from the protuberance. This lead to a re-
lief of this difficulty. Apparently a larger number of grid points
ahead of the protuberance will have to be used.
The second area of interest needing further study is the conver-
gence rate of the iterative scheme. The present ADI numerical algorithm
is conditionally stable and the convergence rate sensitive to the time
step chosen. For the solutions depicted in Figure 4, convergence was
achieved in about 40 iterations (approximately 10 minutes) and this
computer time can be further reduced by either reducing the number of
grid points (presently set at 5000) or parametrically assessing the
optimum value of et. Through these two efforts it would be anticipated
-	
-4 , _.
-10-
that the computer effort could be reduced to the 2 or 3 minute level.
Efforts in this direction will be made during the second half of this
study.
Another point of interest involves the sensitivity of the pre-
sent'sol,ation method to the initialization procedure. At present
the solution method is found to show some influence of initial con-
ditions (time - 0) on the final steady state values, mainly in the
separated r -ion and the source of this anomaly is presently not known.
This problem could be due to either a programming error (a possibility
now being investigated) or a conceptual error. It is possible that
the influence of the confined finite difference mesh disdusse.l above
is causing this problem, but it is not yet clear how this would occur.
Another possibility is that the consistency error of the ADI scheme
(Ref. 10) is producing this difference and this possiblity will also
be considered.
Following this an eddy viscos-ty model which is more appropriate
for separated boundary layers will be incorporated into the algorithm
before comparisons with experimental data for single and multiple
protuberances will be made.
REFERENCES
10	 Polak, A., Werle, M.J., Vatsa, V.N. and Bertke, S.D., "Supersonic
Laminar Boundary-Layer Flow Past a Wavy-Wall with Multiple Separa-
tion Regions," Rept. AFL 74-12-15, Department of Aerospace Engineer-
ing, University of Cincinnati, Cincinnati, Ohio, February 1975.
To appear in Journal of Spacecraft and Rockets, 1976.
2.	 Harris, J.E., "Numerical Solution of the Eq+lations for Compressible
'Laminar, Transitional on Turbulent Boundary Layers and Comparisons
with Experimental Data," NASA TR-R-368, August 1961.
-11-
3. Bertke, S.D., Werle, M.J. and Polak, A., "Finite Difference
Solutions to the Interacting Supersonic Turbulent Boundary
Layer Equations, Including Separation Effects," Rept. AFL 74-4 -91
Department of Aerospace Engineering, University of Cincinnati,
Ci.ncinnati, , Ohio, April 1974.
4. Brandon, H.J., Masek, R.V. and Dunavant, J.C., "Aerodynamic
Heating to Corrugation Stiffened Structures in Thick Turbulent
Boundary Layers," AIAA Paper No, 75-190, presented at the AIAA
13th Aerospace Sciences Meeting, Pasadena, California, January
20-22, 1975.
5. Shang, J.S. and Hankey, W.L., Jr., "Numerical Solution of the
Navier-Stokes Equations for Supersonic Turbulent Flow Over a
Compression. Corner," AIAA Paper No. 75-3, presented at the AIAA
13th Aerospace Sciences Meeting, Pasadena, California, January
20-22, 1975.
6. Alber, I.E., "Similar Solutions for a Family of Separated Turbu-
lent Boundary Layers," AIAA Paper No. 71-203, presented at the
AIAA 9th Aerospace Sciences Meeting, New York, New York, January
25-27, 1971.
7. Davis, R.T., Personal Communication, September, 1975.
8. Blottner, F.G., "Computational Techniques for Boundary Layers,"
AGARA Lecture Series No. 73, Com utational Methods for Inviscid
and Viscous Two-and-Three Dimensional F ow Fie s, Fe ruary 	 5,
pp. 3-1 to - .
9. Dunavant, J.C., Private Communication, August, 1975.
10. Werle, M.J. and Vatsa, V.N., "New Method for Supersonic Boundary
Layer Separations," AIAA Journal, Vol. 12, No. 11, November 1974,
pp. 1491-1497.
t	 i	 f
}	 t f
	
t	
i	
t t ,,r
	
,j
t	 K	 , --	 it }	 }.1	 i_	 1	 i	 _	 t
!	 I	 '	 },!	 j	 21	 i.
	
' i — ....
	
I .
	
' '^ J	 —^	 t	 t	
-	 i	 }
im
is
jk
lk
 rl,
EL
I J11—
J1.1
	
--.	 ;	 -)	 1	 I	 N	
U
41 t	 ,.-,.may;-_.-._.	 ^ _	 ,^..	 _	 ^^.r^'.	 ^.._	 .-}	 1	 -_ 1 . - _^	 .^ t	 •^
4	 ,.41M
t	 t 1	
t	
si 'j, 
LIO
1 t	
^_L
_	 t
i	 !	 t	 '	 i	 J	 J	 --!	 i	 v
/1
if
}
O_
:
,
-
-^
	 -1
:	 I	 I
c
r
>f
tit
I	 t
--	 _—
---
'	 ^ J
,
' l
1
_L
!- -t
l
-72
_	
_	 _	
j -1
i
1
,1 f 1 -rt_1
I _ _
t
,	 I
IL
,
_
s
Q.L
s
^h
t
_6
u
Lt
iz
v
M
Q
1j	 I
_
-
___
-
--	
-1t-,-	 X---^
	
}---	 t	 ?
_ ^ - 
i LIL
_^	 --1-
_
--.'	 -	
^	 Ill
1	
^	 ,
I
-	 -	 -
I_
i	 1	 I


Numerical study of supersonic turbulent flow over small protuberances

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server