2.1 Direct numerical simulation (DNS)

This study focuses on plane Poiseuille flow. Figure 1 shows the geometry of the simulation domain. A constant streamwise ( $x$ -direction) pressure gradient drives the flow between two infinite parallel plates. The periodic boundary condition is applied in the streamwise and spanwise ( $z$ -direction) directions with the period dimensions represented by $L_{x}$ and $L_{z}$ . A no-slip boundary condition is applied to the walls in the $y$ -direction (wall-normal). By default, non-dimensionalization using turbulent outer scales is applied to all variables: i.e. the half-channel height $l$ is used for the scaling of length, the laminar centre-line velocity $U_{c}$ for velocity, $l/U_{c}$ for time and $\unicode[STIX]{x1D70C}U_{c}^{2}$ for pressure (where $\unicode[STIX]{x1D70C}$ denotes the density of fluid). The Reynolds number is thus defined as $Re\equiv \unicode[STIX]{x1D70C}U_{c}l/\unicode[STIX]{x1D702}$ , where $\unicode[STIX]{x1D702}$ is the viscosity of the fluid. Turbulent inner scales are used to report results of near-wall flow statistics and structure, for which the friction velocity $u_{\unicode[STIX]{x1D70F}}\equiv \sqrt{\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}}$ and viscous length (or wall unit) $\unicode[STIX]{x1D6FF}_{v}\equiv \unicode[STIX]{x1D702}/\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F}}$ are used. Quantities so scaled are denoted with a superscript ‘ $+$ ’. Under these definitions, the friction Reynolds number, $Re_{\unicode[STIX]{x1D70F}}\equiv \unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F}}l/\unicode[STIX]{x1D702}$ , can be directly related to $Re$ through $Re_{\unicode[STIX]{x1D70F}}=\sqrt{2Re}$ . The governing equations of momentum and mass balances are

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D66B}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D66B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D66B}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D66B}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D66B}=0. & \displaystyle\end{eqnarray}$$

A Fourier ( $x$ )-Chebyshev ( $y$ )-Fourier ( $z$ ) pseudo-spectral scheme is adopted for spatial discretization while a third-order semi-implicit backward-differentiation Adams–Bashforth scheme (Peyret Reference Peyret2002) is used for time integration. DNSs have been performed at three different $Re$ , i.e.  $3600$ ( $Re_{\unicode[STIX]{x1D70F}}=84.85$ ), $14\,400$ ( $Re_{\unicode[STIX]{x1D70F}}=169.71$ ) and $80\,000$ ( $Re_{\unicode[STIX]{x1D70F}}=400$ ). A summary of the numerical settings for the DNSs of statistical turbulence is provided in table 1. The simulation domain is kept the same in inner units ( $L_{x}^{+}\times L_{z}^{+}$ ; and thus in outer units both $L_{x}$ and $L_{z}$ scale with $1/Re_{\unicode[STIX]{x1D70F}}$ ). Likewise, the grid sizes in the transverse directions $\unicode[STIX]{x1D6FF}_{x}^{+}$ and $\unicode[STIX]{x1D6FF}_{z}^{+}$ are also kept constant in inner units. The number of grid points in the $y$ -direction increases with $Re$ to keep the wall-normal resolution approximately the same in inner units. The numerical solver is implemented in a custom code parallelized with MPI based on the open source ChannelFlow package (Gibson Reference Gibson2012); the code was first reported in Tuckerman et al. (Reference Tuckerman, Kreilos, Schrobsdorff, Schneider and Gibson2014).

Table 1. Summary of the numerical settings for the DNS of statistical turbulence.

2.2 Streak transient growth (STG) simulation

In statistical turbulence, vortices are often irregular in shape, highly concentrated in space and intricately positioned relative to (sometimes partially connected with) one another. Meanwhile, for the initial test of our vortex-tracking algorithm, a benchmark system that enables controllable generation of well-defined three-dimensional vortex structures is required. We adopt the streak transient growth (STG) approach of Schoppa & Hussain (Reference Schoppa and Hussain2002) for this purpose, which controls the vortex configuration by adjusting several parameters of the initial condition. (As another option, one may as well follow the approach of Brandt & de Lange (Reference Brandt and de Lange2008) in which vortices of different configurations are generated from different modes of streak interaction.)

The initial condition for STG is constructed by superposing a base flow with a perturbation velocity. The base flow

(2.3a,b ) $$\begin{eqnarray}U_{b}(y,z)=U_{m}(y)+U_{s}(z)g(y),\quad V_{b}=W_{b}=0\end{eqnarray}$$

is quasi-two-dimensional ( $U_{b}$ , $V_{b}$ and $W_{b}$ are the $x$ -, $y$ - and $z$ -component, respectively) and itself a superposition of the mean velocity profile of statistical turbulence at the same $Re$

(2.4) $$\begin{eqnarray}U_{m}(y)\equiv \int _{0}^{\infty }\int _{0}^{L_{x}}\int _{0}^{L_{z}}v_{x}(x,y,z,t)\,\text{d}z\,\text{d}x\,\text{d}t,\end{eqnarray}$$

(where $v_{x}$ is the instantaneous streamwise velocity component in the statistical turbulence) with a streamwise velocity streak adjustment $U_{s}(z)g(y)$ . The latter is factorized into spanwise and wall-normal dependence terms

(2.5) $$\begin{eqnarray}U_{s}(z)=A_{s}\cos (\unicode[STIX]{x1D6FD}_{s}(z-z_{\unicode[STIX]{x1D6FD}}))\end{eqnarray}$$

and

(2.6) $$\begin{eqnarray}g(y)=y\exp (-\unicode[STIX]{x1D702}y^{2}).\end{eqnarray}$$

Here, $A_{s}$ adjusts the amplitude of the spanwise undulation, $\unicode[STIX]{x1D6FD}_{s}$ adjusts the spanwise streak spacing, $z_{\unicode[STIX]{x1D6FD}}$ is the spanwise phase parameter which is set so that the low-speed streak is aligned to the middle of the domain and $\unicode[STIX]{x1D702}$ is set to align the wall-normal maximum at $y^{+}=20$ . The perturbation velocity

(2.7) $$\begin{eqnarray}v_{x}^{\prime }=v_{y}^{\prime }=0,v_{z}^{\prime }=A_{p}\sin (\unicode[STIX]{x1D6FC}_{p}x)g(y)\end{eqnarray}$$

( $v_{x}^{\prime }$ , $v_{y}^{\prime }$ and $v_{z}^{\prime }$ are the $x$ -, $y$ - and $z$ -component, respectively) adds streamwise dependence to the base flow, without which the instability would not grow (Waleffe Reference Waleffe1997). Here, $A_{p}$ is the perturbation amplitude and $\unicode[STIX]{x1D6FC}_{p}$ is the streamwise wavenumber.

STG parameters used in this study for transient vortex generation, along with the numerical settings of the STG simulations, are listed in table 2. Note that $v_{z,rms}^{\prime +}$ is the root mean square (r.m.s.) magnitude of the spanwise perturbation velocity. A small simulation domain close to a minimal flow unit (MFU) (Jiménez & Moin Reference Jiménez and Moin1991) is used because we only need to focus on a small set of vortex structures for algorithm testing purpose. Vortices are generated by STG only in half of the channel: i.e. both the streak velocity and perturbation velocity are only applied at the $y<0$ side of the domain while for $y>0$ , $g(y)=0$ and the initial velocity is simply $U_{m}(y)$ .

Table 2. Numerical settings and initial condition parameters used for STG simulations.

4.1 Test of VATIP with STG-generated vortices

We start by testing the effectiveness of VATIP in STG flow fields, where vortex generation is controllable by the parameters of initial disturbance (Schoppa & Hussain Reference Schoppa and Hussain2002). Figure 4 shows the time series of the root mean square of $Q$ in our STG simulation (numerical settings given in § 2.2) and vortex configurations of selected moments are shown in figure 5. The initial disturbance flow field ( $t=0$ ) contains strictly streamwise vortices with a spanwise phase shifts between upstream and downstream vortex sets. At the beginning of STG, the $Q_{rms}$ profile starts to grow and reaches the first plateau at around $t=15$ . At this stage, the quasi-streamwise vortices tilt and bend sideways to the spanwise direction but wall-normal lifting up remains small ( $t=20$ in figure 5). After the first plateau, the $Q_{rms}$ profile continuously increases and reaches its peak at $t=75$ . During this period, neighbouring tilted-streamwise vortices lift up and conjoin to form well-defined hairpin vortices ( $t=60$ in figure 5). The value of $Q_{rms}$ gradually decreases after $t=75$ . Vortices in this period have a high lifting tendency despite their lower strength ( $t=120$ in figure 5).

Figure 4. Time series of the root mean square of $Q$ in the STG simulation. Moments of the flow fields shown in figure 5 are marked with red circles.

