12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

DIRECT NUMERICAL SIMULATIONS OF TURBULENT AIR INTAKE IN A
HYPERSONIC FLOW

Youcheng Xi
School of Aerospace Engineering
Tsinghua University
Beijing, China, 100084
xiyc@mail.tsinghua.edu.cn

Song Fu
School of Aerospace Engineering
Tsinghua University
Beijing, China, 100084
fs-dem@tsinghua.edu.cn

ABSTRACT

surface skin friction. Details of the numerical test case are described in Section 2. Results and some discussions are given
in Section 3. The final conclusions are presented in Section 4.

In this study, the numerical simulations of turbulent flowin an air intake are performed to understand the characteristics
of skin friction in hypersonic inner flows. A direct numerical simulation of a boundary layer is performed first to obtainturbulent information at the inlet of an intake, with the full process of transition and the growth of the boundary layer. Thenthe simulation of an air intake is undertaken. Integral analyses
based on Cauchy’s formula have been performed to reveal a
detailed ingredients of skin-friction. It is found that the turbulence induced drag is the most dominant, in this case. The
amplifications of turbulence is found after the interaction of
the shockwaves and turbulent boundary layer, which lead to
obvious increase of skin-friction. The results indicate that to
reduced the drag in a hypersonic air intake, we need to impair
turbulence.

1

2

Methodology

The governing equations for all simulations in this work
are the dimensionless compressible Navier–Stokes equations
for a Newtonian fluid, which can be written as:v∂ Q ∂ Fj ∂ Fj
+
+
= 0,
∂t
∂xj
∂xj

(1)

Q = [ρ, ρu1 , ρu2 , ρu3 , Et ]T ,

(2)




0
ρu j
 ρu1 u j + pδ1 j 


τ1 j

 v 



.
ρu
u
+
pδ
τ
,
F
=
Fj = 
2j  j2j
 2 j ρu3 u j + pδ3 j τ3 j(Et + p) u jτ jk uk − q j(3)

Introduction

The hypersonic scramjet is considered to be one of thebest propulsive choices for hypersonic flight(Urzay, 2018).
Therefore, it is crucial to understand the aerodynamic behaviour for the device. Most air-breathing flights are designed as wave-rider, and the airflow is compressed dynamically through an intake system integrated with the fore-body.
The flow is asked to become turbulent ahead of the combustion
chamber in the air intake for the necessities of complete combustion. However, based on experimental studies, although thesurface area of a scramjet combustion chamber only forms asmall part, the high shear along its walls can account for morethan half of the total engine skin-friction drag(Paull et al.,
1995; Robinson et al., 2006). A reduction of the combustion
chamber skin-friction drag therefore promises relatively large
margins of performance improvement.
Therefore, it becomes a natural question before applying
a method of drag reduction, what features does the flow havein the hypersonic air intake? In this research, direct numerical simulations, which provide a large amount of accurate datathat can be used to analyse, are used to reveal the precise turbulent features in the hypersonic air intakes, focusing on theThroughout this work the coordinates xi , (i = 1, 2, 3) are referred to as x, y, z for the stream-wise, bottom wall-normal andspanwise directions, respectively, with corresponding velocity
components u, v, w.
The total energy Et and the viscous stress τi j are given as,
respectively,



uk uk
T
+
,
2
γ(γ − 1)M∞2
µ∂ ui ∂ u j 2 ∂ uk
τi j =+− δi j.
Re∞ ∂ x j ∂ xi 3 ∂ xk
Et = ρ

1(4)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
The computational domain, mesh number, and mesh resolution in wall unit are shown in Table 2 and 3. For thetransitional boundary layer, 4482 grid points are used along
the stream-wise direction. 140 grid points are used along thewall-normal direction and grid clustering is used to ensure atleast 50 grid points are inside the boundary layer. The spanwise direction z is discretized by 256 grid points. For the
turbulent air intake, 1551, 401 and 256 grid points are usedalong stream-wise, wall-normal and span-wise direction, respectively. Based on inner scale, the stream-wise grid spaceis δ x+ = 6.35 and 6.4. The wall-normal grid space for the
first grid to the wall is unchanged in the whole domain, andδ y+ = 0.47 and 0.5. The wall-normal mesh is gradually coarsened away from wall surfaces. Due to the development of turbulence along the stream-wise direction, spanwise grid spaceis δ z+ = 3.72 in plate boundary layer, and is 4.5 in the intake.

