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

BISPECTRAL ANALYSIS OF ATTACHED EDDY MODELING OF
WALL-BOUNDED TURBULENCE

Oles Dubrovski
Faculty of Mathematics
Technion Israel Institute of Technology
Haifa, 32000 Israel

Guangyao Cui
Faculty of Aerospace Engineering
Technion Israel Institute of Technology
Haifa, 32000 Israel

Ian Jacobi
Faculty of Aerospace Engineering
Technion Israel Institute of Technology
Haifa, 32000 Israel
ijacobi@technion.ac.il

ABSTRACT

stream of a turbulent boundary layer, in a pattern reminiscentof the streamwise skewness profile. The correlation between
the different scales has been variously interpreted as the resultof amplitude or frequency modulation of the envelope of small
scale motions by large scales, or a spatial phase delay between
the scales.
Duvvuri & McKeon (2015) performed a formal analysis of the scale correlation coefficient which showed that despite the somewhat artificial filtering of the velocity signal intolarge- and small-scale components, in reality the correlation
coefficient represented interactions across three triadicallycoupled scales. This observation is a direct consequence of theconvective non-linearity of the incompressible Navier Stokes
equations when written in spectral form, which allows only
non-linear interactions between triadically coupled wavenumbers k0 , k00 , and k = k0 + k00 . In light of this, Cui & Jacobi (2021) proposed that a non-linear diagnostic tool which
accounts for the three distinct wavenumbers would be better
suited to the scale-interaction analysis than the traditional twocomponent, cross-correlation.

Bispectral analysis was applied to the attached eddymodel of wall-turbulence in order to examine whether the linear superposition of attached eddies can exhibit the spectral
signature of non-linear, triadic scale interactions. Traditional
analysis of scale interactions in turbulence has focused on correlations between filtered velocity signals of large- and smallscale motions. In contrast, the bispectrum provides a natural
measure of non-linear scale interactions without filtering. The
magnitude of the bispectrum represents the strength of the nonlinear coupling, while the phase of the bispectrum represents adelay within the triadic interaction. In this study we utilized anumerical implementation of the attached eddy model to generate an ensemble of synthetic velocity fields, and we compared the averaged bispectrum of these fields to a DNS channel flow. The linear superposition of discrete hairpin packets
was found to produce a similar bicoherence signature to thatof velocity fields generated via true non-linear interactions, although the nature of the wavenumber coupling differed significantly. The biphase behavior of the attached eddy model
also differed from the DNS. Tentative interpretations of bothof these differences are offered in the context of the linear superposition mechanism, itself.

Bispectral Analysis
The type of quadratic, non-linear interactions generated
through turbulent convection are naturally described by thebispectrum, B(k0 , k00 ), which is the third-order spectrum of asignal, u(x), defined as:

BACKGROUND
Scale-Interactions In Wall-bounded Turbulence
The interactions between large- and small-scale featuresof wall-bounded turbulence represent the primary mechanism
for energy transfer in turbulence and thus a key target for developing turbulence models and efficient flow control strategies. Bandyopadhyay & Hussain (1984) first proposed a crosscorrelation diagnostic to measure the relationship between thefiltered large- and small-scale signals of a turbulent shear flow,
and showed that the two scales are related in a way that varies
across the width of the shear layer. This approach was later
refined and extended by Mathis et al. (2009) and others whoreported that large- and small-scale motions were positively
correlated near the wall and inversely correlated near the free-B(k0 , k00 ) = ub(k0 )bu(k00 )bu∗ (k0 + k00 ) .

(1)

Just as the sum of the traditional, second-order power spectrum
is related to the statistical variance, hu2 i, the bispectrum is related to the statistical skewness hu3 i by hu3 i = ∑k0 ,k00 B(k0 , k00 ).
However, unlike the real-valued power spectrum, the complex bispectrum preserves phase information between triadically interacting modes in addition to the amplitude information which describes the magnitude of the triadic coupling.
Recently, the bispectrum has received increased interest due to its preservation of triadic interaction information.

