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

DIRECT NUMERICAL SIMULATIONS OF SHOCKLET-CONTAININGTURBULENT CHANNEL COUNTER-FLOWSArash Hamzehloo
Turbulence Simulation Group
Department of Aeronautics
Imperial College London
London SW7 2AZ, UK
a.hamzehloo@imperial.ac.uk

David J. Lusher
Sylvain Laizet
Aerodynamics and Flight Mechanics Group Turbulence Simulation Group
University of Southampton
Department of Aeronautics
Southampton SO16 7QF, UK
Imperial College London
d.lusher@soton.ac.uk
London SW7 2AZ, UK
s.laizet@imperial.ac.uk

Neil D. Sandham
Aerodynamics and Flight Mechanics Group
University of Southampton
Southampton SO16 7QF, UK
n.sandham@soton.ac.uk

ABSTRACT

amenable to Direct Numerical Simulation (DNS), and demonstrated that it could overcome the above-mentioned barriers
associated with free shear layers and Poiseuille/Couette type
flows. Specifically, it retains a statistically stationary onedimensional solution, in common with conventional channel
flows, but contains an inflectional mean flow, representativeof free shear layers. The counter-flow channel has periodic
streamwise and spanwise boundaries and isothermal no-slip
walls and is driven by a mean pressure gradient introduced bya hyperbolic tangent forcing term. It was shown that when
the peak local mean Mach number reached ∼ 0.55, a turbulent Mach number of ∼ 0.6 could be obtained, indicating thatsuch flow configuration could potentially be useful for studying compressibility effects on turbulence.
Our earlier study (Hamzehloo et al., 2021) was limited toa maximum Mach number of M = 0.4 (defined based on a reference velocity deduced from the forcing as discussed in thenext section) since higher values necessitated using a shockcapturing method. Here, by applying an appropriate numerical treatment for flow discontinuities, we extend our previous
DNS study to a Mach number of M = 0.7. This help us investigate, for the first time in such acounter-flow configuration,high compressibility effects, and the formation and evolutionof shocklets, as well as their interactions with turbulence.

Counter-flow or counter-current configurations can maintain high turbulence intensities and exhibit a significant levelof mixing. We have previously introduced a wall-bounded
counter-flow turbulent channel configuration (Physical Review Fluids, 6(9), p.094603.) as an efficient framework to
study compressibility effects on turbulence. Here, we extend
our previous direct numerical simulation study to a relatively
higher Mach number (M = 0.7) to investigate strong compressibility effects (also by reducing the Prandtl number from
Pr = 0.7 to 0.2), and the formation and evolution of unsteadyshocklet structures. It is found that the configuration is able toproduce highly turbulent flows with embedded shocklets andsignificant asymmetry in probability density functions of dilatation. A peak turbulent Mach number close to unity is obtained, for which the contribution of the dilatational dissipation
to total dissipation is nevertheless found to be limited to 6%.

INTRODUCTION
Traditionally, free shear layers (such as jets, wakes andmixing layers) and Poiseuille/Couette type flows (such as
channel flows) have been utilized to study compressible turbulence (Freund et al., 2000). Spatially-developing mixing layer
simulations are computationally expensive and sensitive to farfield and inflow/outflow boundary conditions, but the basic
compressibility effects are captured in temporal simulations(Vreman et al., 1996). A drawback of such simulations is thatthe shear layer thickens as time progresses and large structures
swiftly fill the domain. On the other hand, Poiseuille/Couette
type flows are relatively efficient to compute and can achieve
high Reynolds numbers for relatively low computational costs.
However, they are limited in terms of the turbulent Mach number that can be achieved.
Counter-flows are highly efficient mixers due to thehigh turbulence intensities that can be maintained (Humphrey& Li, 1981; Strykowski & Wilcoxon, 1993; Forliti et al.,2005). Previously, we introduced (Hamzehloo et al., 2021)
a wall-bounded counter-flow turbulent channel configuration,

