2. Data for analysis

The direct numericalsimulation (DNS) data for turbulent channels used in this paper are for $Re_\tau =180$ , 550, 1000 and 5200. Here, $Re_\tau =u_\tau \delta /\nu$ , with $u_\tau =\sqrt {\tau _w/\rho }$ as friction velocity, $\tau _w$ the mean wall shear stress, $\rho$ the density, $\nu$ the viscosity and $\delta$ the half-height of the channel. Our data for the first two Reynolds numbers have been validated in Xie et al. (Reference Xie, He, Bao and Chen2021). For $Re_\tau =1000$ and 5200, the data have been simulated by Lee & Moser (Reference Lee and Moser2015), available from the Johns Hopkins Turbulence Database (Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008).

Table 1 shows the grid spacings in terms of the Kolmogorov length scale $\eta =(\nu ^3 / \langle \epsilon \rangle )^{1/4}$ , where $\langle \epsilon \rangle = \langle \nu \partial _j u_i^\prime \partial _j u_i^\prime \rangle$ is the local turbulent dissipation rate. Here and elsewhere, $\langle \cdot \rangle$ denotes the ensemble average. By nominal standards, this resolution ensures that small scales are resolved. As the dissipation decreases with the wall-normal distance $y$ , $\eta$ increases from its minimum at the wall to its maximum at the centre, which is also why the grid spacings normalised by $\eta$ vary from the wall to the centre.

Table 1. Discretisation of DNS of Navier–Stokes equation. Cases of $Re_\tau =180$ and 550 are our current simulations, and $Re_\tau =1000$ and 5200 are simulations by Lee & Moser (Reference Lee and Moser2015) from Johns Hopkins Turbulence Database. Here, $N_x$ , $N_y$ and $N_z$ are the number of grids in the corresponding directions. Note that $\Delta x$ and $\Delta z$ are uniform grids, while $\Delta y$ stretches in the wall-normal direction in terms of  $\eta$ , which is the local Kolmogorov length scale.

Figure 1. (a) Sketch of the computational domain for channel where $L_x$ and $L_z$ are the streamwise and spanwise sizes, respectively. (b) An illustration of the circulation around a rectangular loop with side lengths $r_x$ and $r_z$ in a wall-parallel plane, top-left for $y^+=Re_\tau$ and bottom-right for $y^+=5$ (superscript ‘+’ denotes normalisation in wall units). The panels show contours of the normalised wall-normal vorticity $\omega _y\boldsymbol {e_y}$ in the domain of 2000 $\eta$ in length and 1000 $\eta$ in width using the DNS data at $Re_\tau =5200$ . Here, $\rm \omega_{\textit{rms}}$ is the root mean square of $\omega_y$ . Spot-like structures appear in the centre plane while streak-like structures are visible near the wall.

Fully developed turbulent channels are homogeneous in wall-parallel $x$ – $z$ planes (figure 1 a), so that circulation in those planes is the focus here. This also enables a potential comparison with the results obtained from HIT. A sketch of the circulation around a rectangular loop in the $x$ – $z$ plane is shown in figure 1(b) for $y^+=Re_\tau$ , i.e. the centre plane, and for $y^+=5$ , a plane in the viscous sublayer. The former plane features vanishing mean shear with blob-like vorticity, whereas the latter experiences strong mean shear with dominant coherent structures, such as velocity streaks (Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967). We thus expect a notable difference in circulation characteristics for planes with increasing wall-normal distance.

