2.1 Numerical method

The turbulent flow data to be used are generated from the DNS of a constant-property fluid driven through a rectangular domain that is bounded by four equal, plane, smooth, and rigid walls and two opposing flow-through sides: see figure 1. The $x$ -axis is normal to the open boundaries, and along it a pressure gradient is applied that maintains the mean flow and ensuing turbulence, which are assumed to be fully developed. Instead of inflow and outflow conditions on the velocity at the open boundaries, the usual assumption of an instantaneously periodic velocity field is made. On the four walls, no-slip and impermeability conditions are applied. The dimensions of the domain in $x$ , $y$ and $z$ are $L_{x}$ , $2h$ and $2h$ , respectively.

Figure 1.

Geometry and axes system.

The dynamical model for this flow is the Navier–Stokes equations constrained by the incompressibility of the fluid. Defining the alternate space coordinates $(x_{1},x_{2},x_{3})\equiv (x,y,z)$ and velocity components $(u_{1},u_{2},u_{3})\equiv (u,v,w)$ , the governing equations are

(2.1) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}t}+\frac{1}{2}\left[u_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{i}u_{j}}{\unicode[STIX]{x2202}x_{j}}\right]=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}u_{i}}{\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}x_{j}}-\frac{\text{d}P}{\text{d}x}\unicode[STIX]{x1D6FF}_{i1} & & \displaystyle\end{eqnarray}$$

and

(2.2) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{j}}=0, & & \displaystyle\end{eqnarray}$$

where the summation convention over repeated subscripts applies with $t$ , $p$ and $\unicode[STIX]{x1D708}$ being the time, kinematic pressure and kinematic viscosity respectively; $\text{d}P/\text{d}x$ is the imposed streamwise pressure gradient, which at every instant is uniform throughout the fluid volume; $\unicode[STIX]{x1D6FF}_{i1}$ is the Kronecker function. The skew-symmetric form of the nonlinear terms shown in (2.1) is coded for its enhanced numerical stability.

The mean velocity components will be denoted by $U$ , $V$ and $W$ and the corresponding fluctuations by $u^{\prime }$ , $v^{\prime }$ and $w^{\prime }$ . The constant bulk flow velocity will be $U_{m}$ , mean centreline velocity $U_{c}$ , and the time-averaged value of $\text{d}P/\text{d}x$ is written $\overline{\text{d}P/\text{d}x}$ . The friction velocity  $u_{\unicode[STIX]{x1D70F}}$ determined by

(2.3) $$\begin{eqnarray}\displaystyle u_{\unicode[STIX]{x1D70F}}^{2}=-\frac{h}{2}\overline{\frac{\text{d}P}{\text{d}x}} & & \displaystyle\end{eqnarray}$$

will be used as one velocity scale and in the definition $Re^{+}=hu_{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D708}$ . The bulk Reynolds number is defined here by $Re_{m}=2hU_{m}/\unicode[STIX]{x1D708}$ . Flow data scaled with $u_{\unicode[STIX]{x1D70F}}$ and $\unicode[STIX]{x1D708}$ are referred to being in wall units/scales as it is common usage, although the suitability of such a descriptive term in the square duct flow is uncertain. A local friction velocity over the duct perimeter can be calculated from the normal derivative of the local wall-parallel velocity with values ranging between zero and about $1.1u_{\unicode[STIX]{x1D70F}}$ . It has a small off- $x$ component which averages to zero over the perimeter, while the $x$ component averages to $u_{\unicode[STIX]{x1D70F}}$ to at least three significant digits. Quantities in wall units will be marked with the plus symbol in the superscript position.

All dependent flow variables are expanded in finite Fourier sums along the homogeneous direction $x$ that are substituted into (2.1) and (2.2). The orthogonality of the Fourier basis is applied to derive a set of equations for the velocity amplitudes of every streamwise mode. The nonlinear products are evaluated in real space using $3/2$ longer vectors in the streamwise direction and are truncated back to their original length after they have been retransformed into Fourier space.