METHODOLOGY
Computational Approach
The dimensionless governing equations of a compressible
Newtonian fluid flow that conserve mass, momentum and energy are solved (Hamzehloo et al., 2021). The latter two equations containing the shear-forcing term are given as:

∂
∂ ρui
+
(ρui u j + pδi j − τi j ) + c j δi j = 0,
∂t
∂xj
∂ ρE
∂
+
(ρEu j + u j p + q j − ui τi j ) + c j u j = 0,
∂t
∂xj

1

(1)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Forcing
term
prole

y
x
y=0

Mean
velocity
prolePeriodic
boundary

Isothermal no-slip wall

Isothermal no-slip wall

Figure 1.

2D schematic of the 3D counter-flow channel configuration.
Table 1.

DNS counter-flow channel cases.

Case

M

Re

Pr

WENO Filter

∆t

1

0.1

400

0.7

No

2 × 10−4

2

0.4

400

0.7

No

3

0.4

400

0.7

Yes

4

0.7

400

0.7

Yes

5

0.7

400

0.2

Yes

5 × 10−5

5 × 10−5

5 × 10−5

5 × 10−5

{u}b

{u} p

⟨ac ⟩b⟨ac ⟩ p⟨T ⟩ p

⟨M⟩ p

Mt p

1.409

2.159

10.29110.479

1.098

0.207

0.213

1.375

2.047

3.669

3.849

2.374

0.556

0.595

1.380

2.062

3.671

3.852

2.376

0.560

0.594

1.358

2.005

2.919

3.132

4.823

0.697

0.758

1.535

2.238

2.382

2.602

3.332

0.973

0.981

Problem Specifications

where ρ represents the density, ui (i = 1, 2, 3) denotes the velocity component (u, v and w, respectively) in the ith direction (x, y and z, respectively), E is the total energy, and p
and δ denote the pressure and the Kronecker delta, respectively. The forcing term c j drives the flow, with a value of
c1 = −c0 tanh(ay) in the x direction and zero in other directions. The maximum value is set as c0 = 1. Since −1H <y < 1H, where H = 1 is the channel half height here, equivalent driving forces are applied in opposite directions to the upper and lower halves of the domain which consequently resultin the formation of a shear-forcing or counter-flow condition.
The coefficient a in the forcing term is positive with a valueof a = 100. The viscous stress tensor (τi j ) and the heat flux∂ujµ ∂ ui2 ∂ uk(q j ) are defined as τi j =+−δi j andRe ∂ x j∂ xi
3 ∂ xk−µ∂T
qj =, respectively. Here, µ denotes the(γ − 1)M 2 PrRe ∂ x j
dynamic viscosity, T is the temperature, γ is the ratio of specific heats with a value of γ = 1.4 here and Pr is the Prandtl
number. Re and M denote the Reynolds and Mach numbers
based pon a reference velocity deduced from the forcing asu∗ref = c0 H/⟨ρ⟩b , where ⟨ρ⟩b is the bulk-averaged density,
together with the channel half height and the wall temperature and viscosity. Angle brackets ⟨⟩ denote averages over the
homogeneous spatial directions (x and z) and time. The additional subscript b here denotes an additional average over y.
The dynamic viscosity is calculated as µ = T 0.7 . The temperature is calculated as T = pγM 2 /ρ. Here the pressure of an
ideal Newtonian fluid is obtained using an equation of state as1
p = (γ − 1)(ρE − ρui ui ).
2