Figure 5. Vortex configurations of selected moments in the STG simulation. The isosurfaces are $Q=0.015$ for $t=0,20$ and $60$ and $Q=0.01$ for $t=120$ . The colour scale maps to the distance from the wall in outer units. Vortex axis lines from (a) a streamwise-only tracking approach (equivalent to the Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997) method) and (b) VATIP are compared (circular marks; different colours are used for different vortices as identified by the method).

Typical vortex configurations in these moments, including the strictly streamwise ( $t=0$ ), titled-streamwise ( $t=20$ ), lifted-up hairpin ( $t=60$ ) and decaying hairpin ( $t=120$ ) vortices, are used as our benchmark systems for vortex tracking. The original method of Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997) is recovered when the algorithm of figure 3(a) is truncated right after subroutine 1 (i.e. no iterative propagation in other directions). This would be sufficient if the target was limited to streamwise ( $t=0$ ) or quasi-streamwise vortices (as in the case of Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997)). Its inadequacy starts to surface in titled-streamwise vortices ( $t=20$ ) where the spanwise segment of the vortex is not fully captured in the axis line obtained (circular markers). For hairpin vortices ( $t=60$ and $120$ ) not only are some of the axis points missing (because they are not $yz$ -plane maxima of $Q$ ), the method also breaks the axis line of a well-defined hairpin into separate pieces. The new VATIP algorithm successfully identified the complete axis lines of vortices of all shapes and correctly grouped axis points of the same vortex into one axis line.

4.2 DNS: flow statistics and visualization

We now give an overview of the DNS results of statistical turbulence at three different $Re$ ( $Re_{\unicode[STIX]{x1D70F}}=84.85$ , $169.71$ and $400$ ) in this section. Application of VATIP to these flow fields will be discussed in § 4.3. The mean velocity $U^{+}$ as a function of $y^{+}$ is plotted in figure 6(a). As $Re$ increases, the profile outside the buffer layer ( $y^{+}\geqslant 30$ (Pope Reference Pope2000)) gradually approaches the von Kármán log law (Kim, Moin & Moser Reference Kim, Moin and Moser1987; Pope Reference Pope2000)

(4.1) $$\begin{eqnarray}U^{+}=2.5\ln y^{+}+5.5.\end{eqnarray}$$

At the lowest $Re_{\unicode[STIX]{x1D70F}}=84.85$ , the profile is slightly higher than the von Kármán asymptote, indicating that the log-law layer is not fully developed. The agreement is much better at the two higher $Re$ and at the highest $Re_{\unicode[STIX]{x1D70F}}=400$ , it nearly completely collapses onto (4.1) for a wide range of $y^{+}$ (until the channel centre).

Figure 6. (a) Mean velocity profiles ( $U^{+}$ versus $y^{+}$ ) and (b) log-law indicator functions ( $y^{+}\text{d}U^{+}/\text{d}y^{+}$ versus $y^{+}$ ) of the statistical turbulence at $Re_{\unicode[STIX]{x1D70F}}=84.85$ , $169.71$ and $400$ .

From a generic logarithmic profile

(4.2) $$\begin{eqnarray}U^{+}=A\ln y^{+}+B,\end{eqnarray}$$

the log-law slope can be expressed as

(4.3) $$\begin{eqnarray}A=y^{+}\frac{\text{d}U^{+}}{\text{d}y^{+}}.\end{eqnarray}$$

When the profile does not strictly follow a logarithmic dependence (4.2), $A$ becomes a function of $y^{+}$ – its variation indicates the departure from the log law. This quantity (4.3), which is thus sometimes referred to as the log-law indicator or diagnostic function (Hoyas & Jiménez Reference Hoyas and Jiménez2006; Marusic et al. Reference Marusic, McKeon, Monkewitz, Nagib, Smits and Sreenivasan2010), is plotted in figure 6(b) for our DNS results. For the lowest $Re_{\unicode[STIX]{x1D70F}}=84.85$ , the function goes nearly straight down with no discernible flat region, indicating the lack of a well-defined log-law layer (despite the fact that the profile is seemingly parallel to the von Kármán asymptote in figure 6 a). For the two higher $Re$ ( $Re_{\unicode[STIX]{x1D70F}}=169.71$ and $400$ ), an inflection point shows up at $y^{+}\approx 50$ with nearly the same value of $2.5$ , which agrees well with the von Kármán log-law slope reported in Kim et al. (Reference Kim, Moin and Moser1987) and is also consistent with the observations of Moser, Kim & Mansour (Reference Moser, Kim and Mansour1999) and Jiménez & Moser (Reference Jiménez and Moser2007). After the inflection point, the profile is not strictly flat but its variation is small for a distinct range of $y^{+}$ ( $50\lesssim y^{+}\lesssim 100$ for $Re_{\unicode[STIX]{x1D70F}}=169.71$ and $50\lesssim y^{+}\lesssim 320$ for $Re_{\unicode[STIX]{x1D70F}}=400$ ), indicating that these $Re$ are sufficiently close to and have already shared some common features with fully developed turbulence (Moser et al. Reference Moser, Kim and Mansour1999; Hoyas & Jiménez Reference Hoyas and Jiménez2006).

Figure 7. Reynolds stress profiles for $Re_{\unicode[STIX]{x1D70F}}=84.85$ , $169.71$ and $400$ .

Figure 7 shows the four components of Reynolds stress, $\langle v_{x}^{\prime +}v_{x}^{\prime +}\rangle$ , $\langle v_{y}^{\prime +}v_{y}^{\prime +}\rangle$ , $\langle v_{z}^{\prime +}v_{z}^{\prime +}\rangle$ and $-\langle v_{x}^{\prime +}v_{y}^{\prime +}\rangle$ , as functions of $y^{+}$ . Consistent with the literature (Moser et al. Reference Moser, Kim and Mansour1999; Abe, Kawamura & Matsuo Reference Abe, Kawamura and Matsuo2001), the profiles of all components rise with $Re$ at $y^{+}$ above ${\approx}30$ while the peak shifts towards the centre of the channel. The $Re$ -dependence is stronger in the transverse components $\langle v_{y}^{\prime +}v_{y}^{\prime +}\rangle$ and $\langle v_{z}^{\prime +}v_{z}^{\prime +}\rangle$ , which reflects increasing energy redistribution (Abe et al. Reference Abe, Kawamura and Matsuo2001), and the dependence in the streamwise component $\langle v_{x}^{\prime +}v_{x}^{\prime +}\rangle$ is much weaker.

Figure 8. Instantaneous vortex structures at (a) $Re_{\unicode[STIX]{x1D70F}}=84.85$ and (b) $Re_{\unicode[STIX]{x1D70F}}=400$ . Isosurfaces are identified by the $Q$ -criterion and in the wall-normal direction only the bottom half and 20 % of the top half of the channel are shown (i.e.  $0<y^{+}<1.2Re_{\unicode[STIX]{x1D70F}}$ ). The colour shade (from light to dark) maps to the distance from the bottom wall in outer units.

Figure 8 shows the vortex structures identified by the $Q$ -criterion in typical snapshots at the lowest ( $Re_{\unicode[STIX]{x1D70F}}=84.85$ ) and the highest ( $Re_{\unicode[STIX]{x1D70F}}=400$ ) Reynolds numbers. In both cases, the flow fields are filled with tube-like vortices. Quasi-streamwise vortices are more prevalent in the vortex field, but hairpin vortices can still be observed. Examples of these hairpins are shown in the enlarged views. Vortices at the higher $Re$ display a high extent of lifting up and many instances of detached vortices are observed. (A vortex becomes detached when its upstream legs leave the wall and become shielded from wall interaction (Perry & Marušić Reference Perry and Marušić1995; Marusic et al. Reference Marusic, McKeon, Monkewitz, Nagib, Smits and Sreenivasan2010).) Meanwhile, most of the vortices at the lower $Re$ remain attached to the wall with a comparatively weaker extent of lifting.

4.3 VATIP application in DNS: vortex classification and conformation

Visualization based on the $Q$ field can only provide a cursory glance of the instantaneous vortex fields and lacks both quantitative precision and statistical certainty. The new VATIP algorithm automatically detects vortices with a variety of shapes without subjective bias. It thus offers a feasible pathway to the statistical analysis of the population and configurations of vortex structures, facilitating the understanding of their roles in turbulent dynamics. In this section, VATIP is applied to the DNS results of statistical turbulence. Another algorithm is also proposed to classify vortices according to the topology of the vortex axis lines identified thereby. Note that a lower $Q$ threshold of $0.4Q_{rms}$ is used for VATIP as discussed previously in § 3.

Figure 9. Measurements of vortices in $x$ and $z$ dimensions.