Spatial approximations in the transverse directions $y$ and $z$ of the equations for each Fourier mode are based on the spectral element method of Patera (Reference Patera1984). The cross-flow area is divided into conforming quadrilateral elements and within each one the velocity and pressure variables are expanded in Lagrange polynomials over the Gauss–Lobatto–Legendre (GLL) quadrature nodes of a chosen order. With this basis, the variational form of the momentum and continuity equations are discretized using the GLL quadrature over the same set of nodes for evaluating the projection integrals. Here, the order of polynomial expansion for the velocity and pressure variables is the same – usually referred to as the $\mathbb{P}_{N}-\mathbb{P}_{N}$ method, where $\mathbb{P}_{N}$ is the space of polynomials of degree $N$ – and therefore boundary conditions for the velocity and pressure fields are necessary. On the walls, the numerical boundary conditions for the velocity are no-slip and no-through flow. The time advancement scheme uses the two-step higher-order splitting method of Leriche et al. (Reference Leriche, Perchat, Labrosse and Deville2006). In the first step, the pressure Laplacian is inverted where the nonlinear source terms and wall-normal pressure gradients are extrapolated by a second-order in time scheme. In the second step, the viscous effects are advanced with an implicit second-order Euler scheme. The overall accuracy is second-order in time and it is free from (other than the constant) spurious pressure modes. The matrices that arise from the pressure and viscous elliptic operators and that are bounded by the rectangular cross-section of the duct are efficiently inverted by diagonalization (see Lynch, Rice & Thomas Reference Lynch, Rice and Thomas1964). The zero wavenumber component of the pressure equation is used to fix the global constant pressure mode by enforcing a zero mean pressure at one corner.

The flow is driven by an adjustable streamwise pressure gradient that maintains a prescribed flow rate to an accuracy better than eleven significant digits (Deschamps Reference Deschamps1988). A laminar flow profile used in the initial conditions is linearly stable, and therefore setting off transition to turbulence requires finite amplitude perturbations. A key stage in this transition process is found to be the formation of laminar streaks – elongated regions of positive and negative streamwise velocity perturbations – present over most of the cross-section and crucially near the walls. An efficient method for attaining such a state is choosing the initial conditions to be the sum of a laminar profile at a low Reynolds number onto which several vortex rings are added whose positions, extents, orientations and strengths are varied using a pseudo-random number generator. The velocity field of each ring is calculated by taking the $\text{curl}$ of a vector potential. The Cartesian velocity components of each ring are computed from

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\times [A\exp (-a^{2}r^{2})\boldsymbol{e}_{R}], & & \displaystyle\end{eqnarray}$$

where $A$ is an amplitude parameter, $a^{-1}$ is a spatial scale, $r^{2}={r_{x}}^{2}+{r_{y}}^{2}+{r_{z}}^{2}$ the distance from the symmetry centre of the ring, and $\boldsymbol{e}_{R}$ a unit vector that defines the ring symmetry axis and thus its orientation. The composite initial field is filtered to remove sharp gradients and integrated forward in time for about five $h^{2}/\unicode[STIX]{x1D708}$ units to ensure that most of the initial vorticity has diffused to the walls. This phase is followed by successive increases of the Reynolds number with matching grid refinements until the final values of the parameters aimed at are reached. These are shown in table 1. With the final set of parameters, the simulations are advanced for a sufficiently long time until a steady mean state in the evolving global statistics is detected. A set of files containing the velocity and pressure fields at equally spaced time intervals is extracted from the later part of the steady-state segment of each simulation. During the time integration segment for storing field data, the incremental changes of the running global averages of the driving pressure gradient, turbulent kinetic energy and dissipation rate are monitored. When changes in these fall below 0.5 % over the most recent $3h/u_{\unicode[STIX]{x1D70F}}$ interval, runs are continued as needed with periodic testing of the updated velocity intensity profiles and secondary flows. The admissible type of convergence in the intensities for $Re^{+}=861$ can be seen in figure 9(a,b), where the marks are the results of averaging over the shorter interval $10.6h/u_{\unicode[STIX]{x1D70F}}$ . The intensities are well converged, especially along the diagonal. The more sensitive secondary flows, which depend on the integrated derivatives of the transverse Reynolds stresses in (6.1), are shown in figure 7(b). The secondary flow peak magnitude near the corner has converged well, while the outer flow nearer the duct centre where the secondary velocity maximum magnitude is ${\sim}0.013U_{c}$ has decreased by less than 3 % with $6h/u_{\unicode[STIX]{x1D70F}}$ or over $50\,\%$ of additional averaging. The difficulty of converging on the mean secondary flows is due to their unusual shape and relatively low magnitude; they are the outcome of averaging typical quasi-circular instantaneous vortices. This weak part of the mean flow field understates the effects of the corner on the turbulent flow which are more pronounced in other statistics such as the turbulent kinetic energy shown in figure 10. The distinctive two-minima structure is also found in separate simulations: here, when $Re^{+}=605$ and in Zhang et al. (Reference Zhang, Trias, Gorobets, Tang and Oliva2015) (figure 7 b,c). The statistics computed from each database use time and all available geometrical symmetries – when applicable – and the use is denoted with an overbar over the relevant fluctuating variables.

