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

DIRECT NUMERICAL SIMULATION OF TURBULENT FLOW AROUND 3D
STEPPED CYLINDER WITH ADAPTIVE MESH REFINEMENT

Daniele Massaro
SimEx/FLOW, Engineering Mechanics
KTH Royal Institute of Technology
SE-100 44 Stockholm, Sweden
dmassaro@kth.seAdam Peplinski
SimEx/FLOW, Engineering Mechanics
KTH Royal Institute of Technology
SE-100 44 Stockholm, Sweden
adam@mech.kth.se

Philipp Schlatter
SimEx/FLOW, Engineering Mechanics
KTH Royal Institute of Technology
SE-100 44 Stockholm, Sweden
pschlatt@mech.kth.se

ABSTRACT

and it consists of two dominating shedding frequencies, labelled fS and fL , corresponding to the small and large cylinder vortex shed respectively. Initially in phase, the two wakecells are connected one by one across the interface. Due tothe different frequencies, they move out of phase and at least
one half-loop connection between oppositely rotating vortices
appears. For D/d > 1.55 the indirect interaction generates another distinct region behind the large cylinder. It was named
modulation zone by Lewis & Gharib (1992) and is characterized by its own frequency fN . According to these three dominating shedding frequencies, Dunn & Tavoularis (2006) identified three main regions: the S and L cells behind the small andlarge cylinder respectively, and the modulation cell N with thelowest shedding frequency ( fS > fL > fN ). The same nomenclature is adopted here. Dunn & Tavoularis (2006) investigated
the interactions among the cells for D/d = 2 at ReD = 150
through experiments and the results were confirmed throughsimulations based on Reynolds Averaged Navier-Stokes equations (RANS). On the one hand, it has been observed that theS–N cell boundary is stable and deflects spanwise into the large
cylinder direction; on the other hand, the N–L cell boundary is
unstable. The modulation region boundaries present complicated dynamics as the vortex dislocation occurs. This concept is adapted from solid mechanics and introduced first by
Williamson (1989). In the laminar vortex shedding regime
Tian et al. (2020) investigated the possible loop structures connection. Their observations seem to hold as ReD increases.
Indeed, the three main cells and the dislocation mechanism
are still detectable. Nevertheless, the literature lacks results as
the Reynolds number increases. Experimentally, Ko & Chan
(1992) measured the main integral quantities at ReD = 80000,
while Morton et al. (2016) presented several results up to
ReD = 5000 performing tomographic particle image velocimetry investigation. Regarding the numerical studies, there areunsteady Reynolds Averaged Navier-Stokes (URANS) simulations which have achieved ReD = 2000, e.g. Morton et al.
(2009). Being said that, just recently a first direct numericalsimulation has been performed at a high turbulent Reynolds
number by Tian et al. (2021), without fully resolving the wakeregion.

In the present study, we investigate the turbulent threedimensional flow around a stepped cylinder, namely two cylinders of different diameters joint at one extremity. We perform
a direct numerical simulation with the spectral element code
Nek5000 that uses a high order spatial discretisation (the polynomial order is p = 7). The adaptive mesh refinement technique is employed in the error-driven meshing procedure, allowing an adequately refined mesh everywhere. We consider
the Reynolds number ReD = 1000, based on the large cylinder
diameter and the uniform inflow velocity. We compare our results with the previous experimental campaign by Morton &Yarusevych (2014b). The results agree very well and we canidentify the three main wake regions: the S, N and L cell witha Strouhal number StS = 0.408, StN = 0.188 and StL = 0.201
respectively. The instantaneous mean flow properties are studied showing that the junction dynamics is more similar to theprevious laminar studies at ReD = 150 rather than at higher
ReD = 3900. Moreover, proper orthogonal decomposition isused to detect the most energetic coherent structures, that resemble the three wake cells.