To begin with, vortex size is measured according to figure 9: the streamwise and spanwise measurements ( $l_{x}^{+}$ and $l_{z}^{+}$ ) are defined as the maximal separation between axis points in these two dimensions, respectively. The statistical distributions of these measurements are presented in figures 10 and 11. (As discussed below, vortices with $l_{x}^{+}<50$ are considered fragments and not included in the statistics.) At all $Re$ , the probability density function (PDF) of $l_{x}^{+}$ monotonically decreases with increasing $l_{x}^{+}$ . The average $l_{x}^{+}$ is approximately $120$ and is nearly independent of $Re$ . This value is comparable with Jeong et al.’s (Reference Jeong, Hussain, Schoppa and Kim1997) $200$ (using their streamwise tracking algorithm) and Panton (Reference Panton2001)’s $100$ (from empirical observation). Sensitivity of this measurement to varying $Q_{threshold}\equiv HQ_{rms}$ is rather small: e.g. for $Re_{\unicode[STIX]{x1D70F}}=84.85$ , increasing $H$ from $0.4$ to $1.6$ (well beyond the percolation level), the average $l_{x}^{+}$ only decreases from $126$ to $104$ . This is consistent with the earlier (§ 3) statement that vortex axis topology is insensitive to the changing $H$ value.

Figure 10. Probability density function of the streamwise measurement of vortices: (a) $Re_{\unicode[STIX]{x1D70F}}=84.85$ , (b) $Re_{\unicode[STIX]{x1D70F}}=169.71$ and (c) $Re_{\unicode[STIX]{x1D70F}}=400$ . The vertical line marks the mean value.

Figure 11. Probability density function of the spanwise–streamwise aspect ratio of vortices: (a) $Re_{\unicode[STIX]{x1D70F}}=84.85$ , (b) $Re_{\unicode[STIX]{x1D70F}}=169.71$ and (c) $Re_{\unicode[STIX]{x1D70F}}=400$ . The vertical line marks the mean value.

By contrast, $Re$ has a much stronger effect on the spanwise vortex measurement. The distribution of the vortex aspect ratio $l_{z}^{+}/l_{x}^{+}$ (figure 11) becomes broader and high $l_{z}^{+}/l_{x}^{+}$ values are more frequently sampled with increasing $Re$ . The average aspect ratio also increases with $Re$ . Because $l_{x}^{+}$ is nearly the same, higher $l_{z}^{+}/l_{x}^{+}$ is solely due to the increasing spanwise measurement of the vortices. There are two major possible contributions to this increase: (i) streamwise vortices becoming increasingly bent and tilted towards the spanwise direction (see the $t=20$ panel of figure 5 for an illustration) and (ii) the increasing occurrence of curved and three-dimensional vortices such as hairpins. Quantitative assessment of these changes requires the statistics of vortices of different topologies.

Figure 12. Flow chart of the vortex classification procedure.

Figure 13. Definitions of the geometric quantities used in the vortex classification procedure of figure 12. Circles and squares represent $x$ - and $z$ -axis-points, respectively.

Figure 14. Classification of three-dimensional vortices: (a) hairpins, (b) hooks, (c) branch type A, (d) branch type B and (e) branch type C. Top row – representative examples from DNS; middle row – schematics of the vortex axis line; bottom row – streamwise profiles of the number of $x$ -axis points $N_{xap}$ and spanwise extent $D_{z}$ .

A new procedure is thus proposed to automatically classify the individual vortex axis lines, obtained from VATIP, according to their dimensions, geometry and topology. A flow chart of the procedure is provided in figure 12; the geometric quantities used in the procedure are shown in figure 13 and typical examples of different types in figure 14. The procedure consists of a series of binary decisions. First, all axis lines identified by VATIP are divided into two groups based on their streamwise measurement: those with $l_{x}^{+}\geqslant 50$ are considered as clear-cut vortices and smaller pieces are identified as fragments. This cutoff is smaller than the $150$ wall units used in Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997) because VATIP considers vortices with three-dimensional curvatures and the streamwise dimension does not necessarily account for the full vortex axis length. Those identified as vortices are further divided into streamwise versus three-dimensional types based on whether a significant spanwise segment can be found in the axis line. Note that any axis line identified by VATIP is formed by connecting axis points in any of the three dimensions (figure 13). The spanwise extent of all segments consisting of spanwise axis points only $l_{z,zap}^{+}$ are measured and if the maximum span $\max (l_{z,zap}^{+})\geqslant 25$ , it is determined that the vortex can no longer treated as a streamwise one.

Non-streamwise (or three-dimensional) vortices are further classified into several types based on the axis topology and geometry. A canonical hairpin is described as a vortex with two largely symmetric streamwise legs conjoining at its downstream head into a spanwise arc (figure 14 a). Many three-dimensional vortices bear some of the key features of a hairpin but significantly depart from its norm in other aspects. The classification procedure relies on two major geometric metrics of the identified vortex axis line (figure 13) to differentiate these different types: (i) the number of $x$ -axis points at a given $x$ position $N_{xap}$ and (ii) the spanwise separation between legs (again) at a given $x$ position $D_{z}$ . (For irregular vortices with more than two legs, e.g. column (c) of figure 14, $D_{z}$ is the spanwise separation between the two closest legs.) Variation of these two metrics with different $x$ positions is sketched for different vortex types in the bottom panels of figure 14. For a canonical hairpin (column (a)), $N_{xap}$ is 2 for the majority of the $x$ range although it may reduce to 1 at the beginning as the legs do not exactly match in length. Its $D_{z}$ starts high near the leg tips and gradually reduces to 0 as the legs fuse.

For any vortices deemed three-dimensional (versus streamwise) from the previous step, the procedure first checks the percentage of $x$ positions with only one $x$ -axis point $P_{x}(N_{xap}=1)$ ( $P_{x}(C)$ is the percentage of $x$ positions where a specific condition $C$ is satisfied) – if this quantity is ${>}80\,\%$ , i.e. for over 80 % of the vortex length it only has one leg, the vortex is a highly asymmetric variant of a hairpin where one of the legs is not clearly developed. This type is termed ‘hooks’ in our taxonomy. The other vortices have at least two legs, but there are various other branching configurations than the canonical hairpin. For example, in vortex packets where vortices are highly entangled and dynamically coalescing with one another, multi-legged – pitchfork-like – vortices are often observed (figure 14 c). If a vortex has more regions with three or more legs than those with two, i.e.  $P_{x}(N_{xap}>2)/P_{x}(N_{xap}>1)>50\,\%$ , it is identified as a branch type-A. Even vortices that only branch into two legs may appear significantly different from a canonical hairpin. For instance, the branch type-B (figure 14 d) looks more like a fusion between a quasi-streamwise vortex with a partial hairpin (or hook). This type of vortex was also reported in Robinson (Reference Robinson1991) and Brooke & Hanratty (Reference Brooke and Hanratty1993) and was believed to result from the spanwise shear dragging a side branch of a quasi-streamwise vortex to form an ‘arch’ on its side (Robinson Reference Robinson1991). The profile of $N_{xap}$ for this type shares some similarity with the hairpins, as both start with two legs which gradually merge. The main difference is that a hairpin ends with the arch where most axis points are counted, whereas in branch type-B the arch is followed by an extended streamwise segment downstream. Here, the $x$ projections of the centre-of-gravity (COG) of all $x$ -axis points $x_{COG,xap}$ and that of all $z$ -axis points $x_{COG,zap}$ are calculated. A canonical hairpin would be much ‘heavier’ at the downstream end, so if both COGs are at the upstream end, i.e.  $x_{COG,xap}<x_{mid}$ and $x_{COG,zap}<x_{mid}$ ( $x_{mid}\equiv (x_{max}-x_{min})/2$ being the $x$ coordinate of the middle point of the vortex axis line – see figure 13), the vortex is classified into branch type-B. In a similar scenario, when a side branch from a quasi-streamwise vortex protrudes towards the channel centre, because of the weaker transverse flows and higher mean velocity, the branch extends substantially downstream before any arch is formed. This is labelled as branch type-C in this study (figure 14 e) and identified by the criterion that $x_{max\text{-}D_{z}}>1.5x_{COG,N_{xap}>1}$ , where $x_{max\text{-}D_{z}}$ is the $x$ coordinate of the maximal branch separation $D_{z}$ and $x_{COG,N_{xap}>1}$ is the $x$ coordinate of the COG of the branched part of the vortex axis line (where $N_{x\text{-}cp}>1$ ). Finally, the remaining vortices – i.e. those predominated by two legs and with no substantial quasi-streamwise downstream segments – are classified as hairpins.