Table 1.

Physical and numerical characteristics of the present simulations. $U_{c}$ is the mean centreline velocity. The entries between parentheses are domain sizes of additional simulations for the same Reynolds numbers using an appropriately reduced number of Fourier modes. $N$ is the degree of the Lagrange polynomial within each element. $N_{x}$ is the number of Fourier modes in $x$ ; $N_{y}$ is the number of independent GLL nodes in $y$ ; the number of independent nodes in $z$ is equal to $N_{y}$ . The minimum $\unicode[STIX]{x0394}y^{+}$ is the normal node spacing closest to the boundaries. The cross-flow node spacing $\unicode[STIX]{x0394}z^{+}$ has an identical distribution to $\unicode[STIX]{x0394}y^{+}$ . Database time spans are in $h/u_{\unicode[STIX]{x1D70F}}$ units.

The characteristics of the simulations that are presented here are chosen around a set of constraints of acceptable node spacings. One additional limit is cited in Karniadakis & Sherwin (Reference Karniadakis and Sherwin1999, p. 283). It is that there should be least ten cosine nodes within nine wall units from each boundary in the wall-normal direction. The least compliant simulations in this respect use nine GLL nodes below the nominal nine wall units and the experience here is that polynomial order and the wall-parallel resolution must also be considered. The elements in the cross-section are generally laid out according to a hyperbolic tangent distribution and are mildly stretched so they cluster near the boundaries, but occasionally the common size of the wall-touching elements is chosen separately. The element sizes, their distribution and the polynomial order are selected for minimizing occasional reflections in the instantaneous dissipation rate peaks when they happen to straddle interelement boundaries.

In terms of the Kolmogorov length $\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D716}_{v})^{1/4}$ , where $\unicode[STIX]{x1D716}_{v}$ is the volume-averaged mean dissipation rate – see definitions in § 3 below – the maximum cross-flow grid spacing is $3\unicode[STIX]{x1D702}$ and the streamwise spacing is about $5\unicode[STIX]{x1D702}$ . The turbulent velocity spectra for the two highest Reynolds numbers in figure 2(a,b) have upturns near the cutoff wavenumbers, indicating insufficient viscous dissipation. The dissipation rate contours of figure 2(c) have small upticks on the inter-element boundaries – derivatives are calculated in the strong form – but only on those boundaries normal to the nearest wall and not those parallel to it. Taken together, the three figures suggest that intermittent steep wall-parallel gradients may require finer node spacing. The difference between the volume and sample averaged input power and the similarly averaged sum of the turbulent and mean-field dissipation rates is less than 3 %.

There is a stability restriction on the time step owing to the explicit temporal extrapolation of the nonlinear terms. Numerical experiments are conducted to establish time step limits for each Reynolds number and spatial resolution. For a fixed distribution of elements, the time step limit scales close to $1/N^{2}$ , but the eventual choice has to allow a safety margin, making it smaller than these estimates, otherwise occasionally very high velocity events can at some point destabilize the computations.

Figure 2.

(a) Streamwise spectra of $u^{\prime }$ (top set) and $v^{\prime }$ (wall-normal, lower set) velocity components at a point near the overall maximum of the dissipation rate $\unicode[STIX]{x1D716}$ : $(y,z)\simeq (0.02h,0.5h)$ ; ——,  $Re^{+}=861$ ; — —,  $Re^{+}=605$ . (b) Streamwise velocity spectra for $Re^{+}=861$ at $(y,z)=(0.516h,0.684h)$ ; ——,  $u^{\prime }$ ; —– —–,  $v^{\prime }$ ; — ⋅ —,  $w^{\prime }$ . (c) Contours of the dissipation rate $\unicode[STIX]{x1D716}^{+}$ . Above the corner bisector: $Re^{+}=605$ ; below the bisector: $Re^{+}=861$ . The round point marks the location of the spectra in (a); the square point marks spectra in (b).