Experiments typically use loop integration as they provide more accurate velocity measurements than velocity gradients (Sreenivasan et al. Reference Sreenivasan, Juneja and Suri1995; Mugundhan & Thoroddsen Reference Mugundhan and Thoroddsen2023), but we calculate the circulation from DNS data via vorticity integration in (1.1). The integrated area is denoted as $A=r_x\times r_z$ and $\omega _y \boldsymbol {e}_y$ indicates the wall-normal vorticity. Both methods were attempted in Iyer et al. (Reference Iyer, Bharadwaj and Sreenivasan2021), obtaining complete equivalence except for loops with side length of $\eta$ . Square loops are first considered, i.e. $r_x=r_z$ as in figure 1(b), while the influence of loop shape will be discussed in § 3. As shown in table 1, the sampling area for averaging at $Re_\tau =5200$ spans from $60\eta ^2$ to $10^8\eta ^2$ in the centre and from $2\eta ^2$ to $10^9\eta ^2$ near the wall, thus offering a wide scale range for examining scaling properties. For $Re_\tau =180$ , the sampling area extends from $5\eta ^2$ to $10^5\eta ^2$ in the centre; although smaller compared with $Re_\tau =5200$ , it is sufficient to verify the scaling by extended self-similarity (ESS) (Benzi et al. Reference Benzi, Biferale, Struglia and Tripiccione1997), as described in § 4.

3. Scaling in the IR for different wall-normal heights

We first focus on the case of $Re_\tau =5200$ for which an IR is clearly present. To capture circulation properties at different shear levels, seven heights are selected, i.e. $y^+=5200\,(Re_\tau )$ , $3900\,(0.75Re_\tau )$ , $1300\,(0.25Re_\tau )$ , $520\,(0.1Re_\tau$ ), $70$ , 15 and 5, marked by different symbols in figures 2(a, b). In the outer region (including the so-called logarithmic layer and the core), where $y^+\gtrsim 0.1Re_\tau$ , the mean shear is negligible; while in the inner region, the flow is shear dominated (figure 2 b) and highly anisotropic. The current dimensionless parameter $S^+=S\nu /u_\tau ^2={\rm d}U^+/{\rm d}y^+$ is equivalent to the definition $S^*=Sk/\epsilon$ , which has been used, among others, most recently by Alhareth et al. (Reference Alhareth, Mugundhan, Langley and Thoroddsen2024), if one adopts the classical wall scaling for the kinetic energy to be $k \sim u_\tau ^2$ , and for the dissipation rate to be $\epsilon \sim u_\tau ^4/\nu$ . Here, $S=dU/dy$ denotes the mean shear rate, and $U$ is the streamwise mean velocity. According to the DNS data used here, we find $S^*\gt 20$ for $y^+\lt 20$ , while $S^*\lt 10$ above the buffer layer.

Figure 2. (a) Mean velocity profile and (b) mean shear distribution for $Re_\tau =5200$ , with dots highlighting the selected heights. The discrete points marked on the profiles correspond to planes on which circulation statistics are presented. Their colours carry over to figure 3. (c–f) Second-order circulation moments (or the variance of circulation) as a function of the loop area $A/\eta ^2$ , with line colours corresponding to the selected heights in (a) and (b). Inset shows the corresponding local exponents, with dash-dotted line indicating the scaling in the IR for each case.

Figure 3. Determination of power-law scaling exponents by fitting circulation moments within the IR at (a)  $y^+=Re_\tau$ and (b) $y^+=5$ . Symbols are simulation data at $Re_\tau =5200$ , and solid lines (colours corresponding to the positions depicted in figures 2 a, b) are the fitting results.

The circulation variances at these selected heights are displayed in figures 2(c–f), all exhibiting unambiguous scaling. The insets further show the local slopes, whose plateau region extends for approximately two decades of IR in the loop area. In particular for the four heights in the outer region (figure 2 c), all curves collapse on each other and closely correspond to the Kolmogorov scaling of $A^{4/3}$ , quite similar to that in HIT. This concurrence indicates that the mean shear effect is, in fact, negligible for the outer circulation. In contrast, closer to the wall, the scaling exponents are not universal, gradually decreasing with smaller $y^+$ : $A^{0.93}$ at $y^+=70$ , $A^{0.87}$ at $y^+=15$ and $A^{0.81}$ at $y^+=5$ .