The criteria used in this classification procedure are mostly empirical. For starters, there is no physical ansatz supporting the classification of three-dimensional vortices into the five particular types listed in figure 14 – they are chosen solely based on empirical observations in our own and previous studies. Likewise, the dividing criteria and cutoff magnitudes used in the procedure (figure 12) are all chosen based on a combination of physical intuition and practical experience. For example, there is no physical basis as to how long a third leg needs to reach for a vortex to be considered a branch type-A (multi-legged) rather than a slightly modified hairpin. Indeed, the question of whether there is any fundamental difference between various types of branches and the canonical hairpin itself cannot be answered. The lack of objective vortex classification criteria is an inevitable consequence of the current limited knowledge of the complex vortex dynamics in wall turbulence. It is for this reason that an algorithm like VATIP is much needed. Future application of VATIP to a wider range of flow systems is anticipated to bring forth better experience and understanding of the characteristics of turbulent vortices, which will lead to a more standardized approach of vortex classification. Finally, we note that any vagueness in the current classification criteria does not affect the validity of any of the following discussion: e.g. the changes of all three-dimensional vortices show similar $Re$ dependence (figure 17) regardless of the further differentiation between hairpins and different branch types. In addition, from our test, changing the cutoff magnitudes by up to 50 % does not affect the comparison of vortex statistics between different $Re$ .

Figure 15. Distribution of vortex axis lines of different classes in a typical snapshot at $Re_{\unicode[STIX]{x1D70F}}=84.85$ : (a) hairpins, (b) hooks, (c) branches (all types), (d) fragments and (e) streamwise vortices. Each marker represents one axis point. Individual vortices are differentiated by colours and marker types.

Figure 16. Distribution of vortex axis lines of different classes in a typical snapshot at $Re_{\unicode[STIX]{x1D70F}}=400$ : (a) hairpins, (b) hooks, (c) branches (all types), (d) fragments and (e) streamwise vortices. Each marker represents one axis point. Individual vortices are differentiated by colours and marker types.

Figures 15 and 16 show the distributions of vortex axis lines, as identified by VATIP, of different classes for one typical snapshot of the lowest ( $Re_{\unicode[STIX]{x1D70F}}=84.85$ ) and the highest $Re$ ( $Re_{\unicode[STIX]{x1D70F}}=400$ ), respectively. Direct visual inspection of these images indicates that the VATIP algorithm together with the vortex classification procedure in figure 12 has successfully identified and extracted all types of vortices and sorted them properly according to their axis topology. This resonates with the earlier tests by STG in figure 8. Comparing different classes of vortices, streamwise ones still dominate at both $Re$ , but the method has no difficulty in finding all types of three-dimensional vortices. Unlike the case of boundary-layer flow where a so-called ‘forest’ of well-organized hairpins was observed (Wu & Moin Reference Wu and Moin2009), in our DNS results clear-cut hairpins are the minority compared with other three-dimensional configurations. In particular, the asymmetric hook type significantly outnumbers all other three-dimensional vortex types, which validates the earlier empirical notion in the literature about the prevalence of incomplete or one-legged hairpins in plane Poiseuille flow (Robinson, Kline & Spalart Reference Robinson, Kline and Spalart1989; Robinson Reference Robinson1991). On the other hand, the frequent appearance of various irregular branch types demonstrates the importance of iterative propagation in all three dimensions – a central element of VATIP.

Figure 17. Fraction of vortices of different types by vortex numbers. Only vortices with streamwise length $l_{x}^{+}>50$ are included.

Comparing between the two $Re$ , three-dimensional vortices (hairpins, hooks and branches) grow larger in size at higher $Re$ . This can be attributed to the increasing thickness of the wall layer (more wall units in the wall-normal direction) which allows these vortices to further lift up and develop to a higher altitude. They also become more populous at higher $Re$ . Indeed, even after factoring in the across-the-board increase of all vortices, the percentage share taken by three-dimensional vortices still steadily climbs. As shown in figure 17, with increasing $Re$ , quasi-streamwise vortices take up a lower percentage (despite a net increase in their number) and their share is replaced by all types of three-dimensional vortices. From $Re_{\unicode[STIX]{x1D70F}}=84.85$ to $400$ , the share of hooks increases by approximately 50 % and those of hairpins and branches more than double. Recall that hooks are often considered as asymmetric or incomplete hairpins, their slower growth (compared with symmetric hairpins and branches) suggests that they are likely the outcome of the insufficient development of hairpins and may become less important at higher $Re$ . Finally, in all $Re$ cases, complete hairpins are significantly outnumbered by their mutants – hooks and branches.

Near-wall vortex growth is often characterized as a lift-up process: the downstream end of the vortex becomes detached from the wall and rises towards the outer layer, where it can further burst and generate new disturbances (Hinze Reference Hinze1975; Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999). Lift-up extent of vortices at different wall layers can now be statistically analysed with the axis lines extracted by VATIP. Figures 18 and 19 show the joint PDF between the wall-normal positions of the heads and tails of all quasi-streamwise and three-dimensional vortices at different $Re$ . The head position $y_{head}^{+}$ is measured as the highest wall-normal position of all axis points, which is normally found at the downstream end; likewise, the tail position $y_{tail}^{+}$ is the lowest position normally found at the upstream end. Obviously, the distribution can only sample the upper-left triangle of the domain. Vortices that have not lifted up are represented by the diagonal where the head position is levelled with the tail and regions closer to the ordinate, i.e.  $y_{head}^{+}\gg y_{tail}^{+}$ , and correspond to highly lifted up vortices.

Figure 18. Joint probability density function (PDF) between the $y^{+}$ positions of the tail and the head of vortices at $Re_{\unicode[STIX]{x1D70F}}=84.85$ : (a) quasi-streamwise and (b) three-dimensional vortices.

Figure 19. Joint probability density function (PDF) between the $y^{+}$ positions of the tail and the head of vortices at $Re_{\unicode[STIX]{x1D70F}}=400$ : (a) quasi-streamwise and (b) three-dimensional vortices.

For the lower $Re_{\unicode[STIX]{x1D70F}}=84.85$ case (figure 18), both the head and tail positions of quasi-streamwise vortices (panel a) concentrate at $10\lesssim y^{+}\lesssim 50$ : i.e. within or near the buffer layer. Three-dimensional vortices (hairpins, hooks and branches) have a higher altitude and their distribution peaks at $(15,80)$ : i.e. the tail (legs) stretches deep into the buffer layer while the head (arc in the case of hairpins) rises up into the log-law layer. In terms of distribution, the tails are concentrated at $y^{+}<25$ whereas the heads are found in a much broader range extending from $y^{+}=25$ to $y^{+}>80$ . These observations can all carry over to the higher $Re_{\unicode[STIX]{x1D70F}}=400$ (figure 19) where, in addition, the larger number of wall units in the $y$ direction allows more room for vortex growth and their lift-up extent is easier to observe. Most quasi-streamwise vortices (figure 19 a) are lying flat in the buffer layer (concentration peak in the lower-left corner) but two more concentration bands can be spotted: one lies along the ordinate up to $y^{+}=100$ , indicating that a small fraction of streamwise vortices can lift up to the log-law layer; the other lies along the diagonal to even higher $y^{+}$ , indicating the existence of flat-lying vortices at higher altitudes. Both bands are also clearly visible in the three-dimensional case (figure 19 b) but the vertical one is stronger over a broader range of $y^{+}$ , meaning that these vortices are more likely to lift up and their heads can reach various altitudes. Vortex activities at $y^{+}>250$ are much weaker and thus not included in figure 19. (Alternatively, following the example of Lozano-Durán et al. (Reference Lozano-Durán, Flores and Jiménez2012), Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014), one may apply non-uniform $Q$ threshold values – with lower thresholds for the bulk – for a complete picture.) Observations from this analysis largely confirm the earlier empirical depiction by Robinson (Reference Robinson1991) that quasi-streamwise vortices dominate the buffer layer and hairpin-like vortices are more likely to be found in the log-law layer and beyond. Robinson (Reference Robinson1991) conceived the log-law layer to be comprised of a mix of streamwise and hairpin vortices, whereas we are able to more clearly show that streamwise vortices are only concentrated in the lower log-law layer ( $y^{+}<50$ ) and three-dimensional vortices can rise up to a variety of altitudes.

4.4 Vortex organization through clustering analysis

Previous observations of LSMs and VLSMs ignited an immense interest among researchers in understanding the organization patterns of coherent structures (Jiménez Reference Jiménez1998; Kim & Adrian Reference Kim and Adrian1999; Lee et al. Reference Lee, Lee, Choi and Sung2014). Given the specific information, available from VATIP, about the location and conformation of axis lines representing individual vortices, we adapt the DBSCAN (density-based spatial clustering of applications with noise) algorithm (Ester et al. Reference Ester, Kriegel, Sander and Xu1996) – a widely used clustering analysis method in data mining and machine learning – to VATIP results for understanding the clustering patterns of vortices. (Structures classified as fragments according to figure 12 are not considered in this analysis.)