The nonlinear products are evaluated in physical space at the GLL nodes and extended Fourier grid which ensures dealiased results for the homogeneous direction modes only. For $Re^{+}>350$ (approximately), the simulations may show a slow build-up of energy at the highest and nearby wavenumbers over a continually broadening interval. With high-order polynomials, simulations can be sustained for tens of thousands of time steps, but they would be increasingly tainted by the unphysical energy accumulation at high wavenumbers. A space–time filtering method based on the two-grid polynomial interpolation of Fischer & Mullen (Reference Fischer and Mullen2001) is used. Every 17 time steps, the three velocity components are interpolated onto the GLL grid for polynomial expansion one less than the simulation and then they are interpolated back onto the original grid. The $C^{0}$ continuity of the velocity field at the interelement boundaries is preserved by this procedure. The weighted averages of the original and filtered fields are used to make one first-order implicit Euler step and then proceed with the second-order scheme. The leading error term of the backward Euler step is diffusive. The filter-Euler step would normally affect the overall time-accuracy of the simulations, but the relatively infrequent use of filtering and the cautious choice of the time step mitigate this effect. Both filtering frequency and averaging weights have been chosen to minimize the energy accumulating near the high wavenumber end of the velocity spectra: see figure 2(a,b).

The simulations are run on a dedicated computer system built around a single socket motherboard (Intel LGA2011, X79 chipset.) The central processor (CPU) is the six-core Intel I7 3930K running at stock speeds – base and turbo – but the bulk of the computational work is carried on two Graphics Processing Units (GPUs) that are used as double-precision floating-point coprocessors. The GPUs, two AMD Radeon R9 280X, are mounted on two 16 lane PCI Express 3.0 (PCIe) slots with the hardware lanes built into the processor, 40, enabling the maximum data transfer rate allowed by this communication interface for both GPUs. Parallelization on the CPU side of the code is built around the OpenMP libraries and is used for launching multiple threads either for doing calculations or controlling data transfers and computing activities associated with each GPU. The code executing on the GPUs, a collection of routines called kernels, conforms to the OpenCL 1.2 standard. The evaluation of the (skew-symmetric) nonlinear terms, Fourier transforms and the sequence of matrix operations needed for finding the pressure and velocity fields are divided among the GPUs, whereas all computations that involve direct stiffness or engage data from the previous two time steps are multithreaded on the CPU. Every simulation is preceded by an optimization and code-checking phase. For the $Re^{+}=605$ simulation with 132 million nodes, the averaged elapsed time for one time step is $0.085~\unicode[STIX]{x03BC}\text{s}$ per grid point with GPU occupancy averaging 41 %, which includes data transfer times over the PCIe interface; the matrix–matrix multiplication GPU kernels in the pressure and velocity calculations, which implement a tiled product algorithm and are batched in groups of 64 wavenumbers, are executing at the sustained rate of $430\times 10^{9}$ double precision floating point operations per second per GPU. The figures corresponding to 213 million points of $Re^{+}=861$ are $0.1~\unicode[STIX]{x03BC}\text{s}$ , 53 % and $490\times 10^{9}$ floating point operations per second per GPU, respectively.

Figure 3.

Two-point velocity correlation coefficients for two Reynolds numbers: (a,b) $Re^{+}=157$ ; (c,d) $Re^{+}=861$ . Figures (a) and (c) are for nodes closest to $y=z=0.2h$ on the diagonal. Figures (b) and (d) are nearest to $(y,z)=(0.12h,0.995h)$ . For points on the diagonal, the two cross-flow coefficients are identical due to symmetry: ——,  $\unicode[STIX]{x1D619}_{11}$ ; — — —,  $\unicode[STIX]{x1D619}_{22}$ ; - - -,  $\unicode[STIX]{x1D619}_{33}$ .

2.2 Two-point correlations

One largely overlooked aspect of duct turbulence surfaced in the precursor simulations in Gavrilakis (Reference Gavrilakis1992). It was the slow decay of the two-point streamwise velocity correlations in the homogeneous direction as $L_{x}$ was increased, which was not evident in the final computation with $L_{x}=20\unicode[STIX]{x03C0}h$ . Some subsequent simulations – Huser & Biringen (Reference Huser and Biringen1993), Sharma (Reference Sharma2004), Fang et al. (Reference Fang, Yang, Wang and Bergstrom2017) – have also encountered persistent correlations of $u^{\prime }$ . But the rationale for the ‘abnormal extent’ of the space correlations in the $u^{\prime }$ field of commonly occurring flows given by Townsend (Reference Townsend1976, p. 111 and §§ 7.20–7.21), could be relevant to duct turbulence. Generalizing, Townsend suggested that fluctuations induced by wall-stress variations ‘whatever the cause’ may persist through the whole flow.

The two-point longitudinal velocity correlation coefficients are defined by

