Do coherent structures organize scalar mixing in a turbulent boundary layer?

Abstract

A scalar emanating from a point source in a turbulent boundary layer does not mix homogeneously, but is organized in large regions with little variation of the concentration: uniform concentration zones. We measure scalar concentration using laser-induced fluorescence and, simultaneously, the three-dimensional velocity field using tomographic particle image velocimetry in a water tunnel boundary layer. We identify uniform concentration zones using both a simple histogram technique, and more advanced cluster analysis. From the complete information on the turbulent velocity field, we compute two candidate velocity structures that may form the boundaries between two uniform concentration zones. One of these structures is related to the rate of point separation along Lagrangian trajectories and the other one involves the magnitude of strong shear in snapshots of the velocity field. Therefore, the first method allows for the history of the flow field to be monitored, while the second method only looks at a snapshot. The separation of fluid parcels in time was measured in two ways: the exponential growth of the separation as time progresses (related to finite-time Lyapunov exponents and unstable manifolds in the theory of dynamical systems), and the exponential growth as time moves backward (stable manifolds). Of these two, a correlation with the edges of uniform concentration zones was found for the past Lyapunov field but not with the time-forward future field. The magnitude of the correlation is comparable to that of the regions of strong shear in the instantaneous velocity field.

1. Introduction

A plume of pollutant that emanates from a point source in a turbulent flow spreads mainly through turbulent diffusion. Turbulence widens its boundary and mixes scalar concentration inside. Surprisingly, this inside mixing is not homogeneous. Inside the plume, large patches remain with little variation of the concentration, separated by sharp boundaries across which the concentration jumps. We call these regions ‘uniform concentration zones’ (UCZs). These regions are reminiscent of the regions of uniform momentum in turbulent boundary layers (Meinhart & Adrian 1995; Hutchins et al. 2012; Eisma et al. 2015). This has inspired our nomenclature, but we will highlight the connection with ramp–cliff structures below.

Recent work suggests that these large-scale motions strongly influence the dispersion of air pollution in urban areas which are immersed in the atmospheric boundary layer (Michioka & Sato 2012; Perret & Savory 2013). For instance, Michioka & Sato (2012) showed that instantaneous large-scale low momentum regions passing over idealized street canyons largely overlap with regions of high pollutant concentrations, while Perret & Savory (2013) related these large-scale motions with the coupling or decoupling of the flow within the canyon and the overlying flow.

The key question then is how zones of uniform concentration are related to structures in the turbulent velocity field. This leads to two problems addressed in this paper: the first problem is how to recognize UCZs, that is, how to devise an algorithm that reproduces zones that can also be recognized by the eye. The second problem is what structures of the velocity field make up the boundaries of UCZs.

UCZs, separated by edges, correspond to plateaus in the concentration field, separated by cliffs. There is an obvious connection with the ramp–cliff structure of the scalar field in shear flows (Warhaft 2000), and in isotropic turbulence with an imposed mean concentration gradient (Holzer & Siggia 1994; Tong & Warhaft 1994) and more recently by Iyer et al. (2018) and Buaria et al. (2021). As these cliffs also exist in isotropic turbulence, large-scale structures, as observed in shear flows, are not necessary to generate scalar cliffs. We use cluster analysis to search for scalar structures. It can find regions where the concentration is (approximately) uniform, but it cannot find ramp–cliff structures.

Ramp–cliff structures have been linked to converging–diverging separatrices of the instantaneous flow field (Warhaft 2000). Actually, this notion predates the concept of local maximum ridges of the finite-time Lyapunov exponent (FTLE) field which are currently thought to organize scalar concentration (Haller 2015), and which are candidate coherent structures in the present article. Lyapunov exponents gauge the exponentially fast spreading of fluid parcels that start close together. In three-dimensional (3-D) flow, the three Lyapunov exponents add up to zero because of incompressibility. The connection between the probability density functions (p.d.f.s) of the smallest Lyapunov exponent (corresponding to compression) and the p.d.f. of local concentration has been discussed by Götzfried et al. (2019), demonstrating the relevance of FTLE for the organization of concentration in sheets. The width of these sheets has been related to the distribution of FTLE (Kushnir, Schumacher & Brandt 2006). However, to the best of our knowledge, no direct connection between concentration sheets and maximum ridges of the FTLE field in 3-D turbulence has been made.

The ramp–cliff structure of concentration is related to intermittency, reflected in the scaling anomaly of high-order structure functions of scalar fluctuations (Iyer et al. 2018; Buaria et al. 2021). This anomalous scaling has been proven in the case of random velocity fields – a remarkable feat of turbulence theory (Shraiman & Siggia 2000).

We have done experiments on dispersion of a passive scalar from a point source in a turbulent boundary layer in a water channel. Information about the complete 3-D velocity field is obtained using tomographic particle image velocimetry (PIV). Simultaneously, the concentration field of the released scalar is registered using laser-induced fluorescence (LIF). Knowledge about the velocity field and its gradients allows the computation of coherent flow structures.

The average p.d.f. of local scalar concentration in turbulence is a smooth near-Gaussian function, but our interest is in snapshots of the concentration field with size comparable to the boundary layer thickness. The associated p.d.f.s are marred by local maxima that signify large regions with near-homogeneous concentration. Identifying these regions starts with these p.d.f.s. We discuss two methods, of which cluster analysis is the preferred one.

In this paper, we also explore the relevance of two different types of flow structures for the organization of the concentration field. The first one is related to the rate at which two close points separate in the velocity field. It involves an integration over a time interval $T$ using the velocity gradient field along a Lagrangian trajectory. A scalar field is constructed from the logarithm of the spreading rate: the FTLE $\varLambda _T(\boldsymbol {x})$. Spreading can also be tracked for backward time, then, the FTLE $\varLambda _{-T}(\boldsymbol {x})$ measures the rate of convergence of two fluid parcels towards $\boldsymbol {x}$. Ridges of local maxima of these scalar fields are known as Lagrangian coherent structures (LCS) (Shadden, Lekien & Marsden 2005). In 2-D time-dependent chaotic flows, LCS are analogous to separatrices, in time-in dependent chaotic fields, where they form impenetrable barriers for scalar transport. In time-dependent two-dimensional flows, LCS form transient barriers of transport, while they are advected almost passively. These objects offer the tantalizing prospect of predicting the spread of contaminants from knowledge of large-scale coherent structures only.

