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

LARGE EDDY SIMULATIONS OF A TURBULENT BOUNDARY LAYERWITH UNSTEADY PRESSURE GRADIENTSFrancesco Ambrogi
Dept. of Mechanical and Materials Eng.
Queen’s University
Kingston, ON Canada
f.ambrogi@queensu.ca

Ugo Piomelli
Dept. of Mechanical and Materials Eng.
Queen’s University
Kingston, ON Canada
ugo@queensu.ca

David E. Rival
Dept. of Mechanical and Materials Eng.
Queen’s University
Kingston, ON Canada
d.e.rival@queensu.ca

ABSTRACT

tially. However, in many cases (e.g., helicopter blades, turbine
blades, pitching airfoils etc.) the pressure gradient varies both
spatially and temporally. In fish locomotion, for instance, unsteady pressure gradients are the norm: fish propel themselves
by waving their tail, and (especially during maneuvers) thechanging body curvature generates unsteadiness in the pressure distribution along the body, and pressure gradients thatare space and time dependent (Domenici et al., 2004).

We use large-eddy simulations to investigate a spatiallydeveloping turbulent boundary layer with a space- and timedependent pressure gradient. The flow oscillates between
times in which a favorable pressure gradient is followed by anadverse one, a zero-pressure gradient, and a period in which
the adverse pressure gradient is followed by a favorable one.
The alternating favorable and adverse pressure gradient causes
the flow to separate periodically and reattach from the wall.
Several cases have been investigated for a range of reduced
frequencies, spanning the range between a very rapid flutterlike motions to a slower quasi-steady flapping. Time- andphase-averaged fields are analyzed, and comparison is madewith steady cases with fixed pressure gradients.

While the steady flow in APG and FPG TBLs has beenstudied extensively (see the papers referenced above, andmany others) fewer investigations of the unsteady case canbe found (Schatzman & Thomas, 2017). In the present paper
we consider an unsteady spatially developing FPG/APG turbulent boundary layer, which was simulated using wall-resolved
large-eddy simulations (LES). By oscillating the freestream
velocity periodically in time, we were able generate a flow
field that can either separate or tend towards re-laminarization,
depending on the phase of the cycle. Various frequencies ofoscillations are considered. The Reynolds number is moderate, but large enough to result in a fully developed equilibrium
TBL if no pressure gradient is applied. In addition, we also
performed steady calculations with time-independent pressure
gradients at selected phases, for comparison with the dynamic
cases.

INTRODUCTION
Almost all turbulent boundary-layer flows of practical relevance are subject to pressure gradients, usually introduced bygeometric changes. The pressure gradient may be favorable(FPG) or adverse (APG) depending on whether the pressure
decreases or increases along the flow direction. Both conditions are very common in engineering and physical systems
such as airfoils or turbine blades, and in geophysical applications.
If the APG is strong enough, the boundary layer may
separate from the wall. Flow separation can significantly affect the efficiency of fluidic devices. A comprehensive review
on the physics of flow separation can be found in Simpson(1989). A strong FPG, on the other hand, may lead to a relaminarization of the turbulent boundary layer (TBL) causedby the predominance of pressure forces over the turbulent
Reynolds stresses (Narasimha & Sreenivasan, 1973). Morecomplete reviews of the flow in FPG and APG boundary layers (both experiments and numerical simulations) can be foundin papers by Fernholz & Warnack (1998) and Wu & Piomelli(2018) for APG TBL, and Bourassa & Thomas (2009); Piomelli & Yuan (2013) for the FPG case.
Most studies of APG and FPG boundary layers consider pressure-gradients that are constant or vary only spa-

Although much effort has been invested to shed light into
the complexity of separated flows, many questions still remainunanswered: (1) How is the reduced frequency affecting theseparation cycle of a flat-plate TBL? (2) What are the underlying physical characteristics behind the phenomenon of dynamic hysteresis and their effects on the flow behavior? In thepresent work, results will be presented to answer these questions.

NUMERICAL METHODOLOGY
In the present work, simulations are performed using theLarge Eddy Simulation (LES) technique. In LES, the incompressible Navier Stokes equations are solved for filtered quan-112th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
tities (here indicated with an overline):
∂ ui
= 0;
∂ xi


∂ τi j∂ ui∂∂pui u j = −++ ν∇2 ui −
(1)
∂t
∂xj
∂ xi
∂xj