(2.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D619}_{ij}(r_{x},y,z)=\frac{\overline{u_{i}^{\prime }(x,y,z)u_{j}^{\prime }(x+r_{x},y,z)}}{\sqrt{\overline{[u_{i}^{\prime }(y,z)]^{2}}\;\;\overline{[u_{j}^{\prime }(y,z)]^{2}}}}, & & \displaystyle\end{eqnarray}$$

where $r_{x}$ is now the streamwise separation. In figure 3, the three diagonal components are shown for two transverse locations and for the highest and lowest Reynolds numbers available here. The $Re^{+}=157$ results with the longer $L_{x}$ of $16\unicode[STIX]{x03C0}h$ show that $\unicode[STIX]{x1D619}_{11}$ remains close to $-$ 7 % near the maximum separation for points on the diagonal (figure 3 a), which supports the use of $L_{x}=20\unicode[STIX]{x03C0}h$ in the earlier simulation employing a comparably low Reynolds number. Persistence in $\unicode[STIX]{x1D619}_{11}$ is found over a substantial fraction of the duct cross-section as well as a decline toward zero in other parts as in figure 3(b). Tests with shorter domains (table 1) show some effects for $Re^{+}$ below 300, but no such effects are detectable above the statistical uncertainty for $Re^{+}=605$ . These can be understood, and perhaps predicted, by referring to the structure of the $u^{\prime }$ component in the instantaneous wall turbulence. At the lowest Reynolds number, the $u^{\prime }$ distribution near the walls is dominated by a single pair of long high- and low-speed streaks (figure 4 a), from which the effects of a shorter domain may be inferred. Even though the mechanism is nonlinear, reducing $L_{x}$ will tend to hinder the lateral displacement of the streaks, thus heightening alignment with the $x$ -axis and therefore increasing $r_{x}$ to the first zero crossing of $\unicode[STIX]{x1D619}_{11}$ , whereas the intensity of $u^{\prime }$ near the wall will appear to increase or decrease depending on the length of the cut domain and position in the cross-section. These effects are petering out beyond $Re^{+}=295$ and figure 4(b) shows why this occurs: a gradual build-up of a substantial sub- $2h$ turbulent component. For example, the $Re^{+}=605$ simulations with $L_{x}=4\unicode[STIX]{x03C0}h$ and $2\unicode[STIX]{x03C0}h$ show a difference in the peak values of $\overline{{u^{\prime }}^{2}}/u_{\unicode[STIX]{x1D70F}}^{2}$ on the wall bisector of about 1 %; assuming that such a trend would continue, a smaller effect would be expected when $Re^{+}=861$ . Overall, for the simulations of table 1 and the statistics used herein, no $L_{x}$ effects are detectable above the sampling uncertainty, and therefore all simulations including the $Re^{+}=861$ (figure 3 c,d) are considered reliable approximations to the decorrelated cases.

Figure 4.

Instantaneous turbulent velocity field component ${u^{\prime }}^{+}$ over a plane near one of the four walls. (a) $Re^{+}=157$ at $y^{+}=16.3$ . The $z$ -axis has been expanded by a factor of 2. Velocity data span the interval $[-11.7,10.8]$ . Greyscale limits: black, minimum; white, maximum. The dashed line traces a high-speed region; it is at an angle of $2.2^{\circ }$ to the $x$ -axis. (b) $Re^{+}=861$ at $y^{+}=22.6$ . True aspect ratio. Velocity range $[-8.3,8.2]$ .

We must still account for $\unicode[STIX]{x1D619}_{11}$ being negative at larger separations, which is also found in Huser & Biringen (Reference Huser and Biringen1993) and Sharma (Reference Sharma2004). Figure 4(a) appears helpful at illustrating the flow underlying a long wavelength oscillatory correlation. High- and low-speed regions are seen extending from wall to wall and exceeding $L_{x}/2$ in downstream length, but at a small angle to the $x$ -axis. The data of figure 4(a) show a streak angle with the $x$ -axis of $2.2^{\circ }$ , and therefore the lateral displacement of the high-speed pattern over the duct half-length will be $z^{+}\sim 150$ . Since the typical separation of adjacent high- and low-speed streaks in boundary layers and plane channels is known to be about 50 wall units, a pair of points aligned on the $x$ -axis are more likely to fall on opposite-signed streaks at (near) maximum $r_{x}$ . For $\unicode[STIX]{x1D619}_{11}$ to remain wholly negative beyond the first zero crossing (figure 3 a), a shallower structure angle must be prevalent at that position. It has been assumed that the streak lateral separation scaling holds for ducts as well and all simulations shown here appear to have this streak spacing. Similarly, the $u^{\prime }$ distribution for $Re^{+}=861$ over planes parallel to the walls (figure 4 b) now has oblique sequences of structures stretching over distances comparable to the duct’s half-length. The angles with the $x$ -axis could be as high as $7^{\circ }$ , but for the instance shown, they are in the opposite sense to that in figure 4(a).