This concept has been applied in 2-D turbulence (Haller & Yuan 2000; Haller 2015) and effectively 2-D geophysical flows; see Bettencourt, Lopez & Hernandez-Garcia (2012) and references therein. Twardos et al. (2008) analysed an experiment on a weakly turbulent 2-D flow, using the related concept of stretching fields (Haller 2015). Line-shaped maxima of the finite-time past Lyapunov field were found to act as barriers of scalar transport. The finite-time Lyapunov field was measured in a turbulent boundary layer in a plane perpendicular to the boundary by Chong, Wang & Zhang (2009). The relatively small main velocity $U_\infty \approx 9\times 10^{-2}\ \textrm {m}\ \textrm {s}^{-1}$ and small ${Re}_\theta = 481$ allowed the tracing of structures up to 2 s, and allowed identification of hairpin-like structures in $\varLambda _{-T}(\boldsymbol {x})$. Wilson, Tutkun & Cal (2013) identify LCS in a plane parallel to the boundary from the 2-D velocity field measured in a turbulent boundary layer. By positioning the plane close to the boundary, the relevant integration time $T$ is the local eddy turnover time, which could be met in the Eulerian frame. The observation time was further extended using Taylor's frozen turbulence hypothesis.

LCS may be a useful concept in fully developed 3-D turbulence. When they organize as sheets, LCS structures may hinder scalar mixing and may thus be associated with boundaries of UCZs. In this paper, we will explore this concept in an experiment of a turbulent boundary layer that provides the full 3-D velocity field. Using this information, we measure the finite-time Lyapunov field in a plane perpendicular to the boundary. In addition, we identify the edges of UCZs, and perform a conditional average of the finite-time Lyapunov field on these edges. The other flow structures studied in this paper are regions of strong shear. These zones organize momentum transport, and correlate with the boundaries of uniform momentum zones (Eisma et al. 2015). The question is whether these regions are also correlated with the edges of the UCZs.

LCS became popular in the context of 2-D (geophysical) flows. The question is how to fruitfully apply this interesting concept to scalar dispersion in fully developed 3-D turbulence.

We will first describe our experimental methods in § 2. Next, we explain in § 3 how we extracted UCZs, while the computation of finite-time Lyapunov fields and shear regions are explained in §§ 4.1 and 4.2, respectively. With these preliminaries, we then describe in § 5 the relation between coherent structures of the velocity field and the edges of UCZs.

2. Experiment

In our experiment we inject dye from a point into a turbulent boundary layer in a water channel. The set-up is illustrated in figure 1. The channel has a 5 m long experiment section and cross-section of $0.6 \times 0.6\ \textrm {m}^2$. Fluorescent dye (Rhodamine B) is injected from a point at the boundary, 0.75 m upstream from the measurement volume, with injection velocity equal to the local mean velocity, $0.5 \, U_{\infty }$, where $U_\infty$ is the free-stream velocity. This ensured minimal perturbation of the velocity field. The measurement volume, located 3.5 m from the boundary layer tripping point, is at one of the sidewalls. The 3-D velocity field in a box $L_x \times L_y \times L_z = 5.72 \times 6.29 \times 0.57 \, \textrm {cm}^3$ was measured using tomographic PIV. In terms of the boundary layer height $\delta _{99} = 3.81\ \textrm {cm}$ the measurement volume is $1.5 \times 1.65 \times 0.15 \, \delta _{99}^3$. Simultaneously, the 3-D concentration field was measured in five slices spanning the $L_z = 0.15 \,\delta _{99}$ spanwise side of the measurement volume. Both fields were registered using a scanning laser sheet at $\lambda = 527\ \textrm {nm}$.

Figure 1. (a) Overview of the experimental set-up indicating the different components: four PCO Dimax PIV cameras to measure the velocity field using tomographic PIV, a Photron Fastcam LIF camera, a high-speed ND:YLF laser and the dye injection needle. The solid green line indicates the path of the laser beam/sheet. (b) Schematic view: cyl indicates a cylinder lens, focal length, $f = -25, 90\, \textrm {mm}$, sph is a spherical lens, $f = 1500\, \textrm {mm}$.

For the simultaneous measurement of the velocity and scalar concentration a scanning tomographic PIV system was combined with scanning LIF. Four high-speed CMOS cameras (Dimax, PCO) are used for tomographic PIV. Two cameras were equipped with $f=200$ mm objectives, whereas the others were equipped with $f=105$ mm objectives, all operating at $f_\# =11$. This optics involved Scheimpflug adapters to allow for large viewing angles while keeping the particle images in focus. The top two PIV cameras were positioned at an angle of $\beta =30^\circ$ with respect to the normal of the light sheet, whereas both side-viewing PIV cameras were looking at an angle of $\beta =45^\circ$. In order to minimize astigmatic aberrations, water-filled prisms were employed for all PIV cameras. LIF is performed using a single high-speed camera (Fastcam, Photron) equipped with a $f=105$ mm objective at $f_\# =8$.

Illumination was provided using a 50 mJ double pulsed Nd:YLF laser (Darwin Duo 80M, Quantronix). The laser beam was focused with a spherical lens ($f = 1500$ mm), and a $90$ mm wide and $1.4$ mm thick light sheet was subsequently formed by two cylindrical lenses with focal lengths of $-25$ and $90$ mm. This light sheet spans the streamwise ($x$) and the wall-normal direction ($y$) of the turbulent boundary layer, while the scanning is performed in the spanwise direction ($z$). The measurement volume is scanned with five thin laser sheets of $1.4$ mm thickness each, with an overlap of approximately 30 %. A scanning frequency $f_s$ of $640$ Hz results in a recording frequency of the full measurement volume of $128$ Hz.