As shown schematically in figure 1, the streamwise andspanwise boundaries of the counter-flow channel configuration are periodic, while isothermal (Tw = 1.0) no-slip walls
are assigned to the boundaries in the normal direction (y). In
order to accurately resolve the near wall region, the grid isstretched in the y direction. We have previously identified an
optimum domain size of 12H × 2H × 6H with a grid resolution of 240 × 151 × 200 (Hamzehloo et al., 2021). This domain size is used with a Reynolds number of Re = 400. Mach
number values of M = 0.1, 0.4 and 0.7 are examined. The
Prandtl number has a value of mainly Pr = 0.7. However, acase with Pr = 0.2 (M = 0.7) is also studied. This reductionin Pr increases the the wall heat transfer, reducing the bulk
temperature and sound speed in the channel, hence increasing
the compressibility effect. The previously-mentioned WENO
filtering is necessary to obtain a stable solution for the caseswith M = 0.7. However, for the case with M = 0.4, results of
two simulations without and with the filtering are provided tomake direct comparisons and examine the filtering effect. A
list of test cases studied here is presented in table 1. It should
be noted that, here, subscript p denotes the peak values of flow
quantities.
Here, the single prime ′ denotes the turbulent fluctuation which for an arbitrary flow quantity (φ ) is defined as
φ ′ = φ − ⟨φ ⟩. Moreover, for the higher Mach number case,
the Favre average is defined as {φ } = ⟨ρφ ⟩/⟨ρ⟩ and the double prime ′′ denotes the turbulent fluctuation with respect to
the Favre average defined as φ ′′ = φ − {φ }. For the Reynolds
stresses, the Favre average is related to the Reynolds averageas ⟨ρ⟩{u′′i u′′j } = ⟨ρui u j ⟩ − ⟨ρ⟩⟨ui ⟩⟨u j ⟩. Also, the mean Machpnumber is defined as ⟨M⟩ = ⟨u⟩2 + ⟨v⟩2 + ⟨w⟩2 /⟨ac ⟩, where
ac is the local speed ofp
sound, while the turbulent Mach number is defined as Mt = ⟨u′ u′ ⟩ + ⟨v′ v′ ⟩ + ⟨w′ w′ ⟩/⟨ac ⟩.

A fourth order finite-difference central scheme is used todiscretise the equations recasted in split skew-symmetric formulations to improve stability. A Lax-Friedrichs Weighted
Essentially Non-Oscillatory (WENO) filter (Yee & Sjögreen,2018) is applied to the flow field after the completion of eachtime step to improve the solution stability in the presence ofshocklets. The density field is corrected if necessary, after applying the filtering to maintain the conservation of mass. A
low-storage three-stage explicit Runge-Kutta scheme is used
to advance the solution in time. The simulations are performed
with the flow solver OpenSBLI (Lusher et al., 2021).

RESULTS AND DISCUSSION
Mean Flow And Turbulence Statistics
Figure 2 shows a direct comparison between the counterflows studied here based on various mean flow quantities, including the streamwise velocity {u}, density ⟨ρ⟩, temperature⟨T ⟩, local speed of sound ⟨ac ⟩ and Mach number ⟨M⟩, and also212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022∂u ∂v ∂w++
, is lower than a negative thresh∂x ∂y ∂z
old i.e. θ < −ζ (Samtaney et al., 2001; Wang et al., 2017).
A value of ζ = 3θ ′ , where θ ′ denotes the dilatation fluctuation (roots mean square of the dilatation magnitude) defined