Figure 4. The IR scaling exponents $\zeta _p$ as a function of moment order $p$ at different heights for $Re_\tau =5200$ . The long-dashed line is the K41 prediction, $\zeta _p=4p/3$ . The solid line indicates the unifractal model for results above the log layer, while the short-dashed lines passing through the origin indicate unifractality below the log layer. Symbols are consistent with figures 2(a,b).

By extending the findings of figure 2 to higher orders, we can assess intermittency effects under varying shear levels. High-order moments are presented in figures 3(a) and 3(b) for the two planes of $y^+=Re_\tau$ and $y^+=5$ , respectively. The best power-laws in the IR, marked by solid lines, are in close agreement with the data. We repeat this procedure for all the selected planes and collect in figure 4 the scaling exponents of circulation moments for orders $p=1{-}10$ , with statistical uncertainties (obtained using the student’s $t$ -distribution with 95 % confidence intervals) subsumed by symbol thicknesses.

Two notable observations should be made about figure 4. First, in the inner layer ( $y^+\lesssim 70$ ) the data for all moment orders are best fitted by straight lines $\zeta _p = k_{y}p$ without intercepts, where $k_y$ increases with wall distance. This property implies that circulation in the inner planes resides on space-filling unifractal sets. Not surprisingly, the slopes of unifractality $k_y$ in figure 4 are identical to the scaling exponents of their variance $\langle |\Gamma _A|^2\rangle \sim A^{k_y}$ in figure 2.

The second observation is that circulation in outer planes ( $y^+\gtrsim 0.1Re_\tau$ ) exhibits a bifractal behaviour. A single linear relation with zero intercept cannot be fitted to the data but can be approximated well by two straight lines, one for $p\lt 2$ , the other for $p \geqslant 2$ . The former is consistent with K41, while the latter can be fitted by the relation

(3.1) \begin{equation} \zeta _p = hp+(2-D), \end{equation}

where $h$ is the Hölder exponent representing the degree of singularity and $D$ represents the corresponding fractal dimension. Note that the Hölder exponent quantifies degree of local singularity of a relevant physical quantity, with the smaller $h$ signifying a stronger singularity, see e.g. Frisch (Reference Frisch1995). Impressively, both $h=1.05$ and the dimension $D=1.38$ are invariant with respect to the wall-normal height in the outer region, suggesting the universality of circulation in this region.

In summary, circulation in the outer region resides on a single bifractal set with known Hölder exponent as well as dimension, similar to HIT. On the other hand, circulation in the inner layer resides in sub-K41 unifractal sets whose dimension varies with the height of the wall. The transition to the outer bifractal behaviour appears quite gradual. The uniform bifractal behaviour in the outer region compared with the unifractality in the inner region highlights the influence on flow structures from different levels of mean shear.

Figure 5. Normalised probability density function $\mathcal {P}$ of $\Gamma _A/(u^\prime _{rms}\delta )$ in (a) inner region (with dashed lines representing the Gaussian distribution) and (b) outer region, both with a loop size of $A/\delta ^2=0.0036$ in the IRs. (c) Circulation flatness $F(A)$ at different heights as a function of loop size for $Re_\tau =5200$ ; dashed line indicates the Gaussian flatness of 3. The grey shading highlights the iR $A/\delta ^2\sim 10^{-3}$ to $10^{-1}$ . Symbols and line colours are consistent with figures 2(a)–2(b).