The standard DBSCAN algorithm groups scattered points in space into clusters based on their spatial proximity and mutual relationship. Two points that are close to each other (within a cutoff distance $\unicode[STIX]{x1D700}$ ) are considered as neighbours. Points inside a cluster are known to have many neighbours. Points with at least $N_{c,min}$ neighbours are thus labelled as ‘core points’ and all interconnected (in the sense of mutually neighbouring) core points are grouped into one cluster. (Both $\unicode[STIX]{x1D700}$ and $N_{c,min}$ are user-specified parameters.) If a point does not qualify as a core point by itself but is neighbour to one or more core points, it is labelled as a ‘border point’ which resides on the surface of a cluster. Border points are grouped to the same cluster as their nearest neighbouring core point. Points that do not neighbour any core points and are not core points themselves are isolated outliers not belonging to any cluster.

Since the VATIP output contains not simple sizeless points but complex axis lines representing vortex geometry and topology, the simple distance criterion used for neighbour identification needs to be adapted. We consider two vortices to be neighbours if the minimum distance between any two axis points – one on each axis line – does not exceed $\unicode[STIX]{x1D700}=4r_{v}$ which is only slightly larger than the detection cone diameter used in VATIP tracking ( $2\times 1.5r_{v}=3r_{v}$ ). We have tested a wide range of $\unicode[STIX]{x1D700}$ and found that for $\unicode[STIX]{x1D700}$ as low as $3.5r_{v}$ nearly all vortices in the domain, from both sides of the channel, are interconnected into the same neighbour network: i.e. for $N_{c,min}=1$ and any $\unicode[STIX]{x1D700}\geqslant 3.5r_{v}$ , the DBSCAN algorithm will identify one supersized cluster that includes nearly all vortices. The fact that a cutoff distance of the same order as the vortex diameter would connect all vortices is not surprising, considering the level of crowdedness found in their distribution (see figures 15 and 16). Ideally, we would also need to test the $\unicode[STIX]{x1D700}$ -dependence at other $N_{c,min}$ levels. However, our priority is to understand the importance of multi-vortex cooperation (instead of inter-vortex distance, which we knew would be close). Therefore, given the limited scope of this investigation, we will focus here on the $N_{c,min}$ -dependence of the clustering results at a constant $\unicode[STIX]{x1D700}=4r_{v}$ .

Figure 20. Dependence of DBSCAN clustering analysis on $N_{c,min}$ ( $Re_{\unicode[STIX]{x1D70F}}=169.71$ ): left/red/circle – number fraction of all vortices grouped into clusters; right/blue/square – number fraction of vortices in the largest identified cluster.

With increasing $N_{c,min}$ , less vortices are qualified as core vortices (counterpart to core points in the standard DBSCAN) and more become isolated outliers. This first leads to the shrinkage of all clusters: for the same total number of vortices $N_{v,tot}$ , the number of vortices assigned to clusters $N_{v,clus}$ decreases (figure 20). At $N_{c,min}=2$ (lowest level shown in figure 20), $N_{v,clus}/N_{v,tot}$ starts at close to $1$ (nearly all vortices are grouped into clusters) and steadily drops afterwards towards $0$ . The number of vortices contained in the largest cluster $N_{max,clus}$ is also calculated. In figure 20, $N_{max,clus}/N_{v,tot}$ starts at ${\approx}0.5$ , because the channel flow geometry has two boundary layers (near each wall) and at the lowest $N_{c,min}$ vortices near each wall are nearly all grouped into one super-cluster. The decline pattern of this profile is very different from that of $N_{v,clus}/N_{v,tot}$ – it drops sharply in a small window of $N_{c,min}=8\sim 14$ with the steepest slope found between $N_{c,min}=10$ and $12$ . This faster decline cannot be solely accounted for by the overall reduction of qualifying core vortices (otherwise $N_{max,clus}/N_{v,tot}$ would have the same slope as $N_{v,clus}/N_{v,tot}$ ). Indeed, the steeper descent indicates a sudden disintegration of the dominant clusters into smaller pieces. As $N_{c,min}$ increases beyond ${\approx}8$ , some vortices in the structure, which are not as highly intertwined as most others in the vortex cluster network, are disqualified as core vortices. Removing those ‘bridge’ vortices dismantles the cluster network into several well-defined and strongly coupled clusters that are much smaller in size. This effect is most clearly seen from the ratio between these two profiles, plotted in figure 21. For all three $Re_{\unicode[STIX]{x1D70F}}$ tested, $N_{max,clus}/N_{v,clus}$ is initially flat at low $N_{c,min}$ , indicating that, within this regime, drops in both profiles in figure 20 are attributed to the overall reduction of clustered vortices. Disintegration of the dominant clusters starts when the $N_{max,clus}/N_{v,clus}$ profile turns downwards, which for $Re_{\unicode[STIX]{x1D70F}}=169.71$ occurs at $N_{c,min}\approx 8$ . The process finishes as the curve reaches its minimum at $N_{max,clus}/N_{v,clus}=0.1\sim 0.15$ : the domain is now populated by $O(10)$ well-defined clusters with comparable size (see figure 22 c for roughly half of the clusters on one side of channel). After the minimum the profile rises again as a result of smaller clusters gradually being eliminated by the increasingly stringent $N_{c,min}$ cutoff.

Figure 21. Number fraction of vortices included in the largest cluster among vortices in all clusters as a function of $N_{c,min}$ .

Figure 22. Distribution of vortex clusters (on one side of the channel and from a typical instantaneous flow image) at $Re_{\unicode[STIX]{x1D70F}}=169.71$ identified by DBSCAN with (a) $N_{c,min}=2$ , (b) $N_{c,min}=12$ , (c) $N_{c,min}=15$ and (d) $N_{c,min}=21$ . Individual clusters are differentiated by colour. Black lines show the contours of $y$ -average Reynolds shear stress $\bar{\unicode[STIX]{x1D70F}}_{xy}$ (4.4) at 11 equispaced levels from $0.5\times 10^{-5}$ to $1.25\times 10^{-4}$ ; higher contour line density corresponds to higher magnitudes.

The vortex cluster configuration during this disintegration process with increasing $N_{c,min}$ is shown in figure 22 for $Re_{\unicode[STIX]{x1D70F}}=169.71$ . Consistent with our earlier analysis, at $N_{c,min}=2$ all vortices on one side of the channel are interconnected (by their neighbouring network) into a super-cluster. Disintegration of the network is observed at $N_{c,min}=12$ ( ${>}8$ where it starts). At $N_{c,min}=15$ (figure 22 c), $N_{max,clus}/N_{v,clus}$ reaches its minimum (figure 21) and the disintegration process has completed with a number of clear vortex clusters remaining unbroken. Much space can be found between the clusters where the turbulent flow field is occupied by unclustered vortices. The Reynolds shear stress, averaged over the wall-normal ( $y$ ) direction,

(4.4) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D70F}}_{xy}=-\int _{0}^{1}v_{x}^{\prime }(x,y,z,t)v_{y}^{\prime }(x,y,z,t)\,\text{d}y\end{eqnarray}$$

is shown with contour lines in the images. Spots with strong $\bar{\unicode[STIX]{x1D70F}}_{xy}$ (dense contour lines) are found within or immediately around these vortex clusters, indicating their strong contribution to the Reynold stress generation. Interestingly, the characteristic length (in the $x$ direction) of these clusters seems to be between 500 and 1500 wall units, which is at the same level as the typical streamwise length scales of LSMs reported in the literature (Adrian Reference Adrian2007; Lee et al. Reference Lee, Lee, Choi and Sung2014). The existence of such clusters consisting of a large number ( $O(10)$ or higher; see figure 23) of vortices strongly intertwined through multi-body interactions ( ${>}N_{c,min}=15$ neighbours, in the case of $Re_{\unicode[STIX]{x1D70F}}=169.71$ , with close contacts between their axis lines for the core vortices) is consistent with the hypothesis that LSMs are results of the cooperative dynamics involving many vortices organized as ‘packets’ (Kim & Adrian Reference Kim and Adrian1999). However, these ‘packets’, as discussed below, are not composed of well-aligned hairpin vortices with their classical shape. In addition, clear evidence for clustering is found in this study for $Re_{\unicode[STIX]{x1D70F}}$ all the way down to below $100$ , suggesting that cooperative dynamics between vortices is a universal feature for wall turbulence not limited to the high- $Re$ regime. Meanwhile, VLSMs are often conjectured to occur at a higher level of organization involving the alignment of multiple LSMs (Kim & Adrian Reference Kim and Adrian1999; Lee et al. Reference Lee, Lee, Choi and Sung2014). This would correspond to the cooperative organization involving multiple vortex clusters in this study. The length scale of VLSMs is comparable to or larger than the current domain size and they were previously studied mostly at much higher $Re_{\unicode[STIX]{x1D70F}}$ ( $O(10^{3})$ ). For these reasons, they are not discussed here. As $N_{c,min}$ further increases to 21, all clusters are now eliminated except the strongest one, which has shrunken in size but still clearly marks the location of strong Reynolds stress activity.