6.1 Reynolds stress form

Like all other duct statistics, the mean vorticity and associated correlations will adapt to the symmetries of the bounding geometry whose basic constituent is an octant. For every pair of adjacent octants, $\unicode[STIX]{x1D6FA}_{x}$ and every term in (6.1) will be antisymmetric with respect to their common boundary on which their values must be zero. In the $\{z\leqslant h;y\leqslant z\}$ – or lower left – octant the vorticity at the centre of the large vortex, as in figure 7(a), will be negative whereas near and at the wall it is positive and on the two internal symmetry lines zero (figure 12 a upper octant but with a reversed sign). The magnitude of the overall maximum of $\unicode[STIX]{x1D6FA}_{x}$ which occurs on the walls, is more than 25 times that at the centre of the secondary vortex in figure 7(a). It is significant that high values of $\unicode[STIX]{x1D6FA}_{x}$ extend over most of the duct wall, which indicates a distributed generation mechanism whose origin must be traced to the effects of the corners. In terms of the mean velocity field, it consists of generating a net velocity component close and parallel to the no-slip wall, which drives the flow within a triangle on whose two internal boundaries free-slip conditions apply.

The distributions of the two turbulent terms over part of the cross-section for $Re^{+}=861$ and $605$ are shown in figure 11. It needs to be remembered here that spectral element computations do not enforce continuity of higher-order velocity derivatives at the inter-element boundaries, but they converge toward such continuity with increasing resolution. With restricted spatial resolution and depending on the local state of the flow, boundary values of such derivatives from one element may differ from those calculated from the adjoining element since each covers a different part of the flow domain. Here, the common-boundary values for the expressions of the turbulent source in (6.1) and $T_{1}$ in definitions (6.2) below are assigned the estimate from whichever of the adjoining elements is closest to the average of their boundary-adjacent nodes.

Figure 12.

Mean streamwise vorticity $\unicode[STIX]{x1D6FA}_{x}$ and the source terms of its equation in the vorticity form (6.2) for $Re^{+}=861$ . (a) Above the diagonal, $h\unicode[STIX]{x1D6FA}_{x}/U_{m}$ ; below, turbulent vorticity shearing $T_{3}$ . (b) Above the diagonal, turbulent transport $T_{1}$ ; below, intensification of $\unicode[STIX]{x1D714}_{x}^{\prime }$ , $T_{2}$ . The $T$ -terms are scaled with the maximum of $|T_{2}|$ . The interval between adjacent $T$ contours is $0.2$ ; the maximum positive contour is $0.8$ . Mean vorticity intervals $0.2U_{m}/h$ with range [ $-3.4,1.4$ ]. The additional contour marked 0.1 is 37 times smaller than the absolute maximum on the wall. Zero contour lines are not shown. Positive values are drawn with solid lines, negative with dotted. All four quantities have the correct sign for the octant they are drawn in; they would change sign on being reflected about the diagonal. The filled symbols show the centre of the secondary flow.

For both $Re^{+}$ , the terms containing the Reynolds stresses attain their peak (absolute) values very close to the corner with marked overlap. The maximum value of the anisotropy term exceeds that of the shear stress term by about 6 % and they can be assumed to both be mostly related to the near-corner part of the secondary flow. In the same octant, the two terms are largely of opposite sign, resulting in substantial cancellation which tends to lessen with increasing Reynolds number. Due to imperfect overlap, the maximum of their sum when $Re^{+}=861$ is about half the peak value of either. This would suggest that the two turbulent terms in (6.1) are unlikely to originate in independent flow processes at least near their peak values. Firm analytical support for this is embedded in the derivation of the ‘stretching’, $T_{2}$ , and ‘tilting’, $T_{3}$ , terms that are used in the next subsection; each one is formed from products of velocity derivatives drawn from both turbulent terms of (6.1). These are statistics of the nominal generation process associated with flow events such as those shown in figure 8(c,d). Farther from the corner, $z>0.3h$ (not shown), the turbulent terms have and add up to low values, leaving their dynamical significance unclear.