where, (x, y, z) are the streamwise, wall-normal, and spanwise
directions, and τi j = ui u j − ui u j is the sub-filter scale stresstensor. In the present study τi j is modelled using the Vreman eddy-viscosity model (Vreman, 2004). We use a wellvalidated finite-difference code that has previously been applied to similar cases (Keating et al., 2004; Yuan & Piomelli,
2015; Wu & Piomelli, 2018). The dimensions of the computational domain are Lx = 600δo∗ , Ly = 64δo∗ , and Lz = 55δo∗ .
δo∗ is the boundary layer displacement thickness at the inflow
plane. The Reynolds number based on δo∗ and the freestream
velocity at the inflow (which is constant in time) is Reδo∗ =1, 000. The grid uses Nx × Ny × Nz = 1536 × 192 × 256 points;a grid-convergence study was performed to verify that thisresolution is adequate. In wall units (defined using the wall
shear velocity uτ at the inflow plane) we have ∆x+ = 18.7,+∆y+min = 0.7, and ∆z = 10. A periodic boundary condition isused in the spanwise direction, a convective boundary condition at the outflow plane, and a velocity obtained from a parallel auxiliary simulation is used at the inflow.
The pressure gradient is generated by imposing a verticalvelocity V∞ (x,t) at the free-stream that changes both in space
and time, V∞ (x,t) = Vo (x) sin(2πt/T ), where T is the oscillation period, and Vo a streamwise distribution of wall-normal
velocity, which was chosen to match the case studied by Na &Moin (1998). The streamwise freestream velocity is obtainedby imposing a zero vorticity condition (Lund et al., 1998; Na& Moin, 1998; Wu & Piomelli, 2018).
The non dimensional parameter that characterizes the unsteadiness is the reduced frequency k = π f LPG /Uo , wheref = 2π/T is the imposed frequency, LPG is a characteristic
length (we use the length over which the Pressure Gradientvaries), and Uo is the free-stream velocity at the inflow plane.
We performed numerical simulations for k = 0.2, 1 and 10.
Leishman (2006) found that a reduced frequency k > 0.05 was
the threshold beyond which the boundary layer was clearly unsteady, and many experimental studies have been carried out
for a wide range of reduced frequency 0.1 < k < 82 (Karlsson,
1959; Brunton & Rowley, 2009). In our case the definition of
the length-scale is not unique, and we chose the length of theregion where the pressure-gradient is significant as the one thatmost closely resembles the chord length of an airfoil. Becauseof this arbitrariness, a one-to-one comparison with the studies
in the literature is difficult; however, we obtained trends thatare in agreement with the results in the literature, as will beshown momentarily.
The oscillation period was divided into 20 equally spaced
phases, and phase averaging (here indicated with ⟨·⟩) has been
carried out for every quantity. As shown in Figure 1, if
φ is the phase angle, φ = 0o and 180o degrees correspondto a zero-pressure gradient (ZPG) case, whereas at 270o an
APG is followed by an FPG, and at 90o the FPG precedes
the APG. We also carried out numerical simulations witha steady pressure gradient corresponding to that imposed at
φ = 0◦ , 54◦ , 90◦ , 270◦ , and 306◦ .

Figure 1. Freestream velocity at four phases in the cycle.
Black arrows denoting the direction of increasing phase angle Φ.

Figure 2. Contours of phase-averaged streamwise velocity⟨u⟩ for the phase Φ = 270◦ . (a) k = 10; (b) k = 1; (c) k = 0.2;(d) Steady calculation.

Case

XS /δo∗

H/δo∗

L/δo∗

XR /δo∗k = 10

251

1.4

124

375

k=1

253

3.2

67

320

k = 0.2

240

25

120

360

Steady

247

23

122

369

Table 1. Dimensions of the four separation bubbles for theunsteady and steady cases at Φ = 270◦ .

phase-averaged steamwise velocity ⟨u⟩ for the extreme phaseΦ = 270◦ . We can observe that flow separation (on average)
occurs in each of the unsteady cases though the size of the separation region greatly changes as the frequency k decreases.
For the high and medium frequency cases (k = 10 and k = 1)
the separation region does not have the time to grow as thickas in the corresponding steady case (Figure 2 d), however it
is very interesting to note that the length of the separation region for the k = 10 case is comparable with the corresponding
steady case, but the length of the separation region for the k = 1
case is significantly lower. This behavior will be explained in
the following section. As the frequency k is reduced to 0.2,
a quasi-steady state is approached, and the dimensions of theseparation bubble (both height and length) match very well thesteady calculation. Table 1 summarizes the dimensions of theseparation bubbles; here XS and XR are the averaged separation
and reattachment locations in the streamwise direction.
Figure 4 shows the behavior of the phase-avergaed
streamwise velocity ⟨u⟩ for the medium frequency case k = 1.
As discussed above, as the frequency decreases, the thickness