the turbulent Mach number Mt . Additionally, figure 3 provides
the Favre Reynolds stresses of the counter-flows. The flow
significantly heats up when the Mach number increases, forinstance, as also provided in table 1, the peak mean temperature is ∼ 2.0 and ∼ 4.4 times higher for the cases with M = 0.4
and M = 0.7 compared to the case with M = 0.1, respectively.
However, reducing the Prandtl number from Pr = 0.7 to 0.2 for
M = 0.7, reduces the peak mean temperature by ∼ 31%. As
shown in figure 2, the mean velocity does not change significantly when changing the Mach number since the latter is governed by the local speed of sound. Mean and turbulent Mach
numbers exhibit significant surges when M increases. Specifcially, the peak turbulent Mach number increases by ∼ 179%and ∼ 256% when the Mach number increases from M = 0.1
to M = 0.4 and 0.7, respectively. Prandtl number reduction
while keeping the Mach number constant at with M = 0.7, further increases the turbulent Mach number by ∼ 29.5% helping
it reach a value of near unity at the channel centreline.
With respect to the normal stresses, as shown in figure 3,
by increasing the compressibility through increasing the Mach
number and/or reducing the Prandtl number, the peak streamwise stress increases and the peaks of the other two normal
stresses reduce. For instance, the peak ⟨ρ⟩{u′′ u′′ } increases
by around ∼ 9.5% when the Mach number increases from
M = 0.4 to 0.7. Reducing the Prandtl number while keeping the Mach number constant at M = 0.7 shows a noticeable
effect on the normal stresses particularly the normal streamwise stress. Specifically, reducing the Prandtl number from
Pr = 0.7 to 0.2 increases the the peak ⟨ρ⟩{u′′ u′′ } stress by∼ 20%. However, the Favre shear stress ⟨ρ⟩{u′′ v′′ } continue
varying linearly with y as the compressibility increases. The
shear stress is required to exhibit such trend by the imposed
forcing term as discussed previously (Hamzehloo et al., 2021).
′
q Figure 4 shows the vorticity fluctuation (ω =defined as θ =

∂u 2
∂v
∂w 2
) + ( )2 + () ⟩, has been used in com∂x∂y∂zpressible decaying turbulence problems to detect shocklets asregions with strong compression rates (in turbulent Mach number values in the range of 0.5 ≤ Mt ≤ 1.0) (Samtaney et al.,
2001; Wang et al., 2017). Also, ζ = θ ′ was used by (Samtaney et al., 2001) to visualise shocklets. Table 2 provides thebulk-averaged values of the dilatation fluctuation (θb′ ) for thecounter-flows studied here.
Following the literature, in the present study, in order todetect and visualise the shocklets, iso-surfaces of the dilatation with various threshold values (iso-values) of ζ = 2θb′ , 3θb′and 3θb′ |M=0.7, Pr=0.7 = 5.655 are used as shown in figure 5
for counter-flows with M ≥ 0.4. The latter iso-value is usedbased on the bulk-averaged dilatation fluctuation of the casewith M = 0.7 and Pr = 0.2 as a fixed reference with the aimof making a direct and relatively fairer comparison between
the counter-flows studied. Based on what is seen from figure 5, and also observations reported in (Wang et al., 2017)a threshold value of ζ = 3θb′ can potentially be used to highlight the differences in high-compression regions of the flow
(i.e. shocklet) as the compressibility increases in the counterflow configuration. However, from subfigure (c) of figure 5,
a fixed threshold may help better understand the flow differences. Specifically, as the compressibility increases, the number and size of the dilatation iso-surfaces increase significantly.
It can be concluded that with M = 0.7 the counter-flow canproduce a shocklet-containing flow with significant numbersof heterogeneously-distributed highly irregular shocklets.
as θ ′ =⟨ωx2 + ωy2 + ωz2 ⟩) and the components of the vorticity
turbulent fluctuations for different Mach and Prandtl numbers.
With Pr = 0.7, by increasing the Mach number, values ofthe vorticity fluctuations (total and all components) reduce
noticeably. This trend is comparable to the relationshipbetween the spanwise and wall-normal Reynolds stresses andthe Mach number as seen in figure 2. The vorticity reduction
is a sign of the three-dimensional structures weakening andturbulence stabilisation. Reducing the Prandtl number from
Pr = 0.7 to 0.2 slightly increases the vorticity fluctuations.
This is attributed to the significant increase in the streamwise
Reynolds stress fluctuation as seen in figure 3 and also an
increase in the mean velocity {u} as shown in figure 2 andtable 1. The peak mean streamwise velocity increases byaround ∼ 13.0% as the Prandtl number reduces.
As shown in figure 2, the counter-flow with M = 0.4
forms a mildly compressible flow with Mt p ≈ 0.6, hence theformation of discontinuities in the form of shocklets is limited.
Figures 2 and 3 show that the cases with M = 0.4 without andwith the WENO filtering exhibit almost identical trends andtable 1 shows less than ∼ 1% changes in their key flow quantities. This suggests that the filter-based shock capturing method
is performing well and the grid resolution is fine enough.