The flow is seeded with $10\,\mathrm {\mu }$m diameter tracer particles (Sphericell) which are nearly neutrally buoyant. The seeding density in the tomographic particle image velocimetry (TPIV) images is around $0.025$ particles per pixel. The excitation of the used dye (Rhodamine B) matches the $527$ nm operating wavelength of the laser. Separation of the PIV and LIF signals is obtained by employing lowpass filters (PIV cameras) and a highpass filter (LIF camera).

Calibration of the distorted LIF and PIV images was performed using a 3-D calibration plate (Type 11, LaVision). The mapping of the distorted images was done using third-order polynomials in the $x$- and $y$-directions, whereas linear interpolation is used in the $z$-direction. Self-calibration, as proposed by Wieneke (2005), was performed to enhance the accuracy of the particle reconstruction. After several refinements, the remaining disparities were typically of the order of 0.1–0.2 pixels.

The particle volumes were reconstructed from the PIV images using five iterations of the fast-MART algorithm (Elsinga et al. 2006). The particle volumes were 3-D cross-correlated using an iterative multi-grid scheme reaching a final interrogation box size of $32^3$ voxels. This yields a spatial resolution of approximately 1.1 mm, corresponding to 34 wall units in each direction. From the relation for the Kolmogorov length $\eta$, $\eta ^+ = (\kappa y^+)^{1/4}$, with $\kappa = 0.41$ (Stanislas, Perret & Foucaut 2008) we estimate $\eta \approx 0.2\, \textrm {mm}$ at our largest distance from the wall. Our velocity resolution is $\approx 5\, \eta$; resolving $\eta$ is probably not possible in TPIV. However, our interest is in large-scale structures. Velocity gradients were obtained by employing a second-order spatial regression filter (Elsinga et al. 2010) with a filter size of 1.5 times the dimension of the interrogation box in each direction.

Care was taken to ensure that the concentration in the LIF images is sufficiently low for the fluid to be optically thin (Crimaldi 2008). In order to convert the local light intensity to a local dye concentration, a calibration procedure was performed. For this purpose, a small container with a known uniform concentration is positioned inside the measurement domain, after which images were recorded. This procedure was repeated for a few different known concentrations, resulting in a calibration curve for each pixel. Linear interpolation was employed to obtain the dye concentration from intensities at each pixel location. After calibration small negative concentrations were present in the free-stream region of the flow, which is a result of the finite accuracy of the LIF calibration. In order to remove these negative concentrations, the mean value of the concentration in the free-stream part of the flow is set to zero by subtracting it from the concentration field.

The turbulent characteristics of the boundary layer were measured with the TPIV set-up, with results shown in figure 2(a,b). Stereo PIV at a sampling rate of $40$ Hz and collecting a total of $6\times 10^3$ samples was used to check the results. Good agreement was found between the two methods. As figure 2 illustrates, the tomographic PIV results also agreed well with the data from DeGraaf & Eaton (2000).

Figure 2. (a) The mean velocity profile in inner scaling units. The dashed line shows the data from DeGraaf & Eaton (2000) at ${Re}_\theta = 2900$. The full line is the law of the wall, $U^+ = ({1}/{\kappa }) \ln y^+ + B$, with the von Kármán constant $\kappa = 0.41$ and offset $B = 5$. (b) The profiles of the root-mean-square fluctuations of the three velocity components, they have been normalized with the free-stream velocity $U_\infty = 0.74\ \textrm {m}\ \textrm {s}^{-1}$, and are also compared with data at ${Re}_\theta = 2900$ from DeGraaf & Eaton (2000) (dashed lines). (c) Mean concentration profile across the boundary layer. The maximum concentration $c_m =0.074\ \textrm {mg}\ \textrm {l}^{-1}$ is used for normalization of the concentration statistics.

At the measurement location the main stream velocity is $U_\infty = 0.74 \, \textrm {m}\ \textrm {s}^{-1}$, the boundary layer thickness is $\delta _{99}=38$ mm. The resulting momentum and displacement thicknesses are $\theta =4.0$ mm and $\delta ^* = 5.2$ mm, respectively, and the friction velocity is $u_\tau = 31 \, \textrm {mm}\ \textrm {s}^{-1}$. The momentum thickness Reynolds number is ${Re}_\theta = 3050$.

The mean scalar concentration $\bar {c}$ profile across the boundary layer is shown in figure 2(c). The maximum concentration is found around $y/\delta _{99}=0.1$. The concentration approaches zero in the free-stream part of the flow, as expected. The maximum value of the mean concentration profile $c_m=0.074\ \textrm {mg}\ \textrm {l}^{-1}$.

The 3-D velocity field in the $L_x \times L_y \times L_z = 1.5 \times 1.65 \times 0.15 \; \delta _{99}^3$ measurement volume was stitched together from the TPIV data in each of the 5 $z$-positions of the 1.4 mm thick laser sheet. The velocity field was corrected for the advection of the mean flow in the time interval ($\approx 1.5\ \textrm {ms}$) between two subsequent laser sheet positions. The quality of the measured velocity field was assessed by computing its divergence. The true velocity field has zero divergence; for the experimental data we compute the normalized divergence computed over the entire velocity field

(2.1)\begin{equation} \zeta = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} / \left[ \boldsymbol{\nabla} \boldsymbol{u} : (\boldsymbol{\nabla} \boldsymbol{u})^{\dagger} \right]^{1/2}, \end{equation}

which has zero mean and mean square variation $\langle \zeta ^2 \rangle ^{1/2} = 0.21$. The value of this quality measure is comparable to that in Casey, Sakakibara & Thoroddsen (2013), and references therein.

The $z$-resolution of the LIF images is determined by the width of the scanned laser sheet. In the present paper we focus on the concentration field in the central $z = 0$ plane. Therefore, this concentration field is an average over the 1.4 mm wide laser sheet. Much as for the scanned velocity field, also the concentration field has been corrected for the streamwise advection by the mean flow during a scan cycle.

Summarizing, the outcome of our experiment is a completely resolved 3-D velocity field together with a concentration field that is approximately resolved in slices of the 3-D measurement volume. This information is illustrated in figure 3. In the following sections we will quantify the relation between coherent structures of the concentration field and those of the velocity field.