RESULTSFLOW EVOLUTION
The first important result is the effect of the reduced frequency k on flow separation. Figure 2 shows contours of2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
in Figure 3) includes phases in the separation side of the cycle,
where the APG precedes the FPG, and the freestream velocity U∞ < 1. The second (black square in Figure 3) includes
the two phases characterized by a ZPG. The third group (redsquare in Figure 3) includes the phases in the acceleration sideof the cycle, where the FPG precedes the APG, and U∞ > 1.
In order to visualize the hysteresis effects we can observe
the behavior of the velocity profile at a selected streamwise
location for every phase in the cycle. If there were no hysteresis effects, the velocity profile for matching phases shouldbe the same. Figure 6 shows results for the high frequency
case k = 10 in the streamwise location x/δo∗ = 300. Due to thehigh frequency (large k), the convective time-scale of the flow
(LPG /Uo ) is dominant over the unsteady imposed one (1/π f ).
The velocity at the matching phases is in good agreement,
except for a very thin layer close to the wall. Hysteresis effects are therefore concentrated very close to the wall, where
viscosity causes a small lag in the response to the freestream
acceleration/deceleration. We can also observe that, at phaseΦ = 198◦ , corresponding to the instant when the freestream
forcing (V∞ velocity profile) turns from blowing-suction tosuction-blowing, the flow separates until it reaches the extreme
phase Φ = 270◦ . At this frequency the flow is synchronizedwith the freestream forcing everywhere except for a very small
near-wall region where friction generates dynamic hysteresis.
For the k = 1 case, Figure 7, at the same streamwise location, hysteresis effects are moving away from the wall towards the center of the separated shear-layer. The convective
time-scale of the flow and the unsteady imposed one are now
comparable, resulting in a significant lag in response of thenear-wall fluid to the freestream forcing. In fact, coming backfrom the acceleration side (close to Φ = 90◦ ) and moving towards separation, although the freestream velocity U∞ is decelerating, the velocity is significantly higher compared to thecorresponding matching phase, and the flow does not separate
until Φ = 234◦ when the magnitude of V∞ reaches 70% of its
maximum value.
As the frequency is further decreased (k = 0.2), at thesame streamwise location, the unsteady imposed time-scale
is now dominant over the convective one, and given the very
long oscillation period T , the flow has time to adapt to thefreestream perturbations. Although the flow is moving towards
a quasi-steady state as shown in Figure 5, hysteresis effects
still play a significant role in the domain. However, a crucial
difference from the other two cases is that, at this frequency,
dynamic hysteresis effects are confined to the separation sideof the cycle whereas all the matching phases in the accelerationside show very good agreement.

Figure 3. Schematic representation of all matching phases in
one complete cycle. Following the black solid line from theextreme phase Φ = 270◦ the cycle reaches the opposite extreme Φ = 90◦ and goes back following the black dashed line
towards completion of the cycle.