Note that the term ‘cluster’ has a different meaning here than that in earlier studies of three-dimensional vortex analysis, such as Del Álamo et al. (Reference Del Álamo, Jiménez, Zandonade and Moser2006) where clusters referred to the interconnected structures with overlapping vortex volumes identified by the scalar identifier ( $\unicode[STIX]{x1D6E5}$ in that study and $Q$ here), regardless of the individual identities of vortices or their conformation and topology. In our analysis, a cluster is defined as individual vortices grouped together based on the existence of a mutually interacting (neighbouring) network between multiple vortex objects rather than a pure spatial-proximity criterion. There are likely close connections between these two interpretations, but at this point, a direct comparison is not possible, because, as further discussed in § 4.6, the current VATIP algorithm can only capture a subset of structures analysed in Del Álamo et al. (Reference Del Álamo, Jiménez, Zandonade and Moser2006) that are directly generated from the lift-up-from-wall process. The strength of the current approach is its access to the information of individual constituting vortices, which we discuss below.

Figure 23. Probability density function of the number of vortices in a single cluster at (a) $Re_{\unicode[STIX]{x1D70F}}=84.85$ , (b) $Re_{\unicode[STIX]{x1D70F}}=169.71$ and (c) $Re_{\unicode[STIX]{x1D70F}}=400$ . The vertical line marks the mean value.

Figure 24. Probability density function of the number of vortices in a single cluster normalized by the total number of vortices in the domain at (a) $Re_{\unicode[STIX]{x1D70F}}=84.85$ , (b) $Re_{\unicode[STIX]{x1D70F}}=169.71$ and (c) $Re_{\unicode[STIX]{x1D70F}}=400$ . The vertical line marks the mean value.

Unless otherwise noted, we pick the $N_{c,min}$ value at the minimum in each $N_{max,clus}/N_{v,clus}$ curve (figure 21) – i.e.  $N_{c,min}=11,15$ and $22$ for $Re_{\unicode[STIX]{x1D70F}}=84.85,169.71$ and $400$ , respectively – for the DBSCAN analysis, which is the point where the percolating super-cluster has been fully disintegrated into unbreakable clusters while most individual clusters are not yet eliminated. (This choice is in the same spirit as the percolation analysis of Lozano-Durán et al. (Reference Lozano-Durán, Flores and Jiménez2012) explained in § 4.5.) The PDF of the number of vortices constituting a single cluster $N_{v,single}$ , shown in figure 23, is clearly skewed to the right with the most probable value at $O(10)$ but in some extreme cases with $O(100)$ vortices in each cluster. The average $N_{v,single}$ increases with $Re$ and is ${\approx}18,26$ and $35$ , respectively, from the lowest to the highest $Re_{\unicode[STIX]{x1D70F}}$ tested. Note, however, that the total number of vortices in the domain $N_{v,tot}$ also increases with $Re$ . Indeed, when $N_{v,single}$ is normalized by $N_{v,tot}$ (figure 24), the distribution profile becomes nearly the same between different $Re$ (mean value at $0.014,0.013$ and $0.012$ for $Re_{\unicode[STIX]{x1D70F}}=84.85,169.71$ and $400$ , respectively). Since the domain size of different $Re$ is kept the same in inner units (table 1), this observation, that the average cluster size in terms of the number fraction of vortices in each cluster (out of all vortices filling the domain), remains roughly constant, suggests that the cluster size is more or less the same in inner units (i.e. streamwise length within $500\sim 1500$ ; see figure 22) within the $Re$ range tested (further demonstrated in figure 27).

Figure 25. Number fraction of vortices of different types grouped into clusters out of the total number of all vortices of the same type in the domain ( $Re_{\unicode[STIX]{x1D70F}}=169.71$ ).

Figure 26. Probability density function of the number fraction of three-dimensional vortices in a single cluster at (a) $Re_{\unicode[STIX]{x1D70F}}=84.85$ , (b) $Re_{\unicode[STIX]{x1D70F}}=169.71$ and (c) $Re_{\unicode[STIX]{x1D70F}}=400$ . The vertical line marks the mean value.

Comparing vortices of different types (figure 25), hairpins and branches are more likely to be included in a cluster than both streamwise and hook vortices. Since hooks are essentially incomplete or asymmetric hairpins that are also highly lifted up, this indicates that the large dimensions of hairpins and branches are likely the key factor determining their higher clustering tendency. In particular, their wide span in the $z$ direction exposes them to vortices from a wider flow region, which enables them to play a central role in stitching more vortices into a cluster. As seen in figure 26, within a single cluster, a significant fraction of the vortices belong to the three-dimensional classes (hairpins, branches or hooks). This fraction increases with $Re$ and at $Re_{\unicode[STIX]{x1D70F}}=400$ , on average nearly half of the vortices in each cluster are three-dimensional ones. However, note that hooks and branches significantly outnumber canonical hairpins (figure 17), the hypothesized picture of packets of clean-cut hairpins (Adrian Reference Adrian2007) forming the LSMs is not seen at least at the current $Re$ range.

Figure 27. Representative images of vortex clusters identified by DBSCAN at (a,b) $Re_{\unicode[STIX]{x1D70F}}=84.85$ , (c,d) $Re_{\unicode[STIX]{x1D70F}}=169.71$ and (e,f) $Re_{\unicode[STIX]{x1D70F}}=400$ . Isosurfaces show all vortices in the viewable region (colour varies from light to dark with $y^{+}$ ); red dots show the axis lines of vortices in the identified cluster.

Direct images of representative vortex clusters are shown in figure 27 where vortices forming the particular cluster are highlighted by explicitly showing their axis lines. Consistent with our earlier observations, these clusters (at different $Re$ ) all have a streamwise length in the range of $500{-}1500$ wall units. Two typical organization configurations are observed. In the first (panels a,f), different vortices forming the cluster have their axis lines braided together along the streamwise direction. These clusters have a shape of twisted doughnuts and they remain slender (narrow in the $z$ direction) while extending downstream for $O(1000)$ wall units. For the second type, which is more frequently observed, other than the downstream twisting, the clusters also expand in the spanwise direction by connecting more vortices through the wider vortex types (hairpins and branches).

Finally, we note that the analysis of vortex clustering and organization in this section is still preliminary and limited in scope. It is intended to provide some first insight into how the axis-line information extracted by VATIP can be used to address some of the most important outstanding questions in turbulent dynamics (Jiménez Reference Jiménez2018). Further research is needed to better connect these observations with the existing conceptual models and results from other structure analysis techniques.

4.5 Determination of parameters and settings in VATIP

After presenting the main results, we are now ready to assess the robustness of VATIP tracking outcomes and discuss the procedure for choosing its parameters and settings. There are two major adjustable parameters in the method: (i) the threshold magnitude $Q_{threshold}$ for vortex identification with the $Q$ -criterion and (ii) the cutoff cone radius $r_{cone}$ used in the axis-line propagation search (figure 2). In addition, sensitivity to grid sizes and the selection of the search starting plane will also be examined.

Figure 28. Vortex disintegration with increasing $Q_{threshold}=HQ_{rms}$ at $Re_{\unicode[STIX]{x1D70F}}=169.71$ : (a) $H=0.2$ , (b) $H=0.4$ and (c) $H=0.7$ . Interconnected vortex tube structures are coded with the same colour. For clarity, only the largest vortices that cumulatively account for 80 % (for a,b) or 60 % (for c) of the total vortex volume are shown. For (b,c), only vortices from the bottom half of the channel are shown.

The choice of the threshold for $Q$ (or any other vortex identifier) has been widely discussed in the literature for the purpose of vortex visualization. It is a common practice to choose a threshold in proportion to its r.m.s. value in the flow field

(4.5) $$\begin{eqnarray}Q_{threshold}\equiv HQ_{rms},\end{eqnarray}$$

where $H$ in this study is chosen based on the percolation analysis, proposed in Lozano-Durán et al. (Reference Lozano-Durán, Flores and Jiménez2012) (from which we also borrowed the notation $H$ ). When $H$ is low, the identified vortex regions interconnect with one another and form a percolating network across the domain (figure 28 a). With increasing $H$ , the ‘necks’ bridging stronger vortex cores gradually break to reveal individual groups of vortices (figure 28 b,c); meanwhile, a higher threshold also erases many weaker vortices from the view. A percolation diagram (figure 29) plots the ratio of the volume occupied by the largest interconnected structure ( $V_{max}$ ) to that of all vortex regions ( $V_{tot}$ ) as a function of $H$ . This value starts at $1$ at the low $H$ end where all structures are interconnected into a complete percolating network. Increasing $H$ reduces both $V_{max}$ and $V_{tot}$ , whereas the decrease of their ratio $V_{max}/V_{tot}$ reflects the disintegration of larger interconnected structures into smaller separate pieces. The latter clearly dominates the window of $0.3\lesssim H\lesssim 0.7$ where the sudden disintegration of the largest structure into multiple objects is reflected in a steep descent in $V_{max}/V_{tot}$ . This window provides the best choices for $H$ , in which individual objects are separate and identifiable whereas the most important vortex structures are not yet erased (as what will happen at higher $H$ ) (Lozano-Durán et al. Reference Lozano-Durán, Flores and Jiménez2012).