INTRODUCTION
In the last years, the flow around a stepped (or step-) cylinder has gained more and more interest as it represents a usefulmodel for several engineering applications, e.g. the offshore
wind turbine towers. It consists of two cylinders with different diameters joint at one extremity and its fluid mechanics isfar from being fully understood. For example, three different
vortex shedding frequencies have been identified in its wake,counter-intuitively to what one could expect. This aspect has
turned out to have relevant side effects on the structural stability of the cylinder itself and further investigations are ongoing.
Unlike the uniform case, for the flow around a single step
cylinder, there is an additional parameter to take into account:the ratio D/d between the diameters of the large and small
cylinders. Lewis & Gharib (1992) reported a direct and indirect vortex interaction mode for the range 1.14 < D/d < 1.76
at 67 < ReD < 200. The direct mode occurs when D/d < 1.25

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
In the light of this, we have decided to perform a direct numerical simulation using the spectral element code Nek5000
comparing with the experimental results by Morton & Yarusevych (2014a,b). Rather than investigating the diameter ratio
and Reynolds numbers effect, we have preferred to focus on asingle case (ReD = 1000 and D/d = 0.5) and provide detailed
description of the flow structures. We use a high order polynomial spatial discretisation, guaranteeing an iso-parametricrepresentation of the cylinder. In addition to the high resolution around the cylinder, we use the adaptive mesh refinement (AMR) technique developed by our group, to guarantee
a DNS-like resolution in the wake as well. We first provide ageneral description through instantaneous flow visualisations
and measuring the main integral quantities. The mean flow
and the turbulent kinetic energy variations along the cylinder
are also described. Eventually, we extract the most relevantcharacteristic features of the flow dynamics using a modal decomposition, i.e. proper orthogonal decomposition (POD).

zx

Figure 1. The ratio between the volumetric grid spacing andKolmogorov length scale lvol /η is shown. We present theplanes xy at z = 0 and xz at y = 0 (zoomed-in window).

act solution to a system of partial differential equations and uN
is an approximate spectral-element solution with polynomial
order N. Let us expand u(x) on a reference element in termsof the Legendre polynomials: u(x) = ∑∞k=0 ûk Lk (x), where ûk
are the associated spectral coefficients and Lk (x) is the Legendre polynomial of order k. The error ε = ∥u − uN ∥L2 consistsof two contributions: a truncation error due to the finite number of coefficients in the spectral expansion and the quadratureerror. The estimated ε results in:

NUMERICAL METHOD
The numerical simulations have been performed using
the open-source CFD code Nek5000, Fischer et al. (2008).
Nek5000 is widely used in the CFD community to perform
high fidelity numerical simulations of transitional and highReynolds incompressible flows. Highly scalable and portable,
the code uses the spectral element method (SEM) by Patera(1984) as spatial discretisation. SEM offers minimal dissipation and dispersion, high accuracy and nearly exponential convergence. It can be seen as a high–order finite element method
where the solution is represented by Lagrangian interpolants
integrated over Gauss-Lobatto-Legendre points (GLL). We use
a staggered grid for the pressure to avoid spurious modes andguarantee the Ladyzhenskaya–Babuška–Brezzi (LBB) condition. Adopting the Pn − Pn−2 formulation, the pressure is discretised on Gauss-Lobatto (GL) points. For time integration,
third-order backward differentiation is applied, treating thenon-linear term explicitly with an extrapolation formula andthe remaining Stokes problem implicitly.

ε=û2N  21
dk
+
2N+1
2k+1

Z ∞ û(k)2
N

2

(1)

2

where we assume an exponential decay for the spectral coefficients of the form û(k) ≈ c exp(−σ k). The parameters c andσ are obtained via interpolating in a linear least-squares sense
the log(ûk ) w.r.t. k for k ≤ N. The online error measurement
ε provides an indication per element at each time step. Afterwards, it is averaged and the meshing procedure stops whenthe error indicator is converged. Once completed the adaptive
mesh design, the final grid is frozen and we start to collectdata.
To assess the quality of the mesh, we can evaluate the ratio between the local grid spacing and the Kolmogorov length
η. This latter is estimated as η = ν 3/4 /ε 1/4 , where the turbulent kinetic energy dissipation ε is computed via statisticalanalysis. With regard to the grid spacing, in each direction
we have: li /η, with lx , ly and lz that are the largest spacing
in x, y and z directions respectively. As an alternative, we canconsider the equivalent volumetric grid spacing lvol /η, wherelvol = (lx · ly · lz )1/3 . In Fig.1, we see how the sei-designed
mesh allows to capture the smallest spatial structures close to
the cylinder walls, up to Kolmogorov scale. We guarantee a
DNS-like resolution in the entire wake region, since at mostlvol /η ≈ 8. Please observe that the front part and the cylinder
boundary layer are laminar, so the η estimation is meaningless.