Figure 3. An illustration of the information gathered in the present experiment. Iso-surfaces of the shear vorticity (defined in (4.5) in § 4.2) are shown in green, with $\omega _{sh} = 2.8\: U_\infty / \delta _{99}$ whereas vortical structures are show by the blue iso-surfaces with constant $Q = 0.95 \, U_\infty ^2 / \delta _{99}^2$. The concentration field is shown together with the in-plane velocity vectors viewed in a frame of reference convecting at $U_c = 0.9\, U_\infty$.

3. Uniform concentration zones

Figure 4 shows three snapshots of the concentration field, taken one large-eddy turnover time $\delta _{99} / U_\infty$ apart. In these snapshots, UCZs can be recognized as large regions where the variation of the concentration is small, bordered by sharp gradients. The concentration fields are from a 1.4 mm thick slice, centred in the spanwise extent of the 0.57 cm thick measurement volume.

Figure 4. Snapshots of the concentration field in the central $z = 0$ plane of the measurement volume. The snapshots are taken one large-eddy turnover time $\delta _{99} / U_\infty$ apart.

The idea of UCZs is intuitively appealing, but methods to find them must be well defined. Here, we follow the well-established methods of cluster analysis (Bezdek 1980), which were recently introduced in turbulence (Fan et al. 2019). We will compare our results with the method used by de Silva, Hutchins & Marusic (2016), Adrian, Meinhart & Tomkins (2000) and Eisma et al. (2015).

The starting point is the histogram $P_{his}(c)$ of concentrations. In the case of $n_{zones}$ clustered concentration values (which may still be spread over the observation plane), this histogram is characterized by $n_{zones}$ well-separated peaks. Each cluster corresponds to a peaked distribution. The concentrations in each of those still vary, but the variation is much smaller than that in the entire snapshot. In practice, clusters are not so well defined and the art becomes estimating the number $n_{zones}$ of clusters and to decompose the histogram into a sum of peaked distributions that are centred on $n_{zones}$ distinct concentration values. The first task is done using the statistical technique of kernel density estimation (Silverman 1986), while the tessellation of the concentration field in zones is done using the fuzzy cluster method (Bezdek 1980).

The histogram $P_{his}(c)$ is a discrete estimation of the underlying p.d.f. $P(c)$. The number of clusters is the number of local maxima of $P(c)$. A standard statistical procedure exists to overcome ambiguity associated with the width and placement of the discrete histogram bins (Silverman 1986). The trick is to endow each pixel value with a Gaussian distribution of concentrations with width $h$, and sum them. A question is the choice of the smoothing factor $h$. It can be proven that, for a Gaussian probability distribution function $P(c)$, the optimum $h$ is $h \approx \sigma _P\: n^{-1/5}$, where $\sigma _P$ is the standard deviation of $P(c)$. Of course, in our case $P(c)$ is not a Gaussian (we are counting its local maxima), but experience teaches us that this choice for $h$ is adequate.

Cluster analysis to find UCZs can be compared with the ad hoc method for finding uniform momentum zones described by de Silva et al. (2016), Adrian et al. (2000) and Eisma et al. (2015). Let us now briefly describe this method. For each concentration snapshot $c(\boldsymbol {x}, t)$ a histogram $P_{his}(c_i)$ is constructed of $n_{his} = 32$ equally spaced concentrations $c_i, i = 1, \ldots , n_{his}$ that span the dynamic range of $c(\boldsymbol {x}, t)$. The number $n_{his}$ is chosen small enough to counter noise, but large enough to distinguish peaks. Then, (i) the maximum $c_i^{max}$ of $P_{his}$ is found. (ii) This maximum is represented by a Gaussian $P_G(c_i) = P_c(c_i^{max}) \exp (-(c_i - c_i^{max})^2 / \sigma ^2)$, while its width $\sigma$ is determined in a least squares procedure. (iii) This peak is deleted from the histogram: $P_{his}(c_i) \rightarrow \textrm {max}(P_{his}(c_i) - P_G(c_i), 0)$, after which steps (i) through (iii) are repeated $n_{zones}$ times. (iv) After sorting the Gaussians with respect to their peak position, the intersections of neighbouring Gaussians are found. These intersections constitute the boundaries of the UCZs. Notice that this assignment of zones differs from that in the cluster analysis: concentration in the tail of a Gaussian is now associated with another cluster. No such ambiguity exists in the cluster analysis, which is also insensitive to the choice of discrete histogram bins. The two methods described will find different zone edges, and will find a different association with flow structures.

The average number of zones from the kernel density estimate is $n_{zones} = 3.7 \pm 0.9$, with 0.9 the root-mean-square variation. Therefore, we fixed the number of zones to four. An estimate of the number of zones from the discrete histogram is fraught with uncertainty. We will always rank the zones according to their area. Thus, the ‘blue sky’ in figure 5 has rank 1, and the neighbouring (green) zone has rank 3.

Figure 5. (a) Snapshot of tracer concentration $c(\boldsymbol {x}, t)$ in the centre plane of the measurement volume. (b) The concentration field of (a), coarse grained into $n_{zones} = 4$ UCZs using the simple histogram method. (c) Same as (b), but using the cluster method. (d) Histogram method: the thick grey line is the histogram of concentrations $c$, with the range of $c$ in (a) binned into 32 concentration levels. The lines are four Gaussian concentration profiles. The grey dashed vertical lines, the intersections of adjacent Gaussians, indicate the boundaries of the uniform concentration zones in (b). The dashed line is the sum of the Gaussians. This method results in clear discrepancies between the measured and reconstructed histograms. (e) Same as (d), but now for the cluster method. The lines are the histograms of the four clusters. The dots are now the reconstructed histogram, they perfectly coincide with the measured histogram. In (b,d) the zones are numbered according to their area. (f) Kernel density estimation: the dots are the measured histogram, the full line is its reconstruction through Gaussian smoothing. The number of local maxima of the smooth curve is $n_{zones} = 4$; the maxima are indicated by the arrows.

The outcome of both methods for a selected concentration snapshot is illustrated in figure 5; both methods find approximately the same zones, but there are differences. Below we will argue that the cluster method is superior, as it can discriminate more acutely between flow structures.