The pressure p and heat flux qi are obtained from:

p=

ρT
,
γM∞2

qi = −µ∂T.
(γ − 1)M∞2 Re∞ Pr ∂ xi

(5)

The viscosity is calculated using the Sutherland lawµ = T 3/2 (T∞ +C)/(T T∞ +C),

(6)

with C = 110.4K. The free-stream Reynolds number Re∞ ,Mach number Ma∞ and Prandtl number Pr are defined as

Re∞ =ρ∞∗ U∞∗ l0∗,µ∞∗Ma∞ = pU∞∗,
γR∗g T∞∗Pr = 0.72,(7)

Table 2. Summary of mesh resolution for numerical simulation of transitional boundary layer

where ρ∞∗ , U∞∗ , T∞∗ and µ∞∗ stand for the freestream density, velocity, temperature and viscosity, respectively. R∗g = 287J/(K ·Kg) represents the gas constant and γ stands for the ratio of
specific heat. The length scale l0∗ is chosen as 1(mm) in thisresearch.
The parallel computational code OPENCFD, developed
by Li et al. (2008), has been used for the direct numerical simulation(DNS). The compressible Navier–Stokes equations are
solved numerically using high-order finite-difference methods.
Convection terms are discretized with Steger–Warming splitting and a seventh-order weighted essentially non-oscillatory
scheme. Viscous terms are discretized with an eighth-order
centred finite-difference scheme. A third-order total variation
diminishing-type Runge–Kutta method is used for the time.
Validations and verifications of the code can be found in a series of papers (Li et al. (2008); Liang et al. (2010); Li et al.
(2010)).
The simulation processes are divided into three steps. At
first, the simulation is conducted over a two-dimensional blunt
plate with the high-order shock-fitting method(Xi et al., 2021).
During the two-dimensional analysis, the inlet profile of theboundary layer is obtained. Then, a transitional flat boundary
layer is calculated with the computational domain inside theshock layer, with the inlet profile of the boundary layer. Wallblowing and suction perturbations are used to trigger the transition as soon as possible, with substantial amplitude. Whenthe turbulent boundary layer is fully developed, the unsteadyflow field at the specific stream-wise station (x = 1400) is
recorded and used as inlet information for the Intake simulation. Finally, with the recorded turbulent information, a simulation of air intake is performed for investigation.

Table 1.

Domain Size

x(m)

y(m)

z(m)

1.33

0.036

0.03

Reτ

C f × 10−3

δ x+

δ y+

δ z+

249 − 582

1.3-1.1

6.35

0.47

3.72

Reθ

Re∞

Nx

Ny

Nz

1424 − 5954

7123

4482

140

256

Table 3. Summary of mesh resolution for the numerical simulation of intake

Domain Size

x(m)

y(m)

z(m)

0.3

0.024

0.03

Rein
τ

−3
Cin
f × 10

δ x+

δ y+

δ z+

552

1.1

6.4

0.5

4.5

Reθ

Re∞

Nx

Ny

Nz

4986 − 9715

7.12 × 106

1551

401

256

Basic Parameters for simulations

Ma∞

T∞ (K)

Re∞

γ

Pr

3

4.8

1314

7123

1.4

0.72

The snapshots of vortices by sio-surface of the λ2 criterion with λ2 = −0.05 are shown in figure 1. The regionx ∈ [1000 − 1100] are enlarged on the bottom right. On the top
left, the Reynolds- and spanwise-averaged van Driest velocity
UV+D (Van Driest, 1951) profiles at x = 1400 is compared withtheoretical results,

Basic parameters for the present simulations are shown in
Table 1. For both simulations, adiabatic surface boundary condition is assumed and periodic boundary conditions are usedalong the spanwise direction. Nonreflecting boundary condition is implemented on the outflow boundary along the streamwise direction.

Result and Discussion

y+ =

2

ρw uτ y +uV D
, uV D =, uτ =µwuτr

τw, uV D =
ρw

Z ur
0