Adaptive mesh refinement
In the spectral element code Nek5000, our group has implemented the adaptive mesh refinement technique (AMR) andextensively used it, see Refs. Offermans et al. (2020, 2022);
Tanarro et al. (2020); Massaro et al. (2022). AMR provides asignificant flexibility improvement to the meshing procedure,
avoiding an over-refinement of the far-field. This results ina significant saving in the number of grid points. As a consequence, on the one hand we can consider a larger computational domain reducing the influence of the boundary conditions, which is critical for incompressible flow simulations.
On the other hand, we can increase the Reynolds number orstudy more complicated geometries. That is the case for the
stepped cylinder when an apriori knowledge of the solution is
not completely given and an error-driven meshing procedure
is pivotal.
In addition to an efficient refinement procedure (here we
perform the h-type of refinement, where the element is split
isotropically), the committed error has to be measured properly. In this regard we have used, the spectral error indicator(sei), which was introduced first by Mavriplis (1989). It measures the truncation and quadrature error for a specified target
variable. In the current case, after an in-depth methodology
study, the velocity solution field has been used. To explain how
it works, let us consider a 1D problem, where u(x) is the ex-

Flow configuration
The stepped cylinder is located in the origin of our systemof reference that is oriented according to the standard notationfor external flow around cylinders, i.e. with the z-axis along thecylinder. The diameter ratio is D/d = 2, with unitary diameter
D. To avoid side and boundary effects the computational box
domain is designed with adequate dimensions: (−10D, 30D)
in x, (−10D, 10D) in y and (−10D, 12D) in z, where x, y, z are
the streamwise, spanwise and vertical directions respectively.

2

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

H3

b)

H2
R

H1
y

c)

B
x

z

y

x

Figure 2. Instantaneous λ2 visualisations at the junction and in the wake. a) The λ2 = −4 structures are visualised on the junction(z = 0) and coloured with the streamwise vorticity (top view). b) The vortex dislocation at S-N interface is presented. The green
and orange vortices are shed at the negative and positive y respectively; the red lines indicate the oblique shedding occurring in themodulation region. c) The λ2 = −0.2 structures identify the symmetric base vortices on the back of the junction.
The inflow at x = −10D is set as a Dirichlet boundary condition with the streamwise velocity being a uniform distribution with unitary value U = u0 = 1. The outflow at x = 30D
consists of natural boundary condition (−pI + ν∇u) · n = 0,
where ν is the kinematic viscosity. Considering an infinity
long cylinder, symmetry boundary conditions are prescribed at
the top/bottom surfaces (u · n = 0 with (∇u · t) · n = 0). For thefront and back boundaries mixed conditions are used, similar
to the open boundary at the outflow, but prescribing zero velocity in non-normal directions. Eventually, for the step cylindersurfaces we have no-slip and impermeable wall.

ing separation generates only H1. This difference is confirmed
by looking at time-averaged streamlines in figure 3-b. That
has a relevant consequence in the H3 formation. Indeed, Tianet al. (2021) report that the edge vortex H3 is a horseshoe vortex caused by the backward flow below the junction vortex
H1. The H4-H1 vortex bridge and the hairpins observed at
ReD = 3900 are identified as responsible for H3. At the current Reynolds number, we assert that the mechanism which
generates the edge vortices agrees with the Dunn & Tavoularis(2006)’s observations at lower ReD = 150. The leading edge
separation with the small cylinder impingement drives the recirculating vortex H1 formation that laterally deviates the incoming flow shaping the H3s. The second type of structure
observed by Tian et al. (2021) is named H2 and is present at
the current ReD = 1000 too. These are horseshoe vortices thatevolve at the base of the small cylinder with cross-flow widths
that continue to increase downstream, as for H1 (see Fig.2-a).
All the three horseshoes (H1, H2, H3) in the forepart of thestep surface (x/D < 0), are stable and with small fluctuations
in time. On the contrary, the back structures are weaker andmore difficult to identify. Differently to Tian et al. (2021),
which uses time-average λ2 structures to observe them, wemanage to visualise the instantaneous base vortices B1 and B2(Fig.2-c). Eventually, Tian et al. (2021) points out a backside
vortex that here appears at a different location. We name it
ring-vortex R and its core is in the back part of the junction at(x/D, y/D, z/D) ∼ (0.39, 0, 0.02). R is formed by the H3 vortices that roll up around the edges and split before diverging
laterally. The vortex segment which remains attached to theedges generates R. The low-pressure region behind the cylinder stabilises R on the back, whereas the strong flow advection
interrupts the ring connection laterally.