of the separation bubble increases but the length of the separation region is significantly reduced (see phase Φ = 270◦
Figure 4). Moreover, the separation region is highly unstable,
and the slow-moving fluid, generated by the flow reversal, isadvected downstream and periodically washed out of the domain (Simpson, 1989). This peculiar flow dynamics has only
been observed for the k = 1 case, which indicates that there isa limiting frequency k beyond which the recirculation vortex
is shed. This behavior is consistent with experimental observations of a pulsating flow over a backward-facing step. A
critical frequency F ∗ = 0.07, defined using the step height hand the mean velocity U, was found, after which the shedding
occurred (Mullin et al., 1980; Simpson, 1989).
Finally Figure 5 shows phase-averaged streamwise velocity ⟨u⟩ profiles for 4 phases in the cycle at 3 different streamwise locations in the domain: x/δo∗ = 270, 300, 450 corresponding to locations upstream, at the center, and downstreamof the pressure gradient region. At Φ = 90◦ , in the acceleration side of the cycle where U∞ > 1, the dynamic cases matchvery well the steady calculations (black dots in Figure 5), withthe best agreement being the lowest frequency k = 0.2. At
Φ = 270◦ on the other hand, in the separation side of the cycle
where U∞ < 1, great differences arise. First, we observe thedifference between the dynamic cases, which is due to the different size and development of the separation bubbles. Second,even though the case k = 0.2 is approaching a quasi-steadystate, displaying a good match with the steady calculation, it ischaracterized by a stronger backflow in the center of the separation region indicating the persistence of transient effects.
The most interesting case however, can be visualized bycomparing Φ = 0◦ and Φ = 180◦ in Figure 5. At these twophases, the boundary layer is experiencing the same pressure
gradient, specifically zero in both cases. The main difference
between the two phases is that in the first case (Φ = 0◦ ) theflow is going towards the acceleration side of the cycle, in thesecond case, the flow is coming from the acceleration side andgoing towards separation. Very evident from Figure 5 is thatwhile k = 10 (very high) and k = 0.2 (very low) frequencies
match reasonably well the steady calculation, k = 1 case displays discrepancies in the velocity profiles. This suggests thatthe advection of the stalled fluid region downstream of the separation bubble is generating significant hysteresis effects.

CONCLUDING REMARKS
Large-eddy simulations of a spatially-developing turbulent boundary layer with space- and time-dependent pressure
gradients have been carried out. The unsteadiness was imposed by oscillating the freestream vertical velocity V∞ periodically. The maximum adverse pressure gradient was strong
enough to generate a large-scale turbulent separation bubble inthe domain. Simulations have been performed for a range ofreduced frequency k = 10, 1, 0.2, and results have been compared with steady calculations at the same instantaneous pressure gradient.
Two main questions have been driving our research: 1)
How is the reduced frequency k affecting the separation cycleof a flat-plate TBL? (2) What are the underlying physical characteristics behind the phenomenon of dynamic hysteresis andtheir effects on the flow behavior?