Averaging is a linear operation, and therefore the unbarred form of the turbulent terms may be used to locate their most intense events within the instantaneous fields. Prominent instances of the turbulent terms in (6.1) within the evolving velocity field are outlined here with an alternative view in terms of the turbulent vorticity to follow. The description is based on the $Re^{+}=861$ data in order to minimize the inclusion of any lingering very low Reynolds number effects in addition to having the largest available inner–outer cross-flow separation. High instantaneous (absolute) values of both turbulent terms are more than an order of magnitude larger than their respective maximum or minimum mean values drawn in figure 11. Unusually high values are rarely seen near the corners despite their averages having their extrema there. They are found mainly closer to the rim of vortices rather than their centres, in shear layers turning under the blocking influence of the walls and on wall-parallel shear layers when detaching. The common traits of the instantaneous velocity field that underlie the turbulent terms in (6.1) are significant instantaneous streamline curvature and large values of $v^{\prime }$ and $w^{\prime }$ , and therefore, where the cross-plane velocity vector is inclined about $45^{\circ }$ to the nearby boundary. These can occur over a broad strip around the duct perimeter. High instantaneous values for these terms have not been observed near the duct centre. They are particularly marked where vortices are generated at the walls including near the corners. The Reynolds stress form of the $\unicode[STIX]{x1D6FA}_{x}$ equation emphasizes statistical aspects of the turbulent velocity field rather than precise dynamical processes, which are difficult to interpret as its two generating terms overlap over the same flow events and which is reflected on their averaged distribution. For the two highest Reynolds numbers, their maxima/minima occur around 26 wall units from the nearest side wall and 3.2 wall units above the nearest ‘horizontal’ one. They are related to vortices whose downflow almost touches one corner wall and impinges on the other as in figure 8(c,d). Both the well-defined position and scaling of these suggest flow instabilities specific to the mean profile close to the corner.

6.2 Turbulent vorticity form

A more incisive way of obtaining the information contained in the terms due to turbulence of (6.1) is to replace their sum by an expression of three terms coupling the turbulent vorticity components $\unicode[STIX]{x1D714}_{x}^{\prime }=\unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}y-\unicode[STIX]{x2202}v^{\prime }/\unicode[STIX]{x2202}z$ , $\unicode[STIX]{x1D714}_{y}^{\prime }=\unicode[STIX]{x2202}u^{\prime }/\unicode[STIX]{x2202}z-\unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}x$ , $\unicode[STIX]{x1D714}_{z}^{\prime }=\unicode[STIX]{x2202}v^{\prime }/\unicode[STIX]{x2202}x-\unicode[STIX]{x2202}u^{\prime }/\unicode[STIX]{x2202}y$ and fluctuating velocity. Using the continuity equation and flow homogeneity in $x$ , the Reynolds stress terms can be replaced by the sum of

(6.2a-c ) $$\begin{eqnarray}\displaystyle T_{1}=-\left(\overline{v^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{x}^{\prime }}{\unicode[STIX]{x2202}y}}+\overline{w^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{x}^{\prime }}{\unicode[STIX]{x2202}z}}\right),\quad T_{2}=2\overline{\unicode[STIX]{x1D714}_{x}^{\prime }\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}},\quad T_{3}=\left(\overline{\unicode[STIX]{x1D714}_{y}^{\prime }\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}y}}+\overline{\unicode[STIX]{x1D714}_{z}^{\prime }\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}z}}\right)\!.\quad & & \displaystyle\end{eqnarray}$$

Unlike the Reynolds stress form, expressions (6.2) represent three distinct vorticity processes where the first two terms show how $\unicode[STIX]{x1D714}_{x}^{\prime }$ can be added to (or subtracted from) $\unicode[STIX]{x1D6FA}_{x}$ : turbulent transport in the presence of gradients and vorticity intensification. The last term, $T_{3}$ , is of the type that characterizes the origin of Prandtl’s secondary flows of the first kind representing the shearing of the two transverse vorticity components. The first term, $T_{1}$ , is in a sense passive; it transports pre-existing vorticity when gradients are present and it would be significant at the point where the generation of $\unicode[STIX]{x1D714}_{x}^{\prime }$ is intense and highly localized in the absence of other competing processes. It appears mostly as a source close to the corner and as a sink term farther toward the mid-wall (figure 12 b upper octant). The term $T_{2}$ , whose dynamically concise and precise form depends on the incompressibility (2.2), is the main source along the smooth duct wall and attains high values close to the boundaries where both direct generation of $\unicode[STIX]{x1D714}_{x}^{\prime }$ and intense normal strains from the deflection of high-speed downflows occur.