INSTANTANEOUS FLOW FIELDSWe qualitatively describe the main structures present in
the flow field using the λ2 criterion by Jeong & Hussain(1995). The convention for the vortices classification follows
Tian et al. (2021).
First, we focus on the junction region. In the front part,the formation of streamwise vortices is similar to the top-end
plate of a finite-length cylinder, see Zdravkovich et al. (1989).
Figure 2-a shows that the impingement of the flow, due tothe presence of the small cylinder, together with the leadingedge separation is mainly responsible for the horseshoe vortex H1 formation. This generates a recirculating pressure
bubble with its core at (x/D, y/D, z/D) ∼ (−0.354, 0, 0.046).
Other typical structures are the edge vortices H3, always visible on the free-end surface and typically named tip-vortices.
The blockage and the H1 rotation induce the incoming flow
to divert laterally, where it spills over the edges of the step.
Here, it rolls up into two streamwise edge-vortices, see Fig.2a. These observations agree with the measurements first madeby Dunn & Tavoularis (2006) at lower ReD = 150. Recently,
Tian et al. (2021) points out some additional structures at
ReD = 3900. The first one is a vortex (H4) which appears
upstream of H3 and rotates in the same direction as H1. At
the current Reynolds (ReD = 1000) the mentioned mechanism
is not present on the junction front. The boundary layer separation occurs downstream and other smaller turbulent eddies
do not appear around the junction. Figure 2-a shows in a clear
way that the horseshoe vortex H4 is not present. The lead-With regards to the coherent structures in the wake region,
we observe several streamwise spaghetti-like structures on topof the typical vertical vortices (whose axis is z-oriented), see
Fig.2-b. In laminar vortex shedding, the z vortices are distinguishable making the vortex dislocations investigation simple. This is not the case for a turbulent wake. Nevertheless,since the S-N boundary is stable, it is possible to observe
when the dislocation occurs. Figure 2-b identifies the initialin phase vortex shed by the cylinder and the consequent dislo-3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Table 1. Strouhal numbers for different probe locations. The
results from the current study and Morton & Yarusevych(2014b) are reported on the first and second rows respectively.

St1

St2

St3

St4

St5

0.408

0.41 − 0.22

0.188

0.19 − 0.20

0.201

0.398

−

0.189

−

0.198

sion of the cell. An analogous observation can be deductedby our study: at the interface Fig.4-b the energy content is diminished substantially compared to that within the single cells
Fig.4-a, c. Nevertheless, a more detailed discussion about theenergy content for each cell is carried out in the POD section.

Integral quantities
In the current section we report the dimensionless dragcoefficient CD , i.e. the horizontal component of the aerodynamic force normalised with 1/2ρu20 DLz , with Lz being thevertical cylinder length. The aerodynamic force is computed
onlineR as the sum of the viscous and pressure contributions
Fa = S −pn + τn dδ S with S the cylinder surface, p the pressure, n the vector normal to the cylinder surface and τ theviscous stress τi j = µ∂ ui /∂ u j . After the normalisation, thedrag coefficient is time averaged to yield CD . For symmetry reasons, the lateral y component of the aerodynamic force,and its normalised formulation Cy , is expected to be zero for
an infinitely long simulation. In the current study they result
CyS = 0.0017 and CyL = 0.00053, assessing the well temporalconvergence for our statistics. The mean drag coefficients are
cD = 1.1084 and cD = 0.9083 for the small and large cylinder
respectively. Both the standard deviations are approximately
σcD ≈ 0.016. As expected in the ReD range 100 < ReD < 105
the drag coefficient is cD ≈ 1. The pressure contribution isdominant (around 90%) w.r.t. the viscous one for both cylinders.