HYSTERESIS EFFECTS
Dynamic hysteresis is observed when a physical quantity
assumes two different values at corresponding phases in a periodic cycle. As shown by the schematic representation in Figure 3, dividing the cycle in 20 equally spaced phases results
in a total of 9 matching phases in which the boundary layer
is experiencing the same pressure gradient twice in one cycle.
The phases are divided into three groups: the first (blue square312th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
favourable pressure gradient and of the reynolds number on
an incompressible axisymmetric turbulent boundary layer.
part 1. the turbulent boundary layer. J. Fluid Mech. 359,
329–356.
Karlsson, S. K. F. 1959 An unsteady turbulent boundary layer.
J. Fluid Mech. 5 (4), 622–636.
Keating, A., Piomelli, U., Bremhorst, K. & Nešić, S. 2004
Large-eddy simulation of heat transfer downstream of a
backward-facing step. J. Turb. 5 (1), 020.
Leishman, G. J. 2006 Principles of helicopter aerodynamics.
Cambridge U. P.
Lund, T. S., Wu, X. & Squires, K. D. 1998 Generation of inflow data for spatially-developing boundary layer simulations. J. Comput. Phys. 140, 233–258.
Mullin, T., Greated, C. A. & Grant, I. 1980 Pulsating flow overa step. Phys. Fluids 23 (4), 669–674.
Na, Y. & Moin, P. 1998 Direct numerical simulation of a separated turbulent boundary layer. J. Fluid Mech. 374, 379–
405.
Narasimha, R. & Sreenivasan, K. R. 1973 Relaminarization inhighly accelerated turbulent boundary layers. J. Fluid Mech.
61 (3), 417–447.
Piomelli, U. & Yuan, J. 2013 Numerical simulations of spatially developing, accelerating boundary layers. Phys. Fluids 25, 101304–1—21.
Schatzman, D. M. & Thomas, F. O. 2017 An experimental investigation of an unsteady adverse pressure gradient turbulent boundary layer: embedded shear layer scaling. J. Fluid
Mech. 815, 592.
Simpson, R. L. 1989 Turbulent boundary-layer separation.
Ann. Rev. Fluid Mech. 21 (1), 205–232.
Vreman, A. W. 2004 An eddy-viscosity subgrid-scale modelfor turbulent shear flow: Algebraic theory and applications.
Phys. Fluids 16 (10), 3670–3681.
Wu, W. & Piomelli, U. 2018 Effects of surface roughness on
a separating turbulent boundary layer. J. Fluid Mech. 841,
552–580.
Yuan, J. & Piomelli, U. 2015 Numerical simulation of a spatially developing accelerating boundary layer over roughness. J. Fluid Mech. 780, 192–214.

The reduced frequency k plays a crucial role in the flow
configuration and the characteristics of the transient flow separation. If k is large enough (k = 10) the separation bubble(formed at the center of the domain) does not have time togrow as thick as in the corresponding steady case, but thelength of the separation region is comparable. For the medium
frequency case, on the other hand, we observed a growth of
the separation bubble in the wall-normal direction, but thelength of the separation region significantly decreased compared to the steady case. At this frequency the separation
region is highly unstable and is periodically advected downstream, and washed out of the domain (Mullin et al., 1980;
Simpson, 1989). Although the lowest frequency case k = 0.2
shows a trend towards a quasi-steady state, transient effects
persist in the separation region.
Dynamic hysteresis was observed in all unsteady cases.
Specifically, when the convective time-scale of the flow dominates over the imposed unsteady one (k = 10), hysteresis effects are confined to the near-wall region, and the flow behavior is synchronized everywhere with the freestream forcing.
As the frequency decreases, hysteresis effects move away fromthe wall towards the center of the separated shear-layer. Whenk = 1, the two time-scales are comparable and the response of
the near-wall fluid lags significantly the freestream forcing. In
the lowest frequency case, the unsteady time scale is now dominant over the convective one, and the flow has time to adapt
to the freestream forcing, however hysteresis effects are stillobserved in the separation side of the cycle.

REFERENCES
Bourassa, C. & Thomas, F. O. 2009 An experimental investigation of a highly accelerated turbulent boundary layer. J.
Fluid Mech. 634, 359–404.
Brunton, S. & Rowley, C. 2009 Modeling the unsteady aerodynamic forces on small-scale wings. AIAA Paper 2009-1127.
Domenici, P., Standen, E. M. & Levine, R. P. 2004 Escapemanoeuvres in the spiny dogfish (squalus acanthias). J. Exp.
Biol. 207 (13), 2339–2349.
Fernholz, H. H. & Warnack, D. 1998 The effects of a4

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

Figure 4. Contours of phase-averaged streamwise velocity ⟨u⟩ for k = 1. Only the main four phases of the cycle are shown: Φ = 0◦ ,
Φ = 90◦ , Φ = 180◦ , and Φ = 270◦ .

Figure 5. Streamwise phase-averaged ⟨u⟩ velocity profile for four phases in the cycle for the different reduced frequencies k (colours)
and streamwise locations (line styles). Comparison is made with steady calculations (symbols) at the same streamwise locations. Eachprofile is shifted by one unit for clarity. Solid line: x/δo∗ = 270 dashed line: x/δo∗ = 300; dotted line: x/δo∗ = 450. Steady case;k = 10;k = 1;
k = 0.2.

Figure 6. Phase-averaged streamwise velocity for the k = 10 case at the streamwise location x/δo∗ = 300. All key phases are displayed.
Black solid and dashed lines represent the extreme phases Φ = 270◦ and Φ = 90◦ respectively, corresponding to APG-FPG and FPGAPG. Colours represent the intermediate phases in the cycle. The solid line represent the phases between 270◦ and 90◦ (APG-ZPG andZPG FPG phases), and symbols represent the matching phases from 90◦ to 270◦ (FPG-ZPG and ZPG-APG phases).
Φ = 270◦ ;
◦
◦
◦
◦
◦
◦
◦
◦
◦
◦
Φ = 90 ;Φ = 288 , 252 ;Φ = 306 , 234 ;Φ = 324 , 216 ;Φ = 342 , 198 ;Φ = 0 , 180◦ ;
◦
◦
◦
◦
◦
◦
◦
◦
Φ = 18 , 162 ;Φ = 36 , 144 ;Φ = 54 , 126 ;Φ = 72 , 108 .

5

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

Figure 7. Phase-averaged streamwise velocity for the k = 1 case at the streamwise location x/δo∗ = 300. All key phases are displayed.
The same notation as in Figure 6 is used.

Figure 8. Phase-averaged streamwise velocity for the k = 0.2 case at the streamwise location x/δo∗ = 300. All key phases are displayed. The same notation as in Figure 6 is used.

6