Figure 29. Percolation diagram for $Re_{\unicode[STIX]{x1D70F}}=169.71$ . The vertical dashed line marks the $HQ=0.4Q_{rms}$ used in this study.

Figure 30. Average number of vortices (excluding fragments) detected by VATIP in each instantaneous flow field as a function of $H$ (with $r_{cone}=1.5r_{v}$ ).

Also as a result of the competition between vortex shrinkage and disintegration, the average number of vortices identified by VATIP in each flow domain $N_{v}$ displays a non-monotonic dependence on $H$ (figure 30). $N_{v}$ initially increases with $H$ , reflecting the splitting of vortex objects. After reaching a maximum at $H\approx 0.2$ , $N_{v}$ starts to decline because the identified vortices shrink in size with increasing $H$ and are increasingly categorized as fragments. This effect becomes more dominant at higher $H$ after the disintegration of the percolating network. Within the acceptable range of $H=0.3\sim 0.7$ – as identified above by the percolation analysis – the drop of $N_{v}$ is relatively mild ( ${\lesssim}10\,\%$ for the two higher $Re_{\unicode[STIX]{x1D70F}}$ ). More importantly, as shown later, all major conclusions from the study remain intact within this range of $H$ . It is worth noting that the peak of $N_{v}$ is found out of this range at a slightly lower $H$ : i.e. the main vortex disintegration events are detected at a slightly lower $H$ using VATIP than the percolation analysis. This is because VATIP is more sensitive to the breakage between vortex structures: neighbouring vortices could well overlap in their shells and be grouped into the same interconnected structure in figure 28 while their axis lines do not have topological connection. This also explains, in part, why $N_{v}$ does not start from $1$ at $H=0$ in figure 30. (Another reason – further discussed in § 4.6 – is that VATIP is designed with wall-generated vortices in mind and may not fully capture the connections between weaker and more isotropic vortex structures in the bulk region, which are only unveiled at very low $H$ .) Taking this into account, the optimal $H$ pick should be slightly lower than that for the steepest descent in figure 29. Therefore, $H=0.4$ is chosen in this study for VATIP tracking (in comparison with $H=0.7$ used in Zhu et al. (Reference Zhu, Schrobsdorff, Schneider and Xi2018) for vortex visualization). Note that in figure 30, $N_{v}$ is nearly constant in the range of $H=0.3\sim 0.4$ .

Figure 31. Average number of vortices (excluding fragments) detected by VATIP in each instantaneous flow field as a function of $\unicode[STIX]{x1D701}$ (with $Q_{threshold}=0.4Q_{rms}$ ).

The cone size is chosen based on the average radius of the vortex tubes (3.1)

(4.6) $$\begin{eqnarray}r_{cone}\equiv \unicode[STIX]{x1D701}r_{v},\end{eqnarray}$$

where $\unicode[STIX]{x1D701}$ is expected to be larger than (but of the same order of magnitude as) $1$ to account for vortex size variations. The average number of vortices identified by VATIP in each flow domain $N_{v}$ is also non-monotonic with increasing $\unicode[STIX]{x1D701}$ (figure 31). For $\unicode[STIX]{x1D701}<1$ , many well-defined vortices are broken into pieces and excluded as fragments. Meanwhile, at $\unicode[STIX]{x1D701}\gg 1$ , false connection between separate vortices becomes more common and $N_{v}$ decreases with $\unicode[STIX]{x1D701}$ . Interestingly, for all $Re_{\unicode[STIX]{x1D70F}}$ tested, $N_{v}$ reaches its maximum at exactly $\unicode[STIX]{x1D701}=1$ , indicating that $r_{v}$ calculated by (3.1) does provide an accurate measurement of the vortex radius. We recommend the range of $\unicode[STIX]{x1D701}=1.2\sim 1.6$ for VATIP where the decline of $N_{v}$ is modest (compared with higher $\unicode[STIX]{x1D701}$ ) and, more importantly, all major physical observations are consistent with changing $\unicode[STIX]{x1D701}$ (shown below). For $\unicode[STIX]{x1D701}=1.5$ used in this study, the resulting $r_{cone}^{+}$ is approximately $15$ , $16$ and $18$ wall units for $Re_{\unicode[STIX]{x1D70F}}=84.85$ , $169.71$ and $400$ , respectively. For comparison, Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997) used $r_{cone}^{+}\approx 10$ for their streamwise-only search at $Re_{\unicode[STIX]{x1D70F}}\approx 180$ , which is equivalent to $\unicode[STIX]{x1D701}\approx 1$ . The larger $\unicode[STIX]{x1D701}$ used in VATIP is necessitated by the expansion of the search to all three spatial dimensions. First, dislocation between successive axis points is typically larger around the bends or turns of the axis line, which does not occur in a unidirectional search along nearly straight lines. Second, inclusion of highly lift-up hairpin-like vortices extends the search deep into the log-law layer, where the vortex diameters are often larger compared with the streamwise vortices in the buffer layer. Lastly, a streamwise search only looks for new axis points in the $yz$ -plane where the numerical grids are typically more refined (than the $x$ direction) in DNS. Searches in other directions need to accommodate axis-point dislocation in the $x$ direction: with the coarser mesh of Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997), $r_{cone}^{+}=10$ covers less than one $x$ -grid spacing – $\unicode[STIX]{x1D6FF}_{x}^{+}=17.7$ : i.e. no dislocation in $x$ would be allowed. Our experience also shows that $\unicode[STIX]{x1D701}=1$ would break well-defined hairpin vortices (such as in STG) into pieces.

Figure 32. Number fraction of different vortex types (excluding fragments) as a function $Re$ identified by VATIP with different $H$ (and a constant $\unicode[STIX]{x1D701}=1.5$ ): (a) streamwise, (b) hairpin, (c) hook and (d) branch.

Figure 33. Number fraction of different vortex types (excluding fragments) as a function $Re$ identified by VATIP with different $\unicode[STIX]{x1D701}$ (and a constant $H=0.4$ ): (a) streamwise, (b) hairpin, (c) hook and (d) branch.

The similarity between $N_{v}$ profiles of different $Re_{\unicode[STIX]{x1D70F}}$ in both figures 30 and 31 suggests the robustness of VATIP at least within the $Re$ range tested. In figures 32 and 33, it is clear that within the recommended ranges of $H=0.3\sim 0.7$ and $\unicode[STIX]{x1D701}=1.2\sim 1.6$ , changes in vortices of different types with increasing $Re$ follow the same consistent trend with different $H$ and $\unicode[STIX]{x1D701}$ . Clear disruption to the trend is only observed in cases well out of these ranges: most notably $H=0.2$ in figure 30 and $\unicode[STIX]{x1D701}=0.8$ and $3.0$ in figure 31. Quantitative magnitudes of the profiles do depend on $H$ and $\unicode[STIX]{x1D701}$ , which is very much expected. As illustrated in figure 34, adjusting these parameters inevitably changes the lengths of vortex branches and legs, to which the classification scheme of figure 12 is very sensitive: missing one axis point at the branch end could result in a vortex being classified as a hook rather than hairpin, or even a fragment rather than a vortex. Nevertheless, quantitative differences between curves are significantly smaller (mostly contained within a few percentage points) in the ranges of $H=0.3\sim 0.7$ and $\unicode[STIX]{x1D701}=1.2\sim 1.6$ , compared with those out of these ranges. The physical observation made in figure 17 are completely robust when acceptable $H$ and $\unicode[STIX]{x1D701}$ are used.

Figure 34. Effects of parameters and settings on VATIP tracking results in a representative flow region at $Re_{\unicode[STIX]{x1D70F}}=169.71$ (note: panel (f) is out of the recommended range of $H$ ). Dots represent axis points identified by VATIP (different colours for different vortices) and grey tubes are the isosurfaces of $Q=HQ_{rms}$ (same $H$ used in VATIP). Partial vortices (with parts extending out of the view box) are not included in VATIP tracking.