cation. The dashed-red lines point out the oblique vortex shedding which creates the intermediate modulation region with atriangular-like shape. Contrary to the upper interface, the N-L
layer is unstable and makes it challenging to study. Nevertheless, we believe that is possible to find the N and L cells’ location and monitor the shedding through a modal decomposition
analysis. In the current study, we rely on the proper orthogonal
decomposition (POD) of the flow to identify the most energetic
structures in the wake and characterize the different shedding
regions. The results are discussed below.

QUANTITATIVE ANALYSISThe quantitative description of the flow begins with measurements of the velocity signal from virtual probes located in
the flow and the evaluation of the relevant integral quantities.
The probes are placed at specific locations, following Morton& Yarusevych (2014b)’s experimental setup. We also measure the drag coefficient of each cylinder and provide a brief
overview of some relevant statistical quantities.

Mean flow quantities
The mean streamwise velocity profiles are shown in
Fig.3-a at different z locations. The cylinder is immersed ina uniform velocity field that forms a laminar boundary layer
and eventually separates. The flow evolves detached from thesurface due to an adverse pressure gradient. The phenomenon
occurs similarly for the large and small cylinders. On the junction, the mean streamwise velocity shows a sharp contractionof the layer before the separation. In a similar way to the freeend uniform cylinder, the downwash region is visible behindthe junction, see Fig.3-b. In the front junction part, the mean
statistics confirm what has been observed via instantaneous λ2
structures: the flow is lifted on the top lead by the leadingedge vortex H1. Then, the impingement at the base of thesmall cylinder pushes the flow downstream on the junction.
The core of the recirculating vortex H1 is located at(x/D, y/D, z/D) = (−0.354, 0, 0.046). In agreement with
Dunn & Tavoularis (2006) and contrary to Tian et al. (2021),
H1 is the only vortex present in the front part, see the mean velocity streamlines in 3-b. The H2 vortex cores are not visible
here because the cross-section is taken for y = 0 and H2 structures are developed on the sides. In the back, the streamlines
confirm that no backside vortex is present. In fact, the ring vortex is visible and located at (x/D, y/D, z/D) = (0.39, 0, 0.02),
see Fig.3-b.
Figure 3-c shows a spatial map of the mean turbulent
kinetic energy k, which is proportional to the trace of the
Reynolds stress tensor. We observe the symmetry w.r.t. to
the streamwise direction supporting the adequacy of the statistical sample. To the authors’ knowledge, this is the first-everk-map for this flow configuration. The white line defines theiso-contour k = 0 and points out the region where the flow islaminar. The “bat”-shape varies with the vertical positions, becoming more irregular at the junction. Since ReD = 2Red weobserve different tails’ elongation, but both have a similar longshape with a more extended laminar region for the large cylin-

Velocity probes measurements
The previous studies indicate that three distinct vortex
shedding cells develop in the wake and eachone is characterised by different a frequency. We use spectral analysis ofvelocity signals to determine the vortex shedding frequency
and investigate how it varies for different vertical locations.
The spectral analysis results are compared with previous
experiments carried out by Morton & Yarusevych (2014b) at
the same ReD = 1000. In this regard, we locate five virtual
velocity probes at the same locations. The time signals are acquired with the sampling frequency fs = 100Hz. It definitely
meets the Nyquist theorem requirements, since the highest frequency of the flow is supposed to be the small cylinder vortex shedding at fs ≈ 0.4Hz. The time-averaging interval islong enough to capture the slowest dynamics exploring over≈ 200 S-cell vortex shedding periods. The virtual probes arelocated at z1 /D = 4, z2 /D = −1, z3 /D = −3, z4 /D = −5 andz5 /D = −9. The streamwise and spanwise locations are thesame for all of them, x/D = 5 and y/D = 0.75 respectively.
The resulting streamwise velocity spectra are presented in Fig.
4.
The Strouhal number is defined as St = f D/u0 and theagreement with the experimental measurements is excellent.
Each wake region is characterized by a distinct shedding frequency and exhibits a single peak: StS = 0.408, StN = 0.188
and StL = 0.201 for the S, N and L-cell respectively (see Fig.4a, b, d). At the boundary between the S and N cells (z2 ), twodominant peaks are related to the Strouhal numbers that correspond to the adjacent cells (Fig.4-c). Similarly, we observea double peak at the N-L boundary. Morton & Yarusevych(2014b) uses the velocity streamwise spectra to analyse how
the energy context varies in the different cells and the exten-

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