1

(a)

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

(c)

(c)

Dimension
Height
Leg Gaussian Vorticity Std.
Leg Effective Radius
Streamwise Spacing
Packet Length
Vertex Angle between Legs
Hairpin Inclination Angle
Packet Inclination Angle
Number of Radial Fibers
Numer of Azimuthal Fibers
Distance Between Radial Fibers
Azimuthal Distance Between Fibers
Fiber Longitudinal Segment Length

β

(d)

(d)

Inner Units
H ≈ 128σ = 0.05H ≈ 6.4
R ≈ 3σ ≈ 19.2
b = 0.4H ≈ 50
L p = 6b ≈ 300
45◦
θ = 80◦α = 10◦nr = 20
nt = 40dr = R/nr ≈ 0.96R(2π/nt ) ≈ 3≈ 5.7

Table 1: Top: The dimensions for the self-similar hairpin geometry corresponding to the smallest hairpin packet, definedby its height, H. Packets are sampled uniformly from the family of packets with heights {1, 2, 4, 8, 16}H, where all length
dimensions in the table are scaled proportionally. Bottom: The
dimensions of the discrete fibers used to represent the hairpin
vorticity distribution for the Biot-Savart calculation.

Figure 1: (a) The streamwise, normalized bispectrum (bicoherence), b(k00 , k0 ), of a turbulent channel flow DNS evaluatedat y+ ≈ 281. Note the high level of bicoherence associatedwith the LSMs circled in the red ellipse. (b) The corresponding
biphase map, β (k00 , k0 ), marked with two different phase locations corresponding to an interaction between (c) two LSMsand (d) an LSM with a smaller scale motion. The quadraticphase-coupled (QPC) interaction and phase delay for the two
scales are illustrated schematically in panels (c) and (d) (corresponding to the circular marks), where the phase delay, β ,is illustrated spatially but in reality represents a more complicated function of the interacting modes.

Perry & Chong (1982) to describe higher order statistics aswell, via spectral analysis. Perry (1987) also applied the attached eddy model numerically by assuming a particular formof the hairpin vortex, calculating the induced velocity field,
and distributing these induced velocities randomly in space according to a reciprocal distribution with respect to eddy size.
de Silva et al. (2016) and Eich et al. (2020) have since refined
this numerical approach to better align the predicted velocity
fields with experiments.
In the present study, we implemented a version of thenumerical attached eddy model, following the hairpin packet
specifications of de Silva et al. (2016). Each individual hairpin was composed of a bundle of discrete, parabolic vortex
fibers, with a Gaussian vorticity distribution with standard deviation, σ across the bundle. The fiber bundles were arranged
to construct the individual hairpins and then a packet of multiple hairpins with maximum height, H = 128 (inner units), as
shown in figure 2. The geometric parameters describing theindividual hairpins and the discrete fibers are detailed in table
1.
The Biot-Savart law was then applied to the vorticity field
described by the fibers in order to generate the induced velocity field from the entire packet. Because the packets weregeometrically self-similar, they could be scaled up in factorsof 2n H to obtain larger packets, where n represents the number of generations of hairpins. The maximum multiple considered here, Hmax = 24 H = 2048 (limited by computer memory
resources) represents the half-height of the simulated channel
flow field and thus also corresponds to its friction Reynolds
number, Reτ .
The velocity fields of the hairpin packets were then uniformly, randomly selected from the geometric set of packet
sizes, and the packet velocity fields were superposed in a uniform, random process to construct an individual volumetricflow snapshot of a channel flow. The channel domain extended
24 half-channel heights in the streamwise direction and 4 in the
spanwise direction, with uniform spatial resolution, ∆x+ ≈ 20,
in all directions. The number of packets applied per snapshot
was fixed to obtain a well-converged, logarithmic mean velocity profile, where the magnitude of the fiber circulation was

Schmidt (2020) defined a bispectral modal decomposition for
low-order representation of turbulent flows. Cui & Jacobi(2021) utilized the normalized amplitude (bicoherence, b) andphase (biphase, β ) of the bispectrum to study the scale interactions in a turbulent channel flow and found that the dominantcoupling between scales was associated with the large-scale
motions (LSMs) of the channel flow interacting with a rangeof other scales, from other LSMs down to smaller scales, as
shown in figure 1(a). They also showed that the biphase between the triadic interactions varied with scale size and could
be interpreted as an interaction delay associated with the scaleinteractions, as illustrated in figure 1(b). However, the precise
meaning of this interaction delay was not fully explored. In
particular, does the biphase represent an interaction delay due
to spatial proximity of the different scale velocity modes, oris it a result of temporal proximity due to variations in theirrespective convection velocities? And is the interaction delay
associated with the particular geometry of the velocity modeshapes? These questions have important implications for modeling efforts involving the superposition of large-scale velocity modes to reconstruct turbulent flow fields. In order to address these questions, we apply the bispectral analysis to velocity fields constructed from a physically-meaningful modelof large-scale motions – the attached eddy model, which involves only linear, spatial-superposition.

METHOD
Simulation via Attached Eddy Model
Townsend’s attached eddy model presents a fundamentally physical approach to explaining the statistical propertiesof wall-bounded turbulence by constructing the turbulent velocity field via the superposition of hierarchies of self-similar
eddies. Although first applied analytically to understanding
the logarithmic mean velocity profile, it was later extended by2

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

x+

kφuu
u2τ

(b)

kφuu
u2τ

α
100

H
y+ 50
θ
0
50

z+

0
-50
-200

-150

-100

-50

0

50

100

150

200

b
Lp
Figure 4: Pre-multiplied, streamwise energy spectra for (a)the attached eddy simulation and (b) the channel DNS. The
wavenumber, k, is normalized by the channel half height.

Figure 2: The smallest hairpin packet, H ≈ 128, illustrating
the discrete vortex fibers used to calculate the packet-induced
velocity via Biot-Savart.

(a)u0 /uτ

(b)

u0 /uτ

lation has a peak around k ≈ 1 and is thus associated with aregion of wavenumbers lower than those corresponding to thepacket sizes themselves (k ≈ 3–43), indicating the streamwisealignment of individual packets to generate large-scale, meandering motions. The early implementations of the numericalattached eddy simulation by Perry (1987) also found a significant energy shift to lower wavenumbers, in comparison toexperiments, although they found a peak energy around k ≈ 8.
The precise shape of the energy spectrum for the attached eddy
simulation is largely determined by the shape of the eddy itself,
as noted by Marusic & Perry (1995); in the present study we
examined only a parabolic, λ -shaped hairpin.

RESULTS
Bicoherence of the Attached Eddy Model

Figure 3: Wall-parallel maps of the fluctuating streamwise velocity, u0 , at y+ ≈ 176 for (a) the attached√eddy simulation andy+ ≈ 281 (b) the channel DNS. y+ = 3.9 Reτ corresponds tothe middle of the log layer for both cases.

Having verified the spectral energetic structure of the attached eddy model, the streamwise bispectrum was calculated
according to (1) and then normalized by the energy spectral
density associated with the different wavenumber triads to obtain the bicoherence, b(k0 , k00 ). Because the attached eddymodel is composed of mostly LSMs, a significant amount of its
spectral energy resides at wavenumbers near the wavenumberof the finite computational domain. This can result in significant spectral leakage, which becomes particularly prominent
in the bispectral calculation, and thus special care was takenwith the removal of the DC component of the velocity signals,
as described in the Appendix. The bicoherence map evaluatedat the middle of the log layer is shown in figure 5.
The bicoherence map for the attached eddy simulationdiffers significantly from that of the channel flow DNS at thesame wall-normal location, shown above in figure 1(a). In
the channel flow, the bicoherence indicated a continuous bandof strong coupling between the LSMs and the full range ofother scales in the flow, which was interpreted in the context of large-scale modulation of the small scales. For the attached eddy simulation, the bicoherence indicates a discrete
set of strong triadic couplings between scales that differ inwavenumber by an order of magnitude, such that k00 ≈ 10 k0 .
The fact that the most prominent bicoherence is constrained
along this line indicates that large-scales interact most strongly
with scales that are substantially smaller, by an order of magnitude, and not with scales at neighboring wavenumbers. Perhaps this is partly explained by the relative lack of small-scale
structures in the attached eddy construction.
The bicoherence across the k00 ≈ 10 k0 slice is shown infigure 6, along with circle-symbols for wavenumbers corresponding to the streamwise extent (L p ) for the different discrete packets. The three largest packets are associated withused as an arbitrary fitting parameter to obtain the desired vonKarman slope (1/κ), where all packets were assumed to advect at the free-stream velocity.
An example of the resulting instantaneous, wall-parallel
velocity field is shown in figure 3(a), with comparison to
the corresponding field obtained from a channel flow DNS at
Reτ = 5200 by Lee & Moser (2015), shown in figure 3(b). The
simulated field exhibits meandering, low- and high-speed velocity regions characteristic of experimentally observed flows,
but clearly composed of much larger, more coherent scales,similar to a low-pass filtering of the DNS field.
The streamwise energy spectra for the simulated velocity
field were calculated by averaging over N = 480 snapshots until statistical convergence was obtained. A map of the premultiplied energy spectra is shown in figure 4 for (a) the attachededdy simulation and (b) the corresponding DNS. We note firstthat the simulated velocity field is limited by the spatial resolution near the wall and thus begins at y+ = 20. However,
the attached eddy model is designed to primarily represent thecoherent structures inhabiting the log-layer, and thus this limitation is not significant. Secondly, we observe that the DNS
spectral map shows the expected two regions of intense spectral energy, associated with the small scale motions involved in
the near wall cycle at y+ ≈ 15, and LSMs inhabiting the outerpeak, at the center of the log layer. The attached eddy simulation captures only this second spectral signature of LSMs in
the center of the log layer.
Moreover, the spectral energy in the attached eddy simu-

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

Figure 6: The bicoherence slice along k00 = 10 k0 in the middleof the log layer taken from figure 5 is shown in the black line.
The sharp peaks appear at wavenumbers in multiples of 2. The
wavenumbers corresponding to the streamwise extent (L p ) of
the individual hairpin packets, k0 = {2.7, 5.4, 10.7, 21.4, 42.9},
are marked with red circles. The lowest-wavenumber, sharp
peak in the bicoherence corresponds to k0 = 1.4, which is twice
the size of the largest hairpin packet. The broader peak in bicoherence at the lowest wavenumbers is the result of spectral
leakage as discussed in the Appendix.

Figure 5: Bicoherence map in the middle of log layer (y+ =1/2
3.9Reτ ≈ 176) for the attached eddy simulation at Reτ =2048. The red ellipse marks the region of high bicoherence associated with the wavenumber triad defined by k0 and
k00 ≈ 10k0 , denoted by the black dashed line.

prominent peaks in bicoherence, whereas the smaller packets
are not. There is also a sharp bicoherence peak at twice thesize of the largest packet. The question about the interpretationof this high bicoherence region thus remains: why the largest
packets don’t seem to interact with the next largest packets inthe packet hierarchy, i.e. why packets separated by a single
generation don’t show significant bicoherence and exhibit ahigh bicoherence along k00 ≈ 2 k0 ?
The precise peak locations in the bispectrum likely reflect the spatial location of the packets in the flow; modifying
the placement of the packets from the entirely random process used here to a constrained random process that controls
the minimum distance between similar packets, as discussed
in de Silva et al. (2016), is expected to alter the observed peakstructure of the bicoherence. Thus, the lack of observed interactions between neighboring packets may be attributed totheir relative spatial orientation in the flow and perhaps due toa lack of overlap, although this question remains unresolved.
The convection velocity of the individual packets is also expected to influence the interaction behavior captured by thebicoherence.
Finally, it is worth noting that the bicoherence is constructed to represent the triadic coupling between scales,
which results from a non-linear interaction. In the attachededdy model, however, the non-linear interaction is simulated
via a purely linear superposition of scales, and nevertheless thebicoherence detects the scale-interactions as if they were due
to non-linear coupling. This provides direct evidence that thelinear superposition process used in attached-eddy type models can be shown to emulate even the non-linear coupling associated with the turbulent convection.

energy to a third. Indeed, we showed previously that the averaged biphase is directly related to the direction of the turbulent
energy cascade. In the case of the DNS channel flow, shown
above in 1(b), the biphase was shown to vary with wavenumber, where nearly all triads exhibited positive interaction delays (β > 0, i.e. all shades of red in the bipshase map), consistent with the classical energy cascade. However, as notedin the introduction, the precise interpretation of the interaction
delays was unclear: were they associated with the geometry
and spatial distance between the interacting modes or a temporal delay in the interaction dynamics? The attached eddymodel removes the non-linear dynamics because the velocity
field is constructed entirely by means of linear superposition of
hairpin packets. Therefore, for the attached eddy biphase calculation, the delay interpretation is more naturally associatedwith spatial distance between scales and the physical geometry of the interacting modes. The biphase map evaluated at the
middle of the log layer is shown in figure 7.
Unlike the DNS biphase which was almost uniformly positive, the attached eddy biphase exhibits both positive and negative delays, depending on the wavenumber triad. Following
along the k00 = 10 k0 slice of significant bicoherence, we seethat the biphase is positive, as in the DNS. But, if instead, we
followed along the k0 = 1 line of the LSM interactions, thenthere is a transition between positive and negative biphase.
Large-scale interactions with neighboring large scales show apositive delay, but large-scale interactions with small-scales
show a negative delay. This suggests that if the high bicoherence region could be shifted or extended by means of modification of the superposition algorithm, then the biphase may
also follow suite, and thus the interaction delay can be tunedby the spatial orientation and placement of the packets.

Biphase of the Attached Eddy Model
The biphase can be interpreted as a measure of the interaction delay between two velocity modes interacting to transfer

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

(a)

(b)

(c)

(d)

Figure 8: Premultiplied streamwise energy spectra (a,c) andbicoherence calculations (b,d) for the DNS (top) and attachededdy (bottom) where the red line indicates unweighted DC
component removal and the blue line indicates the windowweighted DC component removal. Note the significant reduction in DC leakage for the attached eddy bicoherence (d) whenusing the window-weighted mean, in contrast to the minimal
effect in all other calculations. The bicoherence in (b) and (d)is measured along k0 = 2; all spectral measurements are at the
middle of the log layer.

Figure 7: The streamwise biphase map in the middle of log
layer (y+ = 176). The black dashed line is the same line as infigure 5.

CONCLUSIONS
Bispectral calculations were performed for a modest
Reynolds number attached eddy simulation, with five generations of hairpin packets, and contrasted with a DNS channel
flow. The bicoherence showed evidence for non-linear coupling between large- and small-scale motions, even though
the simulation itself involved only linear superposition of hairpin packets, thus providing evidence that the spectral featuresof non-linear dynamics can be modeled using linear models
alone. However, the nature of the scale-interactions for the attached eddy model was much narrower in wavenumber-spacethan that of the DNS, and indicated that only a limited set ofwavenumber triads participated. The sliced bicoherence wasused to identify the specific hairpin packets participating inthese interactions. The biphase, which has been shown torepresent an interaction delay between different wavenumber
modes, and also the direction of average energy transfer due
to the interactions, was shown to exhibit the same sign for theattached eddy model as for the DNS, at least in the region ofhigh bicoherence. Overall, these results suggest that modification of the linear superposition mechanism in the attachededdy model is the key to altering the bispectral signature of thescale-interactions, and thus the key to allowing the attachededdy model to more realistically capture true, non-linear dynamics in wall-bounded flows.

has zero mean. When this happens, the windowing process canactually introduce spurious energy at low frequencies. In order
to avoid this spectral leakage from the DC and near-DC components, Van der Schaaf & Van Hateren (1996) suggest thata ‘window-weighted’ average be removed from the instantaneous signal instead of the traditional, unweighted averaged.
The window-weighted average, hz(t)iw for a given windowfunction, w(t), is defined ashz(t)iw = hz(t)w(t)i/hw(t)i(2)

Once this average is removed, then the resulting signal can betapered and the DFT can be calculated following the standard
procedure.
The difference between the standard DC removal andthe weighted DC removal is shown in figure 8 for boththe DNS data, with predominantly higher wavenumber components, and the attached eddy data, with predominantly
lower wavenumber components. The greatest leakage without weighting is observed in the bicoherence for the attachededdy, and that is the case in which the weighting produces themost significant improvement. The choice of DC removal has
little impact on the other calculations.

APPENDIX
The typical procedure for removing the DC component ofa random signal, z(t), as described by Bendat & Piersol (1986),
prescribes subtracting the mean, hz(t)i from the instantaneous
signal, and then tapering the resulting data with an appropriate window, after which the discrete Fourier transform (DFT)
is calculated, with zero-padding. However, this procedure canintroduce significant artifacts to the resulting energy spectrum(or bispectrum) when the original signal is dominated by very
low frequency (wavenumber) components, in which case mostof the signal may be offset from zero even if, arithmetically, it

REFERENCES
Bandyopadhyay, P. R. & Hussain, A. K. M. F. 1984 The
coupling between scales in shear flows. Physics of Fluids
27 (9), 2221–2228.
Bendat, J.S. & Piersol, A.G. 1986 Random Data: Analysis andMeasurement Procedures, 4th edn. Wiley.
Cui, G. & Jacobi, I. 2021 Biphase as a Diagnostic for Scale
Interactions in Wall-Bounded Turbulence. Physical Review
Fluids 6 (1), 014604.

5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
boundary layers. Journal of Fluid Mechanics 628, 311–337.
Perry, A. E. 1987 Turbulence Modelling Using Coherent
Structures in in Wakes , Plane Mixing Layers and Wall Turbulence. In Perspectives in Turbulence Studies (ed. H. U.
Meier). Springer-Verlag Berlin.
Perry, A. E. & Chong, M. S. 1982 On the mechanism of wall
turbulence. Journal of Fluid Mechanics 119, 173–217.
Schmidt, O. T. 2020 Bispectral mode decomposition of nonlinear flows. Nonlinear Dynamics 102, 2479–2501.
de Silva, C. M., Woodcock, J. D., Hutchins, N. & Marusic, I.
2016 Influence of spatial exclusion on the statistical behavior of attached eddies. Physical Review Fluids 1 (2).
Van der Schaaf, A. & Van Hateren, J. H. 1996 Modelling thepower spectra of natural images: Statistics and information.
Vision Research 36 (17), 2759–2770.

Duvvuri, S. & McKeon, B. J. 2015 Triadic scale interactions
in a turbulent boundary layer. Journal of Fluid Mechanics
767, 1–11.
Eich, F., De Silva, C. M., Marusic, I. & Kähler, C. J. 2020
Towards an improved spatial representation of a boundary
layer from the attached eddy model. Physical Review Fluids
5 (3), 1–13.
Lee, M. & Moser, R. D. 2015 Direct numerical simulation of
turbulent channel flow up to Reτ ≈ 5200. Journal of Fluid
Mechanics 774, 395–415.
Marusic, I. & Perry, A. E. 1995 A wall-wake model for the turbulence structure of boundary layers. Part 2. Further experimental support. Journal of Fluid Mechanics 298, 389–407.
Mathis, R., Hutchins, N. & Marusic, I. 2009 Large-scale amplitude modulation of the small-scale structures in turbulent

6