⟨(Table 2.

Bulk-averaged dilatation fluctuation.

Case

M

Pr

θb′

1

0.1

0.7

0.128

3

0.4

0.7

1.374

4

0.7

0.7

1.885

5

0.7

0.2

2.808

Shocklet Quantification
The total viscous dissipation can be divided into thesolenoidal and dilatational components as ε T = ε S + ε D(Sarkar et al., 1991) where, the solenoidal (ε S ) and dilatational(compressible) (ε D ) dissipations are defined as+#
εS =
" *






1
dw dv 2du dw 2
dv du 2
µ
−
+
−
+
−
and
Re
dz dz
dz dx
dx dy
" *
+2 #
4
du dv dw
εD =
µ
+
+
.
3Re
dx dy dz
Energy dissipation through Mach number-induced changes on
turbulent flow structures (i.e. shocklets) can be linked to thedilatational part of the dissipation (Sarkar et al., 1991). Figure
6 shows the ratio of the dilatational dissipation to the total dissipation for the counter-flows studied here. The contributionof the dilatational part of the total viscous dissipation becomes

Shocklet StructureWith the values of the mean and turbulent Mach numbers
seen in figure 2, regions with instantaneous Mach numbers beyond unity are expected to form in the cases with M = 0.7. For
such transonic values, the formation of shocklets is possible.
Shocklets can be associated with regions where the dilatation,3

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

Figure 2. Profiles of the mean velocity, density, temperature, speed of sound and Mach number, and profile of the turbulent Mach
number (Mt ).

Figure 3.

Profiles of the Favre Reynolds stresses.

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

Figure 4.

Profiles of the vorticity fluctuation and the components of the vorticity turbulent fluctuations.

M=0.4
Pr=0.7

Pr=0.2(a)(c)

Figure 5. Iso-surfaces of the dilatation with iso-values of (a): θiso = 2θb′ , (b): θiso = 3θb′ , and (c): θiso = 3θb′ |M=0.7, Pr=0.7 = 5.655.
Red and blue colours show the flow directions in the positive (left to right) and negative streamwise directions, respectively.

plotted for the cases with M = 0.7 in figure 7. The red curves
show the instantaneous PDFs and the thick black curves show
their average in time. All PDF profiles are skewed towards
the negative dilatation values, a trend expected due to the existence of nonlinear compression waves and shocklets in suchcompressible flows. Instantaneous PDF trends of figure 7 confirms that the dilatation, hence the number and strength of theshocklets, exhibits a significant fluctuating transient behaviour.
However, as expected, the counter-flow with Pr = 0.2 exhibitsa more negatively skewed profile of the time-averaged dilatation PDF consistent with its higher compressibility level.

significantly more important as the compressibility increases
by having a higher Mach number and/or a lower Prandtl number. The peak of the ratio of the dilatational dissipation to
the total dissipation occurs at around |y| ≈ 0.75 for all cases
where the mean streamwise velocity and the mean Mach number also exhibit peak values as shown in figure 2. It is clearthat ε S is directly related to the vorticity which reduces as theMach number increases. In fact, the solenoidal dissipation reduces slightly (not shown here) as the compressibility of thecounter-flow increases. A higher compressibilty results in theformation of more shocklets (as shown in figure 5), and hencea relatively higher dilatational dissipation.
To provide a more quantitative analysis of the distributionof the shocklets, the Probability Density Function (PDF) of thedilation over the entire domain for time intervals of t = 10 are

5

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

M = 0.1, P r = 0.7
M = 0.4, P r = 0.7

less).
Overall, the counter-flow configuration was found to beable to produce highly turbulent flows with embedded shocklets for a relatively modest Mach number. Therefore, the configuration provides a useful framework to study some of thefundamental physics associated with shock-turbulence interactions and could be considerably beneficial to the development of compressible subgrid scale turbulence models as wellas improved representations of compressibility in Reynoldsaveraged Navier-Stokes models.

M = 0.7, P r = 0.7
M = 0.7, P r = 0.2

hǫD i/hǫT i