a)

c)

b)

z
x

y
x

Figure 3. Selection of time-averaged quantities. a) The mean streamwise velocity is presented at different z location, from the top
z/D = 9, z/D = 0 and z/D = −9. b) The mean vertical velocity component is shown. The in-plane streamlines for the mean x velocity
component are illustrated in the front/back junction regions. c) The mean turbulent kinetic energy map is shown at z/D = 9, z/D = 0.01,
z/D = 0 and z/D = −9.

Figure 5. The percentage energy contribution of each k modefor different number of snapshots n.

Figure 4. The power spectral density of the streamwise velocity component at locations specified above in the text: z1 ,z2 , z3 and z5 . The log-log plot reports the Strouhal number
along the x axis.
ICTAM2020+1

Dunn, W. et al., Experimental studies of vortices shed from cylinders with a step-change in diameter. JFM, 2006

7

(1987) to handle the analysis of three-dimensional flows witha big set of grid points.
We collect a large number of snapshots that are sampledaccording to the Nyquist criterion and for a time interval sufficiently long to embed the slowest frequency. The convergenceis observed in Figure 5. According to the adopted convection,mode 0 corresponds to the expected mean value of the flow
and the agreement with the estimated mean value is excellent.
The first pair of travelling modes, namely the modes 1 and 2
which have 90-degree phase shift, correspond to large cylinder modes in the L cell region. They account each for aroundthe ∼ 7% of the total (fluctuating) energy. The second pair,modes 3 and 4, identify the S region behind the small cylinder
contributing each for the ∼ 5.5% of total energy. The N-cellregion (modes 5 and 6) have an energy contribution that is halfof the previous pair ∼ 2.8%. The results agree with a previ-der. At z = 0, we have a triple peak structure that gets evenmore contorted on the junction surface (z = 0.01).

PROPER ORTHOGONAL DECOMPOSITION
The POD performs the modal decomposition of the nonlinear flow minimizing the residual energy between the snapshots and its reduced linear representation. They are constructed to be the orthogonal basis which represents the optimal projection, in energy terms, of the most important dynamical features. It was first introduced in the fluid-mechanics context by Lumley (1970) and the approach is referred as classical
POD. Here we use the snapshot POD, developed by Sirovich

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
REFERENCES
Dunn, W. & Tavoularis, S. 2006 Experimental studies of vortices shed from cylinders with a step-change in diameter. J.
Fluid Mech. 555, 409–437.
Fischer, P., Kruse, J., Mullen, J., Tufo, H., Lottes, J. & Kerkemeier, S. 2008 Nek5000: open source spectral element CFD
solver. https://nek5000.mcs.anl.gov/ .
Jeong, J. & Hussain, F. 1995 On the identification of a vortex.
J. Fluid Mech. 285, 69–94.
Ko, N. W. M. & Chan, A. S. K. 1992 Wakes behind circular cylinders with stepwise change of diameter. Exp. Therm.
Fluid Sc. 5, 182–187.
Lewis, C.G. & Gharib, M. 1992 An exploration of the wakethree dimensionalities caused by a local discontinuity incylinder diameter. Phys. Fluids 4, 104–117.
Lumley, J. L. 1970 Stochastic tools in turbulence. Academic .
Massaro, D., Peplinski, A. & Schlatter, P. 2022 Interface discontinuities in spectral-element simulations with adaptive
mesh refinement. In Spect. High Ord. Meth. Part. Diff. Eq.
- ICOSAHOM 2021 (submitted). Springer.
Mavriplis, C. 1989 Nonconforming discretizations and a posteriori error estimators for adaptive spectral element techniques. PhD thesis, MIT, USA.
Morton, C. & Yarusevych, S. 2014a On vortex shedding from
low aspect ratio dual step cylinders. J. Fluids Struct. 44,
251–269.
Morton, C. & Yarusevych, S. 2014b Vortex dynamics in the
turbulent wake of a single step cylinder. J. Fluids Eng.
136 (3).
Morton, C., Yarusevych, S. & Carvajal-Mariscal, I. 2009
Study of flow over a step cylinder. Appl. Mech. Mater. 15.
Morton, C., Yarusevych, S. & Scarano, F. 2016 A tomographic
particle image velocimetry investigation of the flow development over dual step cylinders. Phys. Fluids 28.
Offermans, N., Massaro, D., Peplinski, A. & Schlatter, P. 2022
Error-driven adaptive mesh refinement for unsteady turbulent flows in spectral-element simulations. Comp. Fluids(submitted) .
Offermans, N., Peplinski, A., Marin, O. & Schlatter, P. 2020
Adaptive mesh refinement for steady flows in nek5000.
Comp. Fluids 197, 104352.
Patera, A.T. 1984 A spectral element method for fluid dynamics: Laminar flow in a channel expansion. J. Comput.
Physics 54 (3), 468–488.
Sirovich, L. 1987 Turbulence and the dynamics of coherent
structures. part i: Coherent structures. Q. Appl. Math. 45.
Tanarro, A., Mallor, F., Offermans, N., Peplinski, A., Vinuesa,
R. & Schlatter, P. 2020 Enabling adaptive mesh refinement
for spectral-element simulations of turbulence around wing
sections. Flow, Turb. Combus. 105, 415–436.
Tian, C., Jiang, F., Pettersen, B. & Andersson, H. I. 2020 Vortex dislocation mechanisms in the near wake of a step cylinder. J. Fluid Mech. 891, A24.
Tian, C., Jiang, F., Pettersen, B. & Andersson, H. I. 2021 Vortex system around a step cylinder in a turbulent flow field.
Phys. Fluids 33.
Williamson, C. H. K. 1989 Oblique and parallel modes ofvortex shedding in the wake of a circular cylinder at low
reynolds number. J. Fluid Mech. 206 (3), 579–627.
Zdravkovich, M., Brand, V., Mathew, G. & Weston, A. 1989
Flow past short circular cylinders with two free ends. J.
Fluid Mech. 203, 557–575.