We now turn to the PDFs of circulation. Clearly, as shown in figure 5(a), the PDFs are essentially Gaussian in the inner region. This is also reflected in the circulation flatness $F(A)=\langle \Gamma _A^4\rangle / \langle \Gamma _A^2\rangle ^2$ , which is approximately $3$ in the scaling region represented by the shaded region in figure 5(c). Based on the Gaussian PDF, one readily has the relation between $p$ th order moments and the variance: $\langle |\Gamma _A| ^p\rangle \propto \langle |\Gamma _A| ^2\rangle ^{p/2}$ . Because $\langle |\Gamma _A|^2\rangle \sim A^{k_y}$ has been validated in figure 2, we have $\langle |\Gamma _A|^p\rangle \sim A^{ k_y p /2}$ , and thus $\zeta _p=k_y p$ . In contrast, the PDFs in the outer region depart strongly from the Gaussian distribution, with approximately stretched exponential functions (figure 5 b), so a bifractal intermittency appears. The difference between inner and outer regions could be attributed to viscous effect, as it can damp out extreme events and suppress intermittency, leading to Gaussian PDFs near the wall. Also note in figure 5(a) the non-monotonic PDFs at $y^+=5$ , $15$ and $70$ . These might be due to the non-monotonic variation of the root mean square of the streamwise velocity fluctuation $u'_{rms}$ , which has its maximum at $y^+ = 15$ . Such a maximum corresponds to prominent velocity streaks, which are induced by wall-normal motions (e.g. sweeps and ejections) of streamwise vortices, so that circulation in wall-parallel planes is less intense and hence has a narrower PDF at $y^+ = 15$ .

Figure 6. (a) Normalised probability density function of normalised circulation around closed loops with a fixed area but varying aspect ratios. Both (a) and (b) are sampling at the channel centre plane for $Re_\tau =5200$ . Solid lines correspond to loops with both sides contained within the IR; they collapse on each other. Dashed lines indicate loops with one side outside the IR, and thus depart from the collapsed curves. (b) Normalised second-order moments of circulation as a function of area. Data for loop ratio $2:1$ are represented by symbols, and for $1:1$ by a solid line. The dash-dotted line indicates the K41’s scaling $\langle \Gamma _A^2 \rangle \sim A^{4/3}$ . Inset shows the relative difference for data for the two loop ratios; it is close to unity shown by the dashed line.

The data considered so far are for square loops. We now explore the impact of the aspect ratio of rectangular loops to assess the area rule, which states that the circulation properties depend only on the loop area (Migdal Reference Migdal1994). This rule was verified conclusively by Iyer et al. (Reference Iyer, Sreenivasan and Yeung2019, Reference Iyer, Bharadwaj and Sreenivasan2021) in HIT at high Reynolds numbers, and somewhat tentatively because of the low Reynolds numbers by Cao et al. (Reference Cao, Chen and Sreenivasan1996) in HIT and by Benzi et al. (Reference Benzi, Biferale, Struglia and Tripiccione1997) in shear flows. For the present case, figure 6(a) shows the PDFs of $\Gamma _A$ for six rectangular loops with the same area $A=1908\eta ^2$ , but varying aspect ratios $r_x : r_y$ from $1:1$ to $9:1$ . The PDFs collapse well when both sides of the rectangle lie within the IR but, as expected, deviations occur when one side of the loop extends outside the IR (for loops with $r_x/\eta \gt 100$ ). Furthermore, figure 6(b) shows the circulation variance as a function of area for two different aspect ratios, that is, $r_x : r_y=1:1$ and $1:2$ . In general, the two curves collapse well with each other, particularly within the IR ( $A/\eta ^2 \lesssim 10^4$ ). Indeed, as shown in the inset of figure 6(b), the differences between the two data sets are less than 2 %. These plots provide clear evidence that circulation properties depend only on the loop area instead of its shape. This is an important conclusion as it suggests the existence of dynamical invariance with respect to variously shaped vortical structures such as hairpins, horseshoes and vortex-packets.

4. Extended self-similarity at moderate Reynolds numbers

We now investigate whether the bifractality of the circulation in the centre plane persists for lower Reynolds numbers. The circulation moments are shown in figure 7(a), with the flatness at different Reynolds numbers shown in figure 7(b), and the relative scaling exponents shown in figure 7(c). Several points are made below.