Apart from the technique of kernel density estimation, we believe that there does not exist an objective way to estimate the number of uniform zones. Moreover, as has been argued by Fan et al. (2019) in the context of uniform momentum zones, adding or subtracting 1 from the number of zones did not change their identification of interfacial layers.

4. Coherent structures of the velocity field

4.1. Finite-time Lyapunov exponents

FTLEs gauge the exponentially fast spreading of nearby points. They were computed from the measured velocity field $\boldsymbol {u}(\boldsymbol {x}, t)$ and strain field $\boldsymbol{\mathsf{A}}(\boldsymbol {x}, t)$, ($\boldsymbol{\mathsf{A}} = \boldsymbol {\nabla } \boldsymbol {u}$, components $A_{i, j} = \partial u_i / \partial x_j$), by integrating the evolution of small separations $\boldsymbol {\delta }(t)$ along a Lagrangian trajectory,

(4.1)\begin{equation} \frac{\,\textrm{d} \boldsymbol{\delta}}{\,\textrm{d} t} = \boldsymbol{\mathsf{A}}(\boldsymbol{x}(t), t) \boldsymbol{\cdot} \boldsymbol{\delta}(t), \quad \text{with} \ \frac{\,\textrm{d}\kern0.7pt \boldsymbol{x}(t)}{\,\textrm{d} t} = \boldsymbol{u}(\boldsymbol{x}, t). \end{equation}

The time integration of (4.1) over a time interval $[t_0, t]$ defines the evolution matrix $\boldsymbol{\mathsf{M}}_{\,t_0}^{\,t}$ as $\boldsymbol {\delta }(t) = \boldsymbol{\mathsf{M}}_{\,t_0}^{\,t} \boldsymbol {\cdot } \boldsymbol {\delta }(t_0)$. This matrix was computed on a 2-D grid of initial points $\boldsymbol {x}_0 = \boldsymbol {x}(t_0)$ in the central $(z = 0)$ plane of the measurement volume. As time progresses, fluid parcels may leave this plane and explore other $z$ values; these parcels are tracked in our 3-D velocity field. Therefore, the full 3-D measured velocity field was used in the computation of $\varLambda _{\pm T}$. Choosing the central plane minimizes the loss of fluid parcels exiting in the spanwise direction. The resulting Lyapunov field at $z = 0$ is compared with the concentration field of the central slice; the information in the other four slices was not used. The Lyapunov exponent is related to the separation of infinitesimally close points. When a computation is done using actual points advected by the velocity field $\boldsymbol {u}(\boldsymbol {x}, t)$, their separation has to remain small, which necessitates their replacement when their separation grows too large. If not, the advected points which started near $\boldsymbol {x}_0$ may move to completely uncorrelated regions of the velocity field and no longer reflect the Lyapunov exponent for the Lagrangian trajectory that started at $\boldsymbol {x}_0$. This is irrelevant if instead the gradient field $\boldsymbol{\mathsf{A}}$ is used, but the statistical noise of $\boldsymbol{\mathsf{A}}$ is larger than that of $\boldsymbol {u}$. The accuracy of measuring gradients of the velocity field in our experiment has been addressed through the divergence error (2.1). The integration of (4.1) was done using a simple forward Euler scheme, and the fields at $\boldsymbol {x}(t)$ were computed from the measured 3-D PIV velocity fields using trilinear interpolation.

The largest eigenvalue $\lambda _3$ of the positive Cauchy–Green tensor

(4.2)\begin{equation} \boldsymbol{\mathsf{C}}_{\,t_0}^{\,t} = \boldsymbol{\mathsf{M}}_{\,t_0}^{\,t} \: \left(\boldsymbol{\mathsf{M}}_{\,t_0}^{\,t} \right)^{\dagger}, \end{equation}

with ${\dagger}$ the adjoint, defines the finite-time future Lyapunov exponent $\varLambda _T$ as

(4.3)\begin{equation} \varLambda_{t_0}^t(\boldsymbol{x}_0) = \frac{1}{2 T} \ln(\lambda_3), \end{equation}

with $T = t - t_0$. Our measurements of the velocity field are Eulerian, which means that the time interval $T$ is restricted by the advection of fluid parcels until they exit the field of view. Consequently, there is a distribution of Lagrangian times $T$ over the field of view. Most of the time, trajectories leave the downstream side of the field of view, but they may also exit a $z$-boundary of the measurement volume. In figure 6 we show these times for one snapshot of $\varLambda _{-T}$. (Throughout we use the shorthand $\varLambda _{T}$ for $\varLambda _{t_0}^{t_0+T}$ and $\varLambda _{-T}$ for $\varLambda _{t_0}^{t_0-T}$.) Backward in time moves us upstream, so that the integration times are shortest for the left part of the frame. With increasing integration times structures of $\varLambda _{\pm T}$ become narrower.

Figure 6. (a) Snapshot of the past Lyapunov field $\varLambda _{-T}$, which is obtained from the measured velocity field and its gradient field by integrating (4.1) backward in time, $t: 0 \rightarrow -T$. (b) Because of the mean velocity in the turbulent boundary layer, there is a distribution of the maximum integration times $T$ available for the computation in (a); it is shown in (b). Since the mean flow is from left to right, the integration time is longest for the rightmost part of the snapshot shown in (a). In the lower right corner, the integration time is limited by fluid parcels exiting through a lateral ($z$) boundary of the measurement volume.

The Lyapunov field $\varLambda _T(\boldsymbol {x})$ measures the spreading of nearby fluid parcels in the future. However, the boundaries of UCZs were formed in the past. The past of fluid structures can be studied by starting Lagrangian trajectories at $\boldsymbol {x}_0$ in the field of view, and integrating equation (4.1) backward in time. The corresponding Lyapunov field $\varLambda _{t_0}^{t_0-T}(\boldsymbol {x})$ is then again defined in terms of the eigenvalue $\lambda _3$.

(4.4)\begin{equation} \varLambda_{{-}T}(\boldsymbol{x}_0) ={-}\frac{1}{2 T} \ln(\lambda_3). \end{equation}