ρ
du,
ρw
(8)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
relationwhere τw and ρw are the wall shear stress and surface density.
As shown in the figure, the result agrees well with all the theoretical values in both the inner layer and the log layer, whichindicate that the turbulent boundary layer is fully developed.
The snapshots for air intake of vortices by iso-surfaceof the λ2 -criterion with λ2 = −0.03 are shown in figure 2.
The Reynolds- and spanwise-averaged skin-friction coefficientC f of the bottom is exhibited at the top left, where C f =2τw / ρ∞∗ U∞∗2 . From the figure, the increases of the coefficient are mainly due to the interaction between the shock
wave and the turbulent boundary layer. The leading shock
wave lead to the obvious increase of the skin friction at aroundx ∈ [60, 80]. Then the coefficient exhibit some decrease aheadof the reflected shock wave at around x ∈ [210, 230]. The shock
wave plays as an amplification which apparently promote theintensity of turbulence through the interaction, which is directly reflected in the skin friction.

y−C f b −Cdf

=
2
1 ∂u
1 ∂u
−
=
µ
µ
Re ∂ y y=yb Re ∂ y y=0
00 u00Z yb∂ ρ ug
1
p
τ
∂
1
∂
xx
−
 dy
+
−
∂x
γMe2 ∂ x Re ∂ x
0
−

Z yb

(13)

yb
∂ ũ
∂ ũ
dy −
ρ ṽ dy.
∂x
∂y
0

Z

ρ ũ
0

At the same time, twofold repeated integration of theequation (9) lead toZ yb Z ξ1

∂ ũ
∂ ũ
+ ρ ṽ dydξ1 =
∂
x
∂y
0
0


Z yb Z ξ1
1 ∂ τ xx ∂ τ xy
1 ∂p
+
+
−
∂y
γMa2∞ ∂ x Re ∂ x
0
0




00
00
00
00
g
g
∂ ρu u∂ ρu v−−dydξ1 .
∂x
∂y
ρ ũ

3.1

Integral analysis

In order to understand the mechanism of the skin friction, a direct connection between the skin friction and flow
field is needed. Integral identities are valuable tools for theinvestigation of wall-bounded flows. By a wall-normal integration, mean-flow quantities varying in the wall-normal direction through the shear layer are reduced to scalar quantities, allowing for a compact characterization of their streamwise evolution. Following the studies of Fukagata et al. (2002)
and Wenzel et al. (2022), we analysis the the behaviour of thex−momentum equation by using integral relations for c f . The
average momental equation along x direction can be written as

Equations (13) and (14), together with the Cauchy’s formula for n repeated integrations
Z yb Z ξ
a



ρ ũ

∂ ũ
∂ ũ
1 ∂p
1 ∂ τ xx ∂ τ xy
+ ρ ṽ
=−
+
+
2
∂x
∂y
∂
x
Re
∂x
∂y
γMa∞




00 u0000 v00∂ ρ ug∂ ρ ug−−∂x∂y

τ xx = µ4 ∂u 2 ∂v−
3 ∂x 3 ∂yτ xy = µ



∂u ∂v
+
∂y ∂x

+ µ0





+ µ0

4 ∂ u0 2 ∂ v0−
3 ∂x 3 ∂y

2τ w
2µ w ∂ u
Cf ==.
Re ∂ yρ∞ u2∞aZ ξn−1af (ξn )dξn · · · dξ2 dξ1 =1(n − 1)!

(9)

Z yb
a

(15)
(yb − y)n−1 f (y)dy,

results in a explicit formula for the skin-friction over the lower
surface

(10)

2
Re yb
|

2
Re yb
|

Z yb

"

2
yb
|

Z yb

−
(11)

2
yb
|

(12)

Z yb 

Cdf =


∂ u0 ∂ v0
+
∂y
∂x

···



where


(14)

0

cMe
f

(y − yb )µ 0

4 ∂ u0 2 ∂ v0−
3 ∂x 3 ∂y

#
dy

{z

}

cVf

0

Z yb
0

0


∂ ū
∂ v̄
µ̄
+ µ̄
dy +
∂y
∂x
{z
}



Z
i
∂ ũ
2 yb h g(yb − y) ṽρ̄dy −
ρ̄ u00 v00 dy +
∂y
yb 0
{z
}
{z
} |
cTf c

cMc
f

"

#
00 u00∂ ũ ∂ ρ̄ ug
1 τxx
1 ∂ p̄
(yb − y) −ũρ̄
−
+
−
dy,
∂x
∂x
Re ∂ x γMe2 ∂ x
{z
}
cDf