To assess the relative contribution of each term, the distribution of $\unicode[STIX]{x1D6FA}_{x}$ over the same area has to be considered. This is drawn in figure 12(a), upper octant, keeping in mind that it and the $T$ -terms are antisymmetric about the diagonal. All comparisons that follow are based on the signs appropriate over the lower octant, where the mean vorticity has a positive maximum at the wall and the secondary flow vortex sense of rotation corresponds to negative vorticity. Numerically, the near-wall vorticity is an order of magnitude larger than the value at the centre of the secondary vortex. Unsurprisingly, the $T$ -terms attain their highest values where $\unicode[STIX]{x1D6FA}_{x}$ is large and away from the corner it is the region of high positive values close to the wall. From figure 12, $T_{1}$ and $T_{3}$ are mainly sink terms and $T_{2}$ is numerically the largest and effectively the only source. It is the sum of the three terms over the region $\unicode[STIX]{x1D6FA}_{x}>0$ that underlies the outer secondary flow, but now the dynamical process of mean vorticity generation has been identified. The highest overall $T_{2}$ value is found close to the corner and it marks the position where it is most likely that the vortices generated near the corner will produce significant positive contributions to it and reflects the more organized state of this part of the flow. It is now possible to directly associate the instantaneous flow field in figure 8(c,d) with the maximum and minimum of $T_{2}$ . The corner extrema are found on the ‘downflow’ side of the larger corner vortex in figure 8(c), once where it impinges on the side wall generating the new small vortex and a second time while the vortex flow is turning due to blockage of the near wall; they re-emphasize the importance of streamline curvature noted earlier. These types of events are expected to be contributing to the generation processes father away from the corner, with the difference that the likelihood of seeing vortices of one sign or the other evens out. These findings are applicable to square ducts with smooth walls and most likely to rectangular ducts of other aspect ratios. However, it is possible that $T_{3}$ will appear as a contributing source of secondary flow vorticity within turbulence near corners of angles other than $90^{\circ }$ with or without graded wall roughness.

The question that remains is how to account for the influence of the corner geometry propagating toward the mid-wall giving rise to the elongated $T_{2}$ contours for $z>0.3h$ with very mild tapering, but abrupt end toward the mid-wall. In addition, there is a noticeable gap around $z=0.25h$ in the $T_{2}$ contours, and directly above it on the diagonal is the area of the maximum between the two negative velocity troughs in figure 7(b,c). This gap could conceivably be attributed to statistical cancellation between near- and far-corner events, but there is no symmetry point or line nearby and similar gaps are found in the remaining $T$ -terms, each representing different vorticity processes. A plausible explanation for both features is suggested by figure 4(a,b), where long oblique structures span a large fraction of the wall. Based on this observation, the twofold structure of the secondary flows for the two highest $Re$ may be constructed as the composite of a near-corner component driven by a mean-field instability producing short-lived relatively localized vortices, and of another with very long and long-living structures most evident in the form of oblique near-wall streaks of $u^{\prime }$ originating well upstream and possibly reaching the opposite side wall for the current $Re$ . The former is well supported by the statistics and instantaneous views of the velocity field while the latter is more of a conjecture that fits the available data well. Within the framework of this explanation, the slight tapering of the $T_{2}$ contours toward the wall centre can be accounted for by symmetry considerations. The mean vorticity generation due to one corner would normally extend beyond the mid-wall, but it encounters a similar but opposite-signed process from the other corner and an antisymmetry constraint about the wall bisector on which $T_{2}=0$ .

The near-corner secondary velocity component has been linked directly to the corner mean velocity profile instability and therefore does not fit well into Prandtl’s twofold classification of secondary flows. It is the dominant form for $Re^{+}$ below approximately 300, but with increasing Reynolds number it becomes confined to ever smaller regions around each corner, with its maximum seemingly at a constant distance from its walls in $\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}}$ units. The confined state of this secondary flow component for the highest $Re$ is also evident in the isolation and compactness of the $T_{2}$ contours close to the corner in figure 12(b).

Sufficiently near the $y$ -wall $\unicode[STIX]{x1D714}_{x}^{\prime }\sim \unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}y$ . The two remaining factors in $T_{2}$ are then also contributors to the turbulent dissipation rate, and therefore, away from the corner, the mean vorticity generation may be viewed as being due to a subset of the dissipative motions of near-wall turbulence. Furthermore, the instability-driven part of the secondary flow extracts its energy from the mean flow field, whereas for the outer part, the proximate source of its energy is turbulent motions, and therefore this process may be termed negative turbulence production.