Note that the smallest eigenvalue of $\boldsymbol{\mathsf{C}}_{\,t_0}^{\,t_0-T}$ is comparable to the forward-time Lyapunov exponent at time $t_0-T$. Local maxima of this field indicate regions in the flow where points were separated in the past, but have converged to $\boldsymbol {x}_0$ at the instant of observation. On the other hand, structures are advected, so that future Lyapunov exponents already existed in the past and then may have organized the concentration distribution. The question then is whether it is the past or future Lyapunov field that is most relevant for boundaries of the UCZs. Another question concerns the time $t'$ at which the scalar field $c(\boldsymbol {x}, t')$ should be compared with structures of $\boldsymbol{\mathsf{C}}_{\,t}^{\,t+T}$ or $\boldsymbol{\mathsf{C}}_{\,t}^{\,t-T}$: should $t'$ be chosen $t$ or a later or earlier time?

When computing statistics, we restrict field coordinates $(x, y)$ to the upstream half of the observation window for past Lyapunov exponents, and the downstream half of the future Lyapunov exponents. This corresponds to integration times $|T| \gtrsim 0.05\; \textrm {s}$. In finding the correlation of $\varLambda _{\pm T}$ with the edges of UCZs we only include points where $\varLambda _{\pm T}$ is convex in the direction of the eigenvector $\boldsymbol {\xi }_3$ corresponding to the largest eigenvalue $\lambda _3$ of $\boldsymbol{\mathsf{C}}_{t}^{t\pm T}$, $\widetilde {\boldsymbol {\xi }}_3\boldsymbol {\cdot } \boldsymbol{\mathsf{H}} \boldsymbol {\cdot } \widetilde {\boldsymbol {\xi }}_3 < 0$, with the Hessian $\boldsymbol{\mathsf{H}}_{\,i j} = \partial ^2 \varLambda _{\pm T}(\boldsymbol {x}) / \partial x_i \partial x_j$, and $\widetilde {\boldsymbol {\xi }}_3$ the projection of $\boldsymbol {\xi }_3$ on the $x y$-plane. More steps to refine the FTLE field to ridges were mentioned by Farazmand & Haller (2012). We finally impose a threshold $\varLambda _{\pm T} > 4 \ \textrm {s}^{-1}$. This threshold corresponds to the noise level of $\varLambda _{-T}$ in the ‘blue sky’ outside the boundary layer (zone 1 in figure 5).

4.2. Shear vorticity $\omega _{sh}$

As a second velocity field candidate for the edges of UCZs, we consider regions of strong shear. In order to isolate those, the 2-D $(u, v)$ projection of the velocity field was decomposed into three parts (Kolár 2007). This decomposition was introduced originally to provide a more robust description of vortices in turbulent flows, but is now used to detect shear layers (Maciel, Robitaille & Rahgozar 2012).

With this aim, the velocity gradient $\boldsymbol {\nabla } \boldsymbol {u}$ is separated into components representing rigid-body rotation, elongation and the desired true shear. The shear vorticity $\omega _{sh}$ is then defined as

(4.5) \begin{align} \omega_{sh} = \left| \frac{\partial \tilde{u}}{\partial \tilde{y}} {\,-\,} \frac{\partial \tilde{v}}{\partial \tilde{x}} {\,-\,} \left\{ \mathrm{sgn} \left(\frac{\partial \tilde{u}}{\partial \tilde{y}}\right) \cdot \mathrm{min} \left( \left|\frac{\partial \tilde{u}}{\partial \tilde{y}}\right|, \left|\frac{\partial \tilde{v}}{\partial \tilde{x}}\right| \right) {\,-\,} \mathrm{sgn} \left(\frac{\partial \tilde{v}}{\partial \tilde{x}}\right) \cdot \mathrm{min} \left( \left|\frac{\partial \tilde{u}}{\partial \tilde{y}}\right|, \left|\frac{\partial \tilde{v}}{\partial \tilde{x}}\right|\right) \right\} \right|, \end{align}

where $\widetilde {\boldsymbol {u}}$ and $\widetilde {\boldsymbol {x}}$ refer to a coordinate system formed by the principal axes of $\boldsymbol {\nabla } \boldsymbol {u}$, rotated over $45^\circ$. Local maxima of this field are the shear layers found by Meinhart & Adrian (1995), which separate regions with nearly uniform velocity (Eisma et al. 2015). A generalization of (4.5) does not exist, but an alternative way to find shear layers, but now using the full 3-D velocity field, has been proposed by Horiuti & Takagi (2005). We also would like to point to Haller & Beron-Vera (2012), who propose quantifying shear using the eigenvectors of the Cauchy–Green tensor (4.2).

Summarizing, our measurement of the space–time velocity field in a 3-D volume allows us to compare three candidate velocity structures. The past Lyapunov field $\varLambda _{-T}(\boldsymbol {x})$, which relates to the past organization of the concentration field and attracting structures, the 2-D shear vorticity $\omega _{sh}(\boldsymbol {x})$, which refers to the instantaneous velocity field and the future Lyapunov field $\varLambda _{T}(\boldsymbol {x})$ which corresponds to repelling structures.

5. Results

Figure 7 shows the concentration, the finite-time Lyapunov exponent, and the shear vorticity in a snapshot of the $x y$ plane centred in the measurement volume. We register time series at a frequency of 128 Hz ($\approx 6 \: U_\infty / \delta _{99}$), lasting $6.5 \, \textrm {s}$, corresponding to $1.2\times 10^2$ large-eddy turnover times $\delta _{99} / U_\infty$. In collecting statistics, a narrow region near the wall, where gradients are large, was excluded. This region is indicated in figure 7; its width is $150\;y^+$ (4.73 mm), which includes the peak in $u'$ in figure 2. Our experiments are in the Eulerian reference frame. Due to limitations on the integration time, the field $\varLambda _{-T}$ in the left side of figure 7(b) is more diffuse than in the right side of the figure.

Figure 7. A snapshot from the measured time series of scalar concentration and 3-D velocity field. (a) Concentration, coarse grained into four UCZs, (b) past Lyapunov field $\varLambda _{-T}$, (c) shear vorticity $\omega _{sh}$. In collecting statistics, a narrow region near the wall (indicated by the grey region) was excluded. The field $\varLambda _{-T}$ has been filtered to show only regions with negative curvature in the direction of the eigenvector $\boldsymbol {\xi }_3$. For clarity, only values $\varLambda _{-T}, \omega _{sh} > 20$ are shown.