Figure 7. (a) The ESS plot of $\langle |\Gamma _A|^p\rangle ^{1/p}$ versus $\langle |\Gamma _A|^2\rangle$ at the channel centre for $Re_\tau =550$ . Symbols are data, and lines are the power-law fits. The top-left inset shows the relative difference between the fits and data; the bottom-right inset shows the centre-plane circulation variance for $Re_\tau =180$ , 550, 1000 and 5200 (from top to bottom), with the dash-dotted lines representing the K41’s scaling. (b) Circulation flatness at the channel centre for all $Re_\tau$ cases (symbols) compared with the Gaussian value $3$ (dashed line). (c) The ESS scaling exponents $\zeta _p/\zeta _2$ at the channel centre for all $Re_\tau$ . The dashed line is K41, and the solid lines are the monofractal fits. Data of three-dimensional HIT (Iyer et al. Reference Iyer, Sreenivasan and Yeung2019), quantum turbulence (Müller et al. Reference Müller, Polanco and Krstulovic2021) and quasi-two-dimensional turbulence (Zhu et al. Reference Zhu, Xie and Xia2023) are also included for comparison.

First, from the bottom-right inset of figure 7(a), there is no clear IR for $Re_\tau \leqslant 1000$ . But ESS (Benzi et al. Reference Benzi, Biferale, Struglia and Tripiccione1997) extends the scaling regime, i.e. $\langle |\Gamma _A|^p\rangle ^{1/p}$ versus $\langle |\Gamma _A|^2\rangle$ , as shown in the main graph of figure 7(a) for $Re_\tau =550$ . The extended scaling range enables a robust determination of the relative scaling exponents as the power law over several decades of the new abscissa. The least-square fit by a power law achieves errors within 3 % over four decades of $\langle |\Gamma _A|^2\rangle$ , shown in the top-left inset of figure 7(a). We have also examined the ESS scaling for $Re_\tau =5200$ , which is closely in agreement with direct measurement in the IR, thus validating the reliability of ESS.

Second, we note that the ESS of $Re_\tau =180$ and 550 are closer to the slope of K41 compared with higher $Re_\tau$ , whose explanations could be found in figure 7(b). That is, for small $Re_\tau$ , the viscous effect is stronger, and the flatness is closer to the Gaussian value $F(A)=3$ , and thus the intermittency of circulation is weaker. However, compared with the Gaussian value, there are still deviations in the IR, e.g. $F(A)\approx 4.5$ at $A/\eta ^2\sim 10^2$ , so bifractality is maintained even for $Re_\tau =180$ and 550.

Finally, the last point of the previous paragraph is expanded here. The relative scaling exponents shown in figure 7(c) depart from the K41 paradigm $\zeta _p/\zeta _2=p/2$ for high order  $p$ , where $\zeta_2$ is the scaling exponents for order 2. Similar to the linear relation of (3.1), the relative exponents for $p\gt 2$ can be fitted by the straight line

(4.1) \begin{equation} \zeta _p/\zeta _2 = [hp+(2-D)]/\zeta _2, \end{equation}

where $\zeta _2 = 8/3$ by taking K41’s scaling as a common normalisation (it does not necessarily mean that K41 is valid for low $Re_\tau$ ). For $Re_\tau =1000$ and 5200, the data collapse with the same $h=1.05$ and $D=1.38$ as shown in figure 4, indicating an asymptotical $Re_\tau$ -invariance. For $Re_\tau =180$ and 550, they are higher than those of $Re_\tau \gtrsim 1000$ with $h=1.18$ and $D=1.54$ , closer to the observation of HIT (Iyer et al. Reference Iyer, Sreenivasan and Yeung2019), quasi-two-dimensional turbulence (Zhu et al. Reference Zhu, Xie and Xia2023) and quantum turbulence (Polanco et al. Reference Polanco, Müller and Krstulovic2021; Müller & Krstulovic Reference Müller and Krstulovic2024). Therefore, the bifractality of circulation is a general property of turbulence despite differences in geometries and Reynolds numbers of turbulent structures.

Note that for lower Reynolds numbers ( $Re_\tau \leq 1000$ ), the verification of area rule is not performed here, as the IR has not been well developed. Also, in planes that are not parallel to the wall, due to the inhomogeneous and anisotropic effect, the shape and orientation of Eulerian loops may affect the statistics of circulation and deserve future studies.