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

THE DYNAMICS OF COHERENT STRUCTURES IN A TURBULENT WAKE
PAST A SPHERE AT RE=3700

Fengrui Zhang
School for Engineering of Matter,Transport and Energy
Arizona State University
fzhang82@asu.edu

Yulia T. Peet
School for Engineering of Matter,Transport and Energy
Arizona State University
ypeet@asu.edu

ABSTRACTon the power spectral analysis in the wake of the flow, but without further linking the frequencies to the shape and dynamics
of coherent structures associated with these frequencies.
In the current paper, we identify the characteristic motions associated with specific frequencies in the wake. Subsequently, with the help of the modal analysis, the characteristic
motions are visualized as coherent structures. The results of
the Proper Orthogonal Decomposition (POD) and Dynamic
Mode Decomposition (DMD) for the identification of thesestructures are compared.

The turbulent wake flow past a sphere at ReD = 3700 is investigated via Direct Numerical Simulation (DNS). The characteristic motions in the wake flow, such as vortex shedding
and bubble pumping are identified by the probes placed in thenear wake with a dominating frequency of St = f u∞ /D = 0.22
and 0.004, respectively. The modal analysis is conducted in
the wake area using Proper Orthogonal Decomposition (POD)
and Dynamic Mode Decomposition (DMD). The vortex shedding and bubble pumping motions are also captured by themodal analysis. The results from POD and DMD show comparable patterns of both characteristic motions. For the bubble
pumping motion, the dominating frequency of the corresponding POD mode is St = 0.004, while the DMD mode that is
directly related to the separation bubble has the frequency of
St = 0.009.

NUMERICAL METHODSIn the current study, the incompressible Navier-Stokes
equations are solved via Direct Numerical Simulations using
the open source code Nek5000 (Fischer et al. (2015)). A
high-order spectral element method (SEM) is employed to discretize the spatial terms. In order to avoid spurious pressuremodes, velocity and pressure are approximated with different
degree polynomials by a PN /PN−2 method. N = 5 is chosen as
the polynomial order in the current study. For the temporal discretization, the convection terms are treated explicitly with thesecond-order extrapolation scheme, while diffusion terms aretreated implicitly with the second-order backward difference
scheme. SEM is also known for its geometrical flexibility andprovides a setup, in which a body-fitted mesh on the sphericalmodel surface is generated with high-order accuracy utilizing
N th -order polynomials.

INTRODUCTION
A turbulent wake forms past a body moving in a fluid flow.
It is a canonical phenomenon, yet it impacts many realistic
applications, particularly in the transport industry. Thus it isbeneficial to understand the physics behind such flows. During the past few decades, turbulent wakes were investigated byboth experimental (Achenbach (1974); Berger et al. (1990))
and computational (Tomboulides & Orszag (2000); Yun et al.
(2006); Rodriguez et al. (2011)) tools. Characteristic motionsof the wake flow such as vortex shedding and bubble pumping were previously identified. In the current work, we intend
to characterize temporal spectra and spatial low-dimensional
flow patterns associated with these motions in a wake of asphere.
A wake flow past a sphere with a sub-critical Reynolds
number ReD = 3700 has been studied by high-fidelity simulations in recent years (Yun et al. (2006); Rodriguez et al.
(2011); Pal et al. (2017)). In this flow, the transition to turbulence occurs in the near wake, which gives rise to rich physical phenomena and interesting flow dynamics. This Reynolds
number is also low enough to allow for the use of Direct Numerical Simulations, making the chosen setup as a model forinvestigation in the current paper. Some of the previous highfidelity simulation studies with similar flow conditions focused

COMPUTATION SET-UP
The DNS of the flow past a sphere at ReD = 3700 is carried out in a cylindrical computational domain with x×r ×ϕ =[−5D, 15D] × [0, 10D] × [0, 2π], where D is the sphere diameter. The center of the sphere is located at the origin of thecoordinate system. The O-grid mesh with 9.8 million degrees
of freedom (DOFs) is utilized for the baseline simulation. The
mesh in the boundary layer around the sphere is refined with aminimum element height ∆rmin = 0.015D.
A Dirichlet boundary condition with uniform freestream
velocity u∞ is imposed at the inflow boundary. The outerboundaries of the cylindrical domain are treated with stabilized1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
outflow boundary conditions (Dong et al. (2014)). The no-slip
wall boundary condition is employed on the sphere surface.
A time step ∆t = 0.0016D/u∞ is used in this case with a
CFL number around 0.88 for the baseline case (Mesh 2). The
initial conditions for the simulation are set up as a uniform velocity and pressure with their corresponding free-stream values. In order to prevent the flow statistics from being affectedby the initial transient, the data for analysis is being collectedstarting at t = tmin = 200D/u∞ . The snapshots are saved every 125 time steps, corresponding to a period between snapshots of ∆tsnap = 0.2D/u∞ . As a result, a total number of 4400
snapshots (spanning over ∆tstat = 880D/u∞ ) are collected untilt = tmax = 1080D/u∞ .