(16)
Mc
where CMef stands for the mean-flow shear term; C f , meanconvection term; CTf c , turbulent-convection term; CDf spatialdevelopment term; CVf viscous-stress-fluctuation term.
Figure 3 exhibits the streamwise evolution of skin friction coefficient and each parts based on equation (16) overthe lower surface for this air intake. From the decomposition,
we can find that turbulence induced drag coefficient CVf +CTf cis the most dominant and the mean-flow shear induced termCMef of the flow field do not have any contribution graphically.

Hereinafter, for an arbitrary variable f Reynolds (ensemble)
averages are denoted by an overbar, f . The density-weighted
Favre averages are denoted by a tilde, f˜ = ρ f /ρand fluctuations around the Reynolds and Favre averages are denoted bysingle and double primes, f 0 = f − f and f 00 = f − f˜, respectively.
Integration of the equation (9), from the lower surface y =0 of air intake to upper surface y = yb , leads to the following

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
Li, Xinliang, Fu, Dexun & Ma, Yanwen 2010 Direct numericalsimulation of hypersonic boundary layer transition over a
blunt cone with a small angle of attack. Physics of Fluids
22 (2), 025105.
Liang, Xian, Li, Xinliang, Fu, Dexun & Ma, Yanwen 2010
Effects of wall temperature on boundary layer stability overa blunt cone at mach 7.99. Computers & Fluids 39 (2), 359–
371.
Paull, A., Stalker, R. J. & Mee, D. J. 1995 Experiments onsupersonic combustion ramjet propulsion in a shock tunnel.
Journal of Fluid Mechanics 296, 159–183.
Robinson, M. J., Mee, D. J. & Paull, A. 2006 Scramjet lift,
thrust and pitching-moment characteristics measured in a
shock tunnel. Journal of Propulsion and Power 22 (1), 85–
95.
Urzay, Javier 2018 Supersonic combustion in air-breathing
propulsion systems for hypersonic flight. Annual Review of
Fluid Mechanics 50 (1), 593–627.
Van Driest, E. R. 1951 Turbulent boundary layer in compressible fluids. Journal of the Aeronautical Sciences 18 (3),145–160.
Wenzel, Christoph, Gibis, Tobias & Kloker, Markus 2022
About the influences of compressibility, heat transfer andpressure gradients in compressible turbulent boundary layers. Journal of Fluid Mechanics 930, A1.
Xi, Youcheng, Ren, Jie & Fu, Song 2021 Hypersonicattachment-line instabilities with large sweep mach numbers. Journal of Fluid Mechanics 915, A44.

These indicate that the total drag inside the intake are mainly
induced by turbulence.
In the inner flow of intake, shockwave turbulent boundary layer interaction also play a very important role. In thiscase, the interaction do not lead to separation. Furthermore,
the shockwave turbulent boundary layer interaction will obviously change the gradients for the variables, which is reflectedin the terms of convection CMc
f and spatial development term
D
C f . Only around the position of interactions, do the amplitudes of convection and development term become large.

4

Conclusion

In this study, direct numerical simulations are performed
to understand the mechanism of drag production inside a hypersonic air intake. An integral analysis is performed and it isfound that the turbulent induced drag is the dominant for thisinner flow. Shockwave turbulent boundary layer interaction
acts as an amplifier which obviously enhances the turbulences
and then leads to obviously increase of drag.

REFERENCES
Fukagata, Koji, Iwamoto, Kaoru & Kasagi, Nobuhide 2002
Contribution of reynolds stress distribution to the skin friction in wall-bounded flows. Physics of Fluids 14 (11), L73–
L76.
Li, Xinliang, Fu, Dexun & Ma, Yanwen 2008 Direct numericalsimulation of hypersonic boundary layer transition over a
blunt cone. AIAA Journal 46 (11), 2899–2913.

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

Figure 1. Snapshots of
pvortices by iso-surface of the λ2 -criterion with λ2 = −0.05 for transitional flat boundary layer. The profiles of
u+ = u/uτ and uV+D = ρ/ρw u+ at x = 1400 are shown on the top left.

Figure 2. Snapshots of vortices by iso-surface of the λ2 -criterion with λ2 = −0.03 for turbulent intake. The average skin-friction
coefficient C f of the lower surface is shown on the top left.

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

Figure 3.

Decomposition of the lower surface skin friction C f .

6