0.06
0.04
0.02
0.00
−1.0

−0.5

0.0

0.5

1.0

y

Figure 6. Profiles of the ratio of the dilatational dissipation
over the total dissipation.

ACKNOWLEDGMENTS
Dr Arash Hamzehloo was funded by the UK Turbulence
Consortium (EPSRC grant EP/R029326/1). Dr David J Lusher
was funded by EPSRC grant EP/L015382/1. The authors acknowledge the use of the Cambridge Tier-2 system operated by
the University of Cambridge Research Computing Service under an EPSRC Tier-2 capital grant (EP/P020259/1). The flow
solver OpenSBLI is available at https://opensbli.github.io.

REFERENCES
Forliti, D.J., Tang, B.A. & Strykowski, P.J. 2005 An experimental investigation of planar countercurrent turbulent
shear layers. Journal of Fluid Mechanics 530, 241.
Freund, J.B., Lele, S.K. & Moin, P. 2000 Compressibility effects in a turbulent annular mixing layer. Part 1. Turbulence
and growth rate. Journal of Fluid Mechanics 421, 229–267.
Hamzehloo, A., Lusher, D.J., Laizet, S. & Sandham, N.D.
2021 Direct numerical simulation of compressible turbulence in a counter-flow channel configuration. Physical Review Fluids 6 (9), 094603.
Humphrey, J.A.C. & Li, S. 1981 Tilting, stretching, pairingand collapse of vortex structures in confined counter-current
flow. Journal of Fluids Engineering 103 (3).
Lusher, D.J., Jammy, S.P. & Sandham, N.D. 2021 OpenSBLI: Automated code-generation for heterogeneous computing architectures applied to compressible fluid dynamics
on structured grids. Computer Physics Communications p.
108063.
Samtaney, R., Pullin, D.I. & Kosović, B. 2001 Direct numerical simulation of decaying compressible turbulence andshocklet statistics. Physics of Fluids 13 (5), 1415–1430.
Sarkar, S., Erlebacher, G., Hussaini, M.Y. & Kreiss, H.O. 1991
The analysis and modelling of dilatational terms in compressible turbulence. Journal of Fluid Mechanics 227, 473–
493.
Strykowski, P.J. & Wilcoxon, R.K. 1993 Mixing enhancement
due to global oscillations in jets with annular counterflow.
AIAA Journal 31 (3), 564–570.
Vreman, A.W., Sandham, N.D. & Luo, K.H. 1996 Compressible mixing layer growth rate and turbulence characteristics.
Journal of Fluid Mechanics 320, 235–258.
Wang, J., Gotoh, T. & Watanabe, T. 2017 Shocklet statistics incompressible isotropic turbulence. Physical Review Fluids
2 (2), 023401.
Yee, H.C. & Sjögreen, B. 2018 Recent developments in accuracy and stability improvement of nonlinear filter methods
for dns and les of compressible flows. Computers Fluids
169, 331–348.

Figure 7. Probability distributions of the dilatation over theentire domain for every t = 10 intervals up to t = 1000.

CONCLUSIONS
Direct numerical simulations of shocklet-containing turbulent flows were conducted using a new counter-flow channel configuration introduced previously by the current authors(Hamzehloo et al., 2021). The simulations were performed forMach number values of M = 0.1, 0.4 and 0.7 with a Reynolds
number of Re = 400. A case with a Prandtl number of Pr = 0.2(reduced from the default value of Pr = 0.7) and M = 0.7 was
also studied in an attempt to boost the compressibility.
It was found that a Mach number as low as M = 0.7 could
produce fluctuating and mean Mach numbers close to 0.7 overa considerable length of the width of the counter-flow channel. Such values produced instantaneous supersonic velocities
that formed relatively strong transient shocklets. Additionally,
fluctuating and mean Mach numbers above 0.95 were achievedby reducing the Prandtl number from Pr = 0.7 to 0.2. Despite
the presence of shocklets the contribution of dilatational dissipation to the total dissipation was found to be small (6% or6