The robustness of VATIP is most clearly demonstrated in figure 34 where different parameters and settings are tested and compared for a same flow region with various vortex configurations. Compared with the standard case (panel a) with $H=0.4$ and $\unicode[STIX]{x1D701}=1.5$ , changing $\unicode[STIX]{x1D701}$ to $1.2$ (panel d) or changing $H$ to $0.6$ (panel e) brings little noticeable difference. In panel (f), $H$ is further increased to $0.8$ (beyond the recommended range), which only causes the identified vortex tubes to shrink in size, and VATIP still faithfully captures their axis lines. Note that the brown vortex at the lower-left corner has changed from a curved shape to a linear shape at $H=0.8$ owing to the erosion of one of its legs, which well illustrates how changing parameters affect the number of vortices classified into each category (including fragments). These changes, however, do not reflect the reliability of VATIP itself. We have also doubled the resolutions in the $x$ - and $z$ -dimensions (panel b; the flow field is interpolated to the finer grid before VATIP analysis) and changed the starting search plane in $x$ and $z$ directions (first steps in subroutines 1 and 2 of figure 3; with no translational symmetry, $y$ -direction searches always start from the walls – also see discussion in § 4.6) from the first planes ( $x=0$ or $z=0$ ) to the middle planes ( $x=L_{x}/2$ or $z=L_{z}/2$ ). Both do not lead to any discernible difference in the tracking outcome. This observation is general: at $Re_{\unicode[STIX]{x1D70F}}=169.71$ , the average number of vortices in each flow domain (excluding fragments) identified by VATIP $N_{v}=2178$ , whereas the high resolution case has $N_{v}=2193$ and the mid-plane start case has $N_{v}=2182$ . In both cases, the difference is way less than 1 %.

Figure 35. Joint PDF between the $y^{+}$ positions of the tail and the head of streamwise vortices at $Re_{\unicode[STIX]{x1D70F}}=169.71$ under different VATIP parameters and settings (note: panel (f) is out of the recommended range of $H$ ).

Figure 36. Joint PDF between the $y^{+}$ positions of the tail and the head of three-dimensional vortices at $Re_{\unicode[STIX]{x1D70F}}=169.71$ under different VATIP parameters and settings (note: panel (f) is out of the recommended range of $H$ ).

Finally, echoing the observations in figures 18 and 19, we examine the effects of VATIP parameters and settings on vortex conformation statistics in figures 35 and 36 – this time at $Re_{\unicode[STIX]{x1D70F}}=169.71$ . Despite the small quantitative differences – which, as discussed above, are inevitable as vortices are classified based on quantitative metrics, the qualitative picture is well preserved for all cases shown, including the $H=0.8$ case which is out of the recommended range. Similar to the earlier observations at other $Re_{\unicode[STIX]{x1D70F}}$ , the distribution of streamwise vortices is highly concentrated in the lower-left corner corresponding to the buffer layer. Weaker, but noticeable, concentration bands extend along both the diagonal and the ordinate, reflecting the flat-lying and lifted-up streamwise vortices, respectively. By contrast, three-dimensional vortices are predominantly lifted up, with their concentration peak found well in the log-law layer. Changing resolution or the starting plane shows little effect on these distributions, whereas adjusting $H$ or $\unicode[STIX]{x1D701}$ more directly affects vortex classification and thus causes some subtle changes in the contour shapes, especially at low density levels.

4.6 Discussion: limitations and future development

Recall from § 3: the algorithm of VATIP is built on the premise that vortices are wall generated, starting with segments or ‘legs’ that align along the $x$ direction (most often in the buffer layer but the algorithm does not impose this restriction) and can lift up to higher- $y^{+}$ layers to bend, curve or branch. VATIP always initiates the propagation points in the $x$ -lying legs and later allows them to move away from the walls (in the $y$ -direction search) and swing sideways (in the $z$ -direction search). For canonical hairpins, the axis lines initiate from both streamwise legs which rise at the downstream end and merge in the middle along the $z$ -direction. Branched vortices are found in a similar manner with one propagation point planted in each streamwise leg and the growing legs (or more appropriately for the branch type – arms) will eventually merge after a limited number of iterations. As demonstrated above (comparing the $Q$ -isosurfaces and VATIP-identified axis lines), this algorithm faithfully captures nearly all vortices identified by the $Q$ -criterion in wall turbulence for $Re_{\unicode[STIX]{x1D70F}}\leqslant 400$ tested in this study.

Recent evidences have indicated that at high $Re$ and large $y^{+}$ , vortices can be generated in the absence of wall interaction (Del Álamo et al. Reference Del Álamo, Jiménez, Zandonade and Moser2006; Jiménez Reference Jiménez2013). These vortices can deviate significantly from this premise: they are nearly isotropic (segments are equally likely to align with any direction) and often highly branched (multiple arms with complex connection topology). The current algorithm would not perform as well on those structures. First, the requirement on initialization in the $x$ -direction search only will undoubtedly bias the resulting axis line to have better sampling of the $x$ -lying segments. For instance, in a strictly $y$ -aligned segment with no connection to any $x$ -segment at the bottom, the current algorithm would still capture a point in the middle that is an $yz$ -planar maximum. It would skip the $x$ -direction propagation (for the lack of other $x$ -axis points) and the next step of $y$ -search would propagate away from the wall only. The end result is missing one part of the segment below the initial point. Second, with only one propagation point in each initial $x$ -segment, the algorithm would struggle with highly branched configurations because the propagation will only pick one of the directions to proceed after each junction. This is not much a problem for wall-generated branched vortices at moderate $Re$ (the focus of this study), because these structures mostly consist of a few conjoining arms, each of which can be traced back to a streamwise leg (i.e. starfish or octopus shaped). The current algorithm has been shown to capture these structures well. However, for detached structures at higher $Re_{\unicode[STIX]{x1D70F}}$ and higher $y^{+}$ (see, e.g. figure 6 of Del Álamo et al. (Reference Del Álamo, Jiménez, Zandonade and Moser2006)) with more complex configurations, it will likely miss some of the bridges connecting different arms, if they do not happen to align in the $x$ direction, and falsely break them into pieces.

There seems to be a easy remedy in sight, that as long as we relax these constraints to allow a truly multidirectional tracking – i.e. axis lines are initiated in all three dimensions and propagation is allowed to follow all branches after each junction – we should be able to capture these isotropic and highly branched structures. The problem, however, is its insurmountable side effect: relaxation of the current constraints will inevitably lead to massive false identification and false connection, making the tracking result next to meaningless. There are two major sources of this problem. First, not all planar $Q$ maxima belong to a vortex axis. Consider a simple linear streamwise vortex as an elongated ellipsoid, the true axis line aligns with its major axis and consists of $yz$ -planar $Q$ maxima. However, the minor axes also contain $Q$ maximum points (in $xy$ - and $xz$ -planes), which do not belong to any vortex axis line. This is further compounded by fluctuations in the $Q$ fields, which may create $Q$ maxima unrelated with any actual vortex. The total number of planar maximum points identified in an $L_{x}^{+}\times L_{z}^{+}=4000\times 800$ flow domain ranges from ${\sim}80\,000$ to ${\sim}350\,000$ (for the lowest and highest $Re_{\unicode[STIX]{x1D70F}}$ tested, respectively). Only $30\sim 40\,\%$ of them are included in the final axis lines (counting both vortices and fragments). Second, in flow fields densely populated by vortices, close encounters between axis lines of separate vortices are common: spatial proximity does not necessarily indicate connectivity. The current algorithm takes advantage of the fact that, despite the overall complexity of vortex configuration and distribution, wall-generated vortices can be traced back to the near-wall region where their legs are regularly aligned (largely in parallel) in the streamwise direction. Regularity in their distribution pattern makes tracking easier and connection is usually unambiguous. This is the rationale behind the choice of axis-line initiation in the $x$ -search. Continued propagation in other ( $y$  and $z$ ) directions minimizes false inclusion of points and false connection with other vortices by requiring the new segments to be natural extensions from the growing axis line. By contrast, a general multi-initiation and multi-directional algorithm would indiscriminately connect any points in the vicinity of a growing end. Extension of the VATIP algorithm to detached vortex structures with no preference to the $x$ direction and more complex branch configurations calls for new physics-based constraints to be incorporated, which is a focus of our future research.

Another potential challenge of extending the method to higher $Re$ is the determination of parameters. Section 4.5 thoroughly examined the effects of the adjustable parameters $H$ and $\unicode[STIX]{x1D701}$ and proposed their recommended ranges of use based on the balance between minimizing false connections and minimizing false disintegrations or truncations of vortex axis lines. For higher $Re$ , structures in the bulk regions (higher $y^{+}$ ) become non-trivial. Del Álamo et al. (Reference Del Álamo, Jiménez, Zandonade and Moser2006) showed that a lower $H$ is required to capture these detached structures because of the overall lower turbulent intensity in those regions. A $y^{+}$ -dependent $H$ was thus proposed for vortex identification in that study. It is likely that, for VATIP, a similar approach needs to be taken for both $H$ and $\unicode[STIX]{x1D701}$ . Determining the dependence of these parameters on $y^{+}$ will require trial and error. Moreover, whether a ‘sweet-spot’ range still exists for these parameters in the high- $Re$ and high- $y^{+}$ regime remains to be seen. Finally, we note that, as a brand new method, its future application and testing in broader parameter regimes and systems will be essential for its continued improvement and generalization. In this sense, the development of VATIP itself is an ‘iterative’ process that requires the experience and feedback from its application.