(a)

RESULTSWe first demonstrate the grid convergence of the current
simulations and present the results of the validation against
the data available in the literature for the flow over a sphereat comparable Reynold numbers. For the grid convergencetest, two more meshes were generated: Mesh 1 with the lower
number of DOFs (2.8 mln.), and Mesh 3 with the higher number of DOFs (36.7 mln.), as compared to the baseline mesh
(Mesh 2 with 9.8 mln. DOFs). The results for the threemeshes were compared with the experimental data from Kim& Durbin (1988), and with the numerical data from Seidl et al.
(1998); Yun et al. (2006); Rodriguez et al. (2011); Bazilevs
et al. (2014); and Dorschner et al. (2016). All these referencesused ReD = 3700 as in the current simulations, except for Seidl
et al. (1998) who used ReD = 5000. Table 1 presents a comparison of the vortex shedding frequency (via St = f u∞ /D), thelength of the separation bubble (Lsep /D) , the drag coefficienton the sphere (Cd ), and the angle of the boundary layer separation along the sphere surface (ϕsep ). Figure 1 documents the
distribution of the pressure coefficient and the normalized wall
shear stress along the sphere surface, and Figure 2 presents themean streamwise velocity profiles at several spatial locationsin the sphere wake. It can be seen that Nek5000 results areconverged for Mesh 2, and a good agreement with the previous data is obtained for all the presented quantities.
We now proceed with analyzing the results of Mesh
2, having demonstrated the grid convergence. Temporallyaveraged and instantaneous contours of streamwise velocity
at a streamwise-radial symmetry plane at z = 0 are shown in
Figure 3. A recirculation region (also known as a separation
bubble) is formed in the near wake area. In the temporallyaveraged plot, the wake is symmetric about the centerline,
while in the instantaneous plot the wake shows a wavy pattern.
Two probes are placed in the near wake to further investigate the dynamics of this flow (Figure 3b). Probe 1, at(x, r, ϕ) = (3.07D, 0, 0), is located close to the edge of therecirculation region to capture the bubble pumping motion,
manifesting itself as a periodic shrinkage and expansion of
the separation bubble (Rodriguez et al. (2011)). Probe 2,
at (x, r, ϕ) = (3.07D, 0.54D, 0), is located in the shear layer,
where the vortex shedding occurs. In order to identify the characteristic frequency of these two motions, the power spectral
density is defined as
Z tmaxPSD (St) =tmin



−i 2π(St)u∞ t
φ (t) exp
D(b)
Figure 1: Distribution of the pressure coefficient and normalized wall shear stress along the sphere surface incomparison with available data.

The time series and the corresponding power spectral density of the two probes are plotted in Figure 4. Here φ is u0 , thefluctuating streamwise velocity. For probe 1, the bubble pumping motion is dominated by a low frequency St = 0.004. However, we should note that for such a low frequency motion, weonly capture 3.5 periods of it in the current simulations. The
next strongest peak is at St = 0.009 , which is approximately
half of St = 0.0178 identified by Rodriguez et al. (2011). The
reason for disagreement can be either an insufficient samplingtime for this low frequency motion (for St = 0.009, the current
DNS captures around 8 periods, while for St = 0.0178, DNS ofRodriguez et al. (2011) captures around 6 periods), or it can be
related to a difference in the length of the separation bubble,
see Table 1. A longer separation bubble of the current DNS“breathes” with a lower frequency than the shorter bubble in
Rodriguez et al. (2011). For probe 2, in the high frequency
regime, St = 0.22 is the dominating frequency, which is thevortex shedding frequency, and it agrees well with the previous studies (Rodriguez et al. (2011); Bazilevs et al. (2014);
Kim & Durbin (1988)). It is interesting to note that the Probe
2 location also feels the low-frequency motion at St = 0.003,
presumably associated with the bubble pumping.
The coherent structures are illustrated by the contours ofthe λci criterion in Figure 5. It shows that a helical structure
is formed in the wake, particularly downstream of the location where vortices are separated from the shear layer. We observe that, as time proceeds (from top to bottom in Figure 5),
the vortical structure indicated by the yellow arrow is advecting downstream without rotating about the centerline. Similarfindings are also reported in Yun et al. (2006) and Rodriguez
et al. (2011). At the same time, in order to form the helicalstructure, vortices should shed at different azimuthal locations,

2



dt .

(1)

We can further define the normalized PSD as
Normalized PSD(St) =PSD.
maxSt PSD(St)(2)

212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
Table 1: Comparison of St, Lsep /D, Cd , and ϕsep , for a flow over a sphere at ReD = 3700 between Nek5000 and otheravailable data. Dash indicates that the data is not available. ∗ DNS of Seidl et al. (1998) is performed with ReD = 5000.
All other data is for ReD = 3700. LES: Large Eddy Simulations; DNS: Direct Numerical Simulations; LBM: Lattice
Boltzmann Method.
StudiesSt

Lsep /DCd

ϕsep (◦ )Nek5000 Mesh 1

0.229

2.51

0.381

90.2

Nek5000 Mesh 2

0.219

2.67

0.377

89.9

Nek5000 Mesh 3

0.220

2.70

0.375

89.8

Kim & Durbin (1988) (exp)

0.225

-

-

-

Yun et al. (2006) (LES)

0.21

2.62

0.355

90.0

Rodriguez et al. (2011) (DNS)

0.215

2.28

0.394

89.4

Bazilevs et al. (2014) (DNS)

0.221

2.28

0.392

89.4

Dorschner et al. (2016) (LBM)

-

2.51

0.383

90.0

Seidl et al. (1998) (DNS∗ )

-

2.1

0.38

89.5

DMD mode with a frequency of St = 0.009 is more asymmetric, but both POD and DMD modes are able to capture themain pattern around the separation bubble. It suggests thatboth frequencies, St = 0.004 and St = 0.009, are involved in
the bubble pumping motion. The fact that one frequency is essentially twice as high as the other indicates that they may berelated as sub- and super-harmonics of the same process.

which must be induced by some rotating mechanism.
In order to ”extract” the characteristic motions, such asvortex shedding and bubble pumping, the modal analysis isconducted using the collected snapshots. The modal analysis isperformed in a box region identified by a dash line in Figure 3,
x × y × z = [0.5D, 10D] × [−1.5D, 1.5D] × [−1.5D, 1.5D], tofocus our attention on the wake area and to lower the computational cost. For that, the flow information is interpolated onto
the 101 × 51 × 51 = 262701 uniformly distributed grid points
in this box. We first build the matrix X ∈ Rn×m (Taira et al.
(2017)), where n = 788103 is the number of grid points times
3 for the number of dimensions, and m = 4400 is the numberof snapshots. In this matrix, we have all the three componentsof the fluctuating velocity at each location in each snapshot.
Based on this matrix, we perform POD and DMD decomposition using the techniques outlined in Taira et al. (2017), tostudy the characteristic motions in the wake flow.
The visualization of POD modes and the power spectral
analysis of their temporal coefficients is shown in Figure 6.
The first two modes are related to the vortex shedding motion with the same dominating frequency of St = 0.22. Fromthe spatial pattern of this mode pair, it is clear that the vortex shedding doesn’t start from the trailing edge of the sphere;
instead, it occurs more than 1.5D downstream of the trailingedge. The fifth mode is linked to the motions of the separation
bubble with a dominating frequency of St = 0.004. The mostenergetic portion of this mode is located around the separation
bubble. Thus, this mode reflects the motion of the whole separation bubble. The third and fourth modes (not shown here)
are also dominated by the vortex shedding frequency, but witha noisier pattern.
DMD modes with vortex shedding frequency and bubble
pumping frequency are depicted using the fluctuating streamwise velocity in Figure 7. The visualization of vortex shedding
motion shows almost identical patterns with the results fromthe first mode of POD. Figure 7(a) only presents the real partof this DMD mode, the imaginary part depicts similar patterns
with the second POD mode. It is because the first two POD
modes are dominated by the vortex shedding frequency with arelatively clean spectrum, while a DMD mode is only linked
with a single frequency. For the bubble pumping motion, the

CONCLUSIONS
Direct numerical simulations are carried out to investigate
the wake flow past a sphere at ReD = 3700. The main focusof this paper is to characterize the dynamics of coherent structures that are related to the characteristic motions of the wakeflow.
Before conducting analysis of this wake flow, we demonstrate the grid convergence of results and validate them against
both experimental and numerical data available in the literature. An overall good agreement is achieved.
For the analysis, we first look at the temporally-averaged
and instantaneous contours of the streamwise velocity. Twoprobes are placed in the near wake region to capture the vortex shedding and bubble pumping motions. A low frequency
St = 0.004 is identified at probe 1 as the dominating frequencyof the bubble pumping motion. For probe 2, the dominating
frequency is St = 0.22, which is the vortex shedding frequency.
The helical structure of the wake is illustrated by the contour of λci . A vortical structure is traced as the time advances. No obvious rotation about the geometric centerline
is observed; instead, the structure is merely translating downstream.
To visualize and further investigate the coherent structures that are related to the characteristic motions, we carry out
modal analysis using POD and DMD. Both POD and DMD
are able to capture the patterns of vortex shedding and bubble
pumping. The vortex shedding patterns from both methods arealmost identical. For the bubble pumping, the dominating frequency of the corresponding POD mode is St = 0.004, and thesecond peak is at at a higher frequency of St = 0.009, while
the DMD mode that is related to the bubble pumping has thefrequency of St = 0.009.

3

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

(a)

(a)

(b)
(b)

Figure 3: Temporally-averaged (top) and instantaneous(bottom) streamwise velocity contours on the plane z =0.

pressible flow simulations on severely-truncated unboundeddomains. Journal of Computational Physics 261, 83–105.
Dorschner, B., Frapolli, N., Chikatamarla, S. S. & Karlin, I. V.
2016 Grid refinement for entropic lattice Boltzmann models. Physical Review E 94 (5), 053311.
Fischer, P., Lottes, J., Kerkemeier, S., Marin, O., Heisey,K., Obabko, A., Merzari, E. & Peet, Y. 2015 Nek5000
User’s manual. Technical Report ANL/MCS-TM-351, Argonne National Laboratory, http://nek5000.mcs.anl.gov.
Kim, H. J. & Durbin, P. A. 1988 Observations of the frequencies in a sphere wake and of drag increase by acoustic excitation. Physics of fluids 31 (11), 3260–3265.
Pal, A., Sarkar, S., Posa, A. & Balaras, E. 2017 Direct numerical simulation of stratified flow past a sphere at a subcritical Reynolds number of 3700 and moderate Froude number.
Journal of Fluid Mechanics 826, 5–31.
Rodriguez, I., Borell, R., Lemkuhl, O. & Oliva, A. 2011 Direct
numerical simulation of the flow over a sphere at Re = 3700.
Journal of Fluid Mechanics 679, 263–287.
Seidl, V., Muzaferija, S. & Peric, M. 1998 Parallel DNS withlocal grid refinement. Appl. Sci. Res. 59, 379–394.
Taira, K., Brunton, S. L., Dawson, S. T. M., Rowley, C. W.,
Colonius, T., McKeon, B. J., Schmidt, O. T., Gordeyev, S.,
Theofilis, V. & Ukeiley, L. S. 2017 Modal analysis of fluid
flows: An overview. AIAA Journal 55, N. 12.
Tomboulides, A. G. & Orszag, S. A. 2000 Numerical investigation of transitional and weak turbulent flow past a sphere.
Journal of Fluid Mechanics 416, 45–73.
Yun, G., Kim, D. & Choi, H. 2006 Vortical structures behinda sphere at subcritical Reynolds numbers. Physics of Fluids
18, 015102.

(c)
Figure 2: Time-averaged streamwise velocity profiles at(a) x = 0.2D, (b) x = 1.6D, (c) x = 3D, in comparison
with available data.
Acknowledgements
This work is supported by NSF CBET grants # 1707075
and 1944568.

REFERENCES
Achenbach, E. 1974 Vortex shedding from spheres. Journal of
Fluid Mechanics 62, 209–221.
Bazilevs, Y., Yan, J., De Stadler, M. & Sarkar, S. 2014 Computation of the flow over a sphere at Re= 3700: A comparisonof uniform and turbulent inflow conditions. Journal of Applied Mechanics 81 (12), 121003.
Berger, E., Scholz, D. & Schumm, M. 1990 Coherent vortex
structures in the wake of a sphere and a circular disk at restand under forced vibrations. Journal of Fluids and Structures 4, 231–257.
Dong, S., Karniadakis, G. E. & Chrissostomidis, C. 2014 A
robust and accurate outflow boundary condition for incom-

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

(a)

(b)

(c)

(d)

Figure 4: Time series of u0 (left) and its power spectral analysis (right) at probe 1 (top) and probe 2 (bottom).

(a)

(b)

(c)

(d)

(e)

(f)
2

Figure 5: Contour of λci with the value of 0.5 (u∞ /D) from (left) view of the streamwise-radial plane; (right) view fromupstream. Time advances from top to bottom. Yellow arrow tracks a specific hairpin.

5

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

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: The first (top), second (middle) and fifth (bottom) POD modes on the z = 0 plane (left), and power spectral

analysis of each mode’s temporal coefficient (right). The leftmost edge of each snapshot is at the location of the trailingedge of the sphere (x/D = 0.5).

(a)

(b)
Figure 7: DMD modes with St = 0.22 (a) and St = 0.009 (b) on the z = 0 plane. The leftmost edge of each snapshot is at

the location of the trailing edge of the sphere (x/D = 0.5).

6