Figure 6. Velocity magnitude iso-surfaces for the first (orange), third (blue) and fifth (green) POD mode.

ous study by Morton & Yarusevych (2014a) where the POD
has been performed based on planar velocity data obtained via
Particle Imagine Velocimetry (PIV). For the mode 5 we do notsee any interactions between S and N cells, which are characterized by a stable interface; unlike the N-L layer. Moreover,the POD helps to identify the modulation region extension andits shape. This latter agrees well with the visualisations carried out by Dunn & Tavoularis (2006), showing a triangle-likeshape that is persistent in the turbulent regime. Finally, we
would like to comment on the single mode 7. It has an energy
contribution larger than the 1% and seems related to the downwash phenomenon which occurs behind the junction. It has afrequency that agrees with the previous estimation by Morton& Yarusevych (2014b). As we are aware that each POD mode
does not correspond to a single frequency in time, further investigations are ongoing through other modal decompositions,such as dynamic mode decomposition (DMD), which allows
having every mode with a single temporal frequency.

CONCLUSIONSThe study performs a DNS of a turbulent flow around thethree-dimensional stepped cylinder at ReD = 1000 by meansof adaptive mesh refinement (AMR). The AMR enables tohave an accurate solution for the entire domain, with a volumetric grid spacing less than 10 Kolmogorov length scales
everywhere. The results are compared with the experiments
by Morton & Yarusevych (2014b) showing an excellent agreement. The three S, N and L wake regions are confirmed in
the turbulent regime with a stable interaction at the S-N interface, but unstable at the N-L boundary. The proper orthogonal
decomposition (POD) identifies these three coherent regions
and allows for assessing their energy contributions. In addition, even if less energetically relevant, the downwash mode
is observed. We also discuss some of the mean flow properties. The mean turbulent kinetic energy k shows how the flow
evolves differently around the cylinder according to the z locations. In the future, a more insightful and detailed statisticaldescription, e.g. considering the k budgets, might help to improve our understanding of such a flow configuration.

ACKNOWLEDGEMENTS
This study was funded by Knut and Alice Wallenberg Foundation and the Swedish e-Science Research Centre(SeRC); enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC).

6