We will now discuss the relation between the $\varLambda _{-T}$ and $\omega _{sh}$ fields, and finally present our main result, that is the relation between maxima of the finite-time Lyapunov field and edges of UCZs.

5.1. Statistics of $\varLambda _{\pm T}$ and $\omega _{sh}$

Loosely speaking, $\varLambda _{\pm T}$ is the time average of the largest stretching induced by the strain $\boldsymbol{\mathsf{A}}$ over a Lagrangian trajectory, while $\omega _{sh}$ is the instantaneous value of the strongest shear that it induces. The p.d.f.s of $\varLambda _{-T}$ and $\varLambda _{T}$ in figure 8(a) are remarkably similar, but that of $\omega _{sh}$ has a long tail. Local maxima of all three fields are candidate structures organizing the distribution of passive scalar. Still, their nature is very different. Not surprising, the fields $\varLambda _{-T}$ and $\omega _{sh}$ are not strongly correlated, as is illustrated in figure 8: for each $\omega _{sh}$, $\varLambda _{-T}$ has a broad distribution, with the average $\varLambda _{-T}$ increasing very slightly with increasing $\omega _{sh}$.

Figure 8. (a) Full lines: probability density function of $\varLambda _{-T}$ (black), $\varLambda _{T}$ (grey), and $\omega _{sh}$ (dashed). (b) Joint histogram of $\varLambda _{-T}$ and $\omega _{sh}$; the colour scale of the histogram is logarithmic. The full line indicates $\langle \varLambda _{-T} \rangle$ as a function of $\omega _{sh}$.

The regions of large $\omega _{sh}$ in figure 7 are diffuse blobs, so are those of $\varLambda _{-T}$. Sharp line-like structures of $\varLambda _{\pm T}$ are only expected for long enough integration times $T$, which should be of the order of the large-scale turnover time. On average, our maximum integration time $T \approx 0.1\, \textrm {s}$, while the typical eddy turnover time $\delta _{99} / u_{rms} \approx 1\, \textrm {s}$. Clearly, sharply defined structure of $\varLambda _{\pm T}$ needs a true Lagrangian measurement. Nevertheless, it is possible to measure the statistics of blobs with large $\varLambda _{\pm T}$ and $\omega _{sh}$. By computing the second-order moment matrix of these regions and finding the eigenvector with the largest eigenvalue, the orientation of these structures can be found. As figure 9 illustrates, patches of large $\varLambda _{\pm T} \ge 20\ \textrm {s}^{-1}$ are inclined at an angle of $\varphi \approx 20^\circ$, which agrees with Chong et al. (2009) who identified these regions with hairpin vortices. Patches of large $\varLambda _{T}$ do not have a preferential orientation whereas the shear field $\omega _{sh}$ appears to be oriented parallel to the boundary. As will be argued below, edges of concentration zones correspond to $\varLambda _{-T}$, and much less to $\varLambda _{T}$ and the association of $\varLambda _{-T}$ (and not $\varLambda _{T}$) with genuine structures is tempting.

Figure 9. The orientation $\varphi$ of (unconnected) regions $\varLambda _{\pm T} \ge 20\ \textrm {s}^{-1}$ and $\omega _{sh} \ge 20\ \textrm {s}^{-1}$. (a) The angle $\varphi$ is defined as the angle of the longest principal axis of the region with the horizontal. (b,c) The p.d.f.s of $\varphi$. Regions of $\omega _{sh}$ are mostly oriented parallel to the boundary; those of $\varLambda _{-T}$ are inclined at an angle of $\varphi \approx 20^\circ$.

5.2. Conditional averages

In finding a connection between the edges of uniform scalar concentration zones and regions of large $\varLambda _{\pm T}$ we follow the procedure of Bisset, Hunt & Rogers (2002), but with an important caveat. The procedure is to compute the conditionally averaged function $\langle \varLambda _{\pm T}(x, y-y_0; t) \rangle _{x, t}$, with $y_0$ the vertical (wall-normal) coordinate of an edge of a uniform concentration zone. The caveat is that we compute conditional averages relative to those where $y_0$ is chosen randomly. Specifically, we show the ratio of the average function $\langle \varLambda _{\pm T}(x, y-y_0; t) \rangle _{x, t}$, with $y_0$ the true zone edge over the average function with $y_0$ taken randomly. When this ratio is one, no distinction can be made between the true and random $y_0$, and zone edges are not correlated with $\varLambda _{\pm T}$. The random choice of $y_0$ is such that its p.d.f. is the same as that of the true $y_0$. We will denote the relative averages as $\widetilde {\varLambda }_{-T, T}$ and $\widetilde {\omega }_{sh}$. Further details are in the Appendix.

Conditional averages are shown in figure 10. We have computed relative conditionally averaged vertical profiles of $\varLambda _{\pm T}$ and $\omega _{sh}$ on the boundaries of the second largest UCZ shown in figure 10(d), which corresponds the green zone in figures 5(b,c) and 7(a), but the conclusions are comparable to those involving all boundaries. In addition, we have removed small isolated patches (relative area $< 5\times 10^{-3}$ of the full frame) of this zone and zones of other rank. Finally, we recall that we only consider $\varLambda _{-T}$ for the right half of the frame, and vice versa for $\varLambda _{T}$.

The central result of this paper, shown in figure 10, involves the conditional statistics of $\varLambda _{\pm T}$ and $\omega _{sh}$ on zone edges. Those edges are determined in two ways: using the cluster method (figure 10b,c,e,f,g) and using the simple histogram method (figure 10d). In the case of the cluster method, the correlation of the past Lyapunov field $\varLambda _{-T}$ with concentration edges is much larger than that of the future $\varLambda _{T}$ (figure 10b,c). This agrees with the result of Twardos et al. (2008) in weakly turbulent 2-D flow who found that contours of an advected scalar field follow the past stretching field, and that of Kushnir et al. (2006) in numerical turbulence (${Re}_\lambda = 24$) who conclude that gradient sheets of the dissipation field are formed by the contracting eigenvalue of $\boldsymbol{\mathsf{M}}_{\,T}$. A similar correlation with edges of the concentration zones is seen for the relative shear vorticity $\widetilde {\omega }_{sh}$ in figure 10(e).

Figure 10. Conditional averages of the finite-time Lyapunov fields $\langle \varLambda _{\pm T}(x, y-y_0; t) \rangle _{x, t}$, and shear vorticity $\langle \omega _{sh}(x, y-y_0; t) \rangle _{x, t}$ on the vertical (wall-normal) coordinates $y_0$ of UCZs, relative to those conditioned on randomly picked coordinates $y_0$. (a) Illustration of conditional averaging on the edge of the zones with ranks 2 and 3 (zone 3 is the green zone in figures 5b,c and 7a). For the past Lyapunov field $\varLambda _{-T}$, averages are restricted to the right half of a frame, since fluid parcels there have the longest observation history. The left half (grey region) is used for $\varLambda _{T}$. (b,c) Relative conditional averages $\widetilde {\varLambda }_{-T, T}$ as a function of the time delay $\Delta t$ of the concentration field; $\Delta t$ increases from $\Delta t = -0.023\ \textrm {s}$ to $\Delta t = 0 \, \textrm {s}$. The edges of the concentration field were found with the cluster method. (d) Same as (b,c) but for $\Delta t = 0$, and edges found with the histogram method. (e) Relative shear vorticity $\widetilde {\omega }_{sh}$, conditionally averaged on concentration edges found with the cluster method. (f) Averages $\widetilde {\varLambda }_{-T}$ for $\Delta t = 0$, with zone edges found from the cluster method, with $n_{zones} = 3, 4$ and 5. (g) Values of $\widetilde {\varLambda }_{-T}$ and $\widetilde {\omega }_{sh}$ conditionally averaged on edges of uniform momentum zones, found with the cluster method. Notice the change of the vertical scale.

The increase at a zone edge of the past FTLE field $\varLambda _{-T}$ is a mere 0.05 (5 %) above the random background $\widetilde {\varLambda }_{-T} = 1$, and a similar amount above the future $\widetilde {\varLambda }_{T}$. The significance of this enhancement should be compared with the concentration of $\widetilde {\omega }_{sh}$ on edges of uniform momentum zones (regions of uniform $u(\boldsymbol {x}, t)$), $\widetilde {\omega }_{sh} \approx 1.3$ (30 %), shown in figure 10(g). It was derived from the present data using the same cluster analysis. The linkage of these two structures – both derived from the velocity field – has been discussed in Eisma et al. (2015), Christensen & Adrian (2001) and Elsinga & Marusic (2010).

The significance of our result is further illustrated in figure 10(b,c) where the correlation is seen to depend sensitively on the time delay $\Delta t$ between the FTLE field and the concentration zone edges. Specifically, we correlate edges of the concentration field $c(\boldsymbol {x}, t + \Delta t)$ with the FTLE fields $\widetilde {\varLambda }_t^{t + T}$ and $\widetilde {\varLambda }_t^{t - T}$ for $\Delta t$ increasing from $-0.023\, \textrm {s}$ to 0 s. Although the field $\widetilde {\varLambda }_t^{t - T}$ explores the past, the conditional average peaks at the present.

The cluster method needs the number of zones $n_{zones}$ as input. Throughout, we have fixed it to the averaged kernel density estimate, $n_{zones} = 4$, both for the cluster and the histogram method. Figure 10(f) shows the conditional average $\widetilde {\varLambda }_t^{-T}$ for $n_{zones} = 3, 4$ and 5. There is no significant effect of setting $n_{zones} = 5$, but there is for $n_{zones} = 3$. We always count the blue sky as a zone, and take conditional averages on the edges of the zone with rank 2. Therefore, fixing $n_{zones} = 3$ almost always amounts to averaging with respect to the turbulent–non-turbulent interface.

For the edges found with the histogram method (figure 10d) there is almost no distinction between past and future Lyapunov exponents; perhaps illustrating the superiority of the cluster method.

Summarizing, we find a small but significant correlation between the edges of UCZs, and both the past-time Lyapunov field and the shear vorticity field.

6. Conclusion

We have quantified the large-scale organization of scalar concentration in a turbulent boundary layer. In analogy with uniform momentum zones in turbulent boundary layers (Meinhart & Adrian 1995; Hutchins et al. 2012; Eisma et al. 2015) we have used the term ‘uniform concentration zones’, but we emphasize the strong connection with ramp–cliff structures (Warhaft 2000; Iyer et al. 2018; Buaria et al. 2021). Two different procedures were described for the identification of UCZs, with superior results being obtained by the cluster method. Still, it can miss zones, or can find zones which are not obvious to the eye. After all, the concentration inside a zone is not strictly uniform. Ramp–cliff structures share the sudden change of concentration (cliffs) with the edges of UCZs. However, cluster analysis cannot identify ramp–cliff structures.

Next, we considered structures of the velocity field as candidate zone edges. From the measured 3-D velocity field and its gradients we compute the intersection of the 3-D finite-time Lyapunov field $\varLambda _{\pm T}$ with the central plane of the measurement volume. Due to the advection by the mean velocity, the longest time interval over which trajectories could be integrated is $T \approx 2\: \delta _{99} / U_\infty$. This may be too short to see a clear development of a Lagrangian coherent structure. Further refinement of the $\varLambda _{-T}$ field to its local maximum ridges may decrease the random background of the measured correlations.

Concentration edges found with the cluster method only correlate with the past finite-time Lyapunov field $\varLambda _{-T}$. They select converging rather than diverging flow structures. No such distinction between past and future can be made for concentration edges found with the histogram method. Perhaps this points to the superiority of the cluster method.

The concentration of both $\varLambda _{-T}$ and $\omega _{sh}$ fields on the edges of UCZs is quite small compared with the random background. Using our techniques, it is a mere 5 %. This value can be compared with the concentration of shear vorticity on edges of uniform momentum zones, which is well established in the literature; it is 30 % using the same techniques.

Our conclusion is that the question as to which large-scale velocity structures organize the scalar concentration is still open.