Turbulent transport of momentum and passive scalar around a rectangular cylinder

Abstract

The turbulent transport of momentum, energy and passive scalar is investigated in the flow around a rectangular cylinder of aspect ratio 5 : 1 – a geometry representative of separating and reattaching flows from sharp-edged bodies. The study is based on direct numerical simulation (DNS) conducted at Reynolds numbers up to ${\textit{Re}} = 14\,000$, based on the cylinder thickness, with Schmidt number fixed at ${\textit{Sc}} = 0.71$. At this Reynolds number, the flow exhibits features of asymptotic high-${\textit{Re}}$ behaviour. Budgets of mean momentum, Reynolds stresses, mean scalar and scalar fluxes provide a detailed view of the underlying transport mechanisms. The mean momentum balance elucidates the role of turbulence in entraining free stream fluid, promoting shear-layer reattachment, sustaining backflow in the recirculation region and regulating wake dynamics through large-scale vortex shedding. The leading-edge shear layer is the main site of turbulence production, with energy injected into streamwise fluctuations and redistributed to cross-flow components by pressure–strain interactions. As ${\textit{Re}}$ increases, vertical fluctuations increasingly return energy to the mean upward flow, stabilising the separation bubble height. Turbulent transport dominates scalar redistribution. Scalar fluxes are primarily generated by interactions between Reynolds stresses and scalar gradient, and modulated by pressure-scalar gradient effects. An a priori evaluation of eddy-viscosity and diffusivity models quantifies the misalignment between modelled and DNS-resolved stress and flux tensors, as well as the inhomogeneity of eddy transport coefficients. This analysis deepens the understanding of transport phenomena in bluff-body flows approaching the asymptotic regime, and underpins the validation and improvement of turbulence models for separating and reattaching flows.

1. Introduction

Flow separation and reattachment substantially affect the transport of mass, momentum and energy (both kinetic and thermal) around bluff bodies, with implications for performance, efficiency and stability in numerous engineering systems. Examples include the aerodynamic drag of ground vehicles, wind loading on buildings, vortex-induced vibrations in structures such as bridges and chimneys, fuel–air mixing in combustors, and pollutant dispersion in urban and natural environments. Despite their importance and considerable advances in recent years, many aspects of separated and reattaching flows remain poorly understood due to their inherent complexity. Turbulent processes, such as momentum redistribution, energy exchange between mean and fluctuating motions, and heat transfer across separation and reattachment regions, are difficult to measure experimentally and predict computationally.

These difficulties stem from the unsteady, multi-scale and three-dimensional nature of such flows. Current turbulence models often fail in capturing the intricate interactions between separation, reattachment and turbulent fluctuations. On the other side, minor disturbances in experimental set-ups can significantly alter separation and reattachment dynamics, hindering data repeatability and accuracy. Overcoming these obstacles requires high-fidelity simulations, advanced experimental techniques and innovative modelling approaches capable of elucidating the underlying mechanisms governing these flows.

Separating and reattaching flows have been extensively studied using different nominally two-dimensional configurations. Sharp-edged geometries, where the separation location is fixed and independent of the Reynolds number, are particularly useful for isolating the mechanisms of separation and reattachment without the added complexity of time-dependent separation points. Commonly studied geometries include forward- and backward-facing steps (Armaly et al. 1983; Adams & Johnston 1988b ; Hattori & Nagano 2010), sudden expansions in plane channels (Durst, Pereira & Tropea 1993; Casarsa & Giannattasio 2008), blunt flat plates (Kiya & Sasaki 1983; Cherry, Hillier & Latour 1984), surface-mounted obstacles like ribs and fences (Castro 1981; Fang, Tachie & Dow 2022), and bluff bodies with rectangular cross-sections of various aspect ratios (Nakamura, Ohya & Tsuruta 1991; Moore, Letchford & Amitay 2019). These configurations share fundamental features while offering distinct insights into separation and reattachment phenomena. Comprehensive reviews of these geometries have been provided by Simpson (1989) and Ota (2000).

This study investigates the flow around a rectangular cylinder of chord-to-thickness ratio equal to 5 : 1 in a uniform stream. This geometry provides a simple yet representative model for examining the fluid dynamic characteristics of real-world separating and reattaching flows. An overview of the main features of this configuration, along with a discussion of separated flows, is provided in § 3.

The turbulent flow around the rectangular 5 : 1 cylinder has gained significant attention as part of the BARC (Benchmark on the Aerodynamics of a Rectangular 5 : 1 Cylinder) initiative (Bruno, Salvetti & Ricciardelli 2014; Mannini et al. 2018; Bartoli et al. 2020; Chiarini & Quadrio 2022; Corsini et al. 2022; Crivellini et al. 2022; Mariotti, Lunghi & Salvetti 2024). Established in 2008, BARC provides a collaborative framework for comparing numerical simulations and wind tunnel experiments. The initiative seeks to establish best practices, ensure consistency across research methodologies, and promote the integration of experimental and computational results to improve the understanding of separating and reattaching flows. In this context, direct numerical simulation (DNS) studies have proven crucial in providing high-fidelity insights into the aerodynamic characteristics and separation-reattachment mechanisms. Recent DNS studies (Cimarelli et al. 2018a ; Chiarini & Quadrio 2021, 2022; Corsini et al. 2022; Cimarelli, Corsini & Stalio 2024; Li et al. 2024b ) have started to examine this flow configuration over a range of Reynolds numbers representative of turbulent regimes, spanning $1000 \lt {\textit{Re}} \lt 14\,000$ , where ${\textit{Re}} = U_0 D / \nu$ is based on the cylinder thickness $D$ and the free stream velocity $U_0$ .

Currently, features of the turbulent transport around the rectangular 5 : 1 cylinder are well documented by various authors for Reynolds numbers up to ${\textit{Re}} = 3000$ , with results derived from DNS employing different numerical methodologies. At ${\textit{Re}} = 3000$ , Cimarelli et al. (2018a ) and Cimarelli et al. (2018b ) presented the distributions of Reynolds stresses (turbulent momentum fluxes) and turbulent kinetic energy around the rectangular 5 : 1 cylinder. Based on the same DNS dataset, Cimarelli et al. (2019) investigated turbulence production statistics, revealing the occurrence of negative production phenomena in the shear layer, where energy is transferred from turbulent fluctuations back to the mean flow. For the same Reynolds number, Chiarini & Quadrio (2021) conducted a detailed analysis of Reynolds-stress budget terms, while Chiarini et al. (2022a ) explored the role of flow structures in producing and redistributing large- and small-scale velocity fluctuations using the anisotropic generalised Kolmogorov equations. At a lower Reynolds number ( ${\textit{Re}} = 1000$ ), Li et al. (2024a ) recently examined the connection between transport mechanisms and vortical structures, extending their analysis to include rectangular cylinders of different aspect ratios.

While several experimental studies have investigated the flow around the rectangular 5 : 1 cylinder at high Reynolds numbers, detailed analyses of turbulent momentum and energy transfer for ${\textit{Re}} \gt 10^4$ remain limited. The works of Lander et al. (2018) and Moore et al. (2019) examined the dynamics of the separated shear layer over the range $13\,400 \leqslant {\textit{Re}} \leqslant 118\,000$ , showing that turbulent kinetic energy production increases with Reynolds number due to intensified instability interactions and spatial amplification of disturbances. Recently, Kumahor & Tachie (2022) employed particle image velocimetry (PIV) at ${\textit{Re}} = 16\,200$ to quantify in-plane turbulent fluctuations and the terms of the turbulent kinetic energy budget, providing insight into Reynolds stress transport by turbulence.

The present study aims to bridge the gap between low- ${\textit{Re}}$ DNS results and higher- ${\textit{Re}}$ experiments by analysing turbulent transport mechanisms up to ${\textit{Re}} = 14\,000$ . At sufficiently high Reynolds numbers, the fluxes of mass, momentum and energy in turbulent flows are expected to become independent of molecular viscosity and diffusivity (Frisch 1995). Based on the same dataset employed here, Cimarelli et al. (2024) recently showed that several integral quantities around the 5 : 1 cylinder exhibit only a weak dependence on Reynolds number for ${\textit{Re}} \gt 10^4$ , suggesting that the flow is approaching an asymptotic high- ${\textit{Re}}$ regime. Accordingly, while results at ${\textit{Re}} = 3000$ and $8000$ indicate how turbulence evolves with Reynolds number, the analysis at ${\textit{Re}} = 14\,000$ can be considered representative of the near-asymptotic behaviour. As compared with the cited work by Cimarelli et al. (2024), which focuses on mean flow topology, unsteadiness and entrainment processes, here we examine momentum and energy transfer through the full budgets of the mean momentum and Reynolds stresses at ${\textit{Re}} = 14\,000$ , and analyse passive scalar transport with comparable detail, including the budgets of turbulent scalar fluxes. These results elucidate the role of turbulence in governing shear-layer reattachment, sustaining the backflow region and modulating wake dynamics. Furthermore, this analysis enables an a priori assessment of turbulence models based on eddy viscosity and eddy diffusivity, highlighting their limitations in separated and reattaching flows, characterised by strong anisotropy and non-equilibrium effects.

Studying scalar transport gives additional insight into how turbulence affects mixing and redistribution phenomena. In applied contexts, it provides a clearer view of how turbulence spreads quantities like heat or concentration in separating and reattaching flows. In such flows, scalar transport is significantly affected by the structure and dynamics of the separation bubble, reattachment region and wake. Although scalar transport has been studied across different geometries (see the review by Ota 2000), investigations specific to the rectangular 5 : 1 cylinder are missing. The only recent contribution in this regard is by Cimarelli et al. (2024), who reported mean and variance of the scalar field. Earlier studies on blunt plates (Ota & Kon 1974; Cooper, Sheridan & Flood 1986) described typical scalar transfer distributions characterised by a local minimum downstream of separation, followed by a peak near the mean reattachment location. However, the relation between maximum wall scalar transport and reattachment dynamics is still unclear (Sparrow, Kang & Chuck 1987; Yanaoka, Yoshikawa & Ota 2003).

To summarise, this paper investigates turbulent momentum and scalar transport in the flow around a rectangular 5 : 1 bluff body by means of DNS at Reynolds numbers up to $14\,000$ and fixed Schmidt number, ${\textit{Sc}} = 0.71$ . The main objectives are to: (1) analyse the full budgets of mean momentum and Reynolds stresses to characterise turbulent transport mechanisms in separating and reattaching flows at ${\textit{Re}} \sim 10^4$ ; (2) investigate for the first time passive scalar transport in this configuration, including turbulent scalar fluxes and their budgets; (3) perform an a priori assessment of eddy viscosity and diffusivity-based turbulence models, identifying and quantifying their limitations in highly anisotropic, non-equilibrium regions.

The remainder of the paper is structured as follows. Section 2 outlines the numerical methodology and details of the simulations. Section 3 provides an overview of the key flow features around the rectangular cylinder. The turbulent transport of momentum and its influence on flow reattachment are analysed in § 4 through the examination of mean velocity fields, Reynolds stress distributions and their respective budgets. Section 5 extends the analysis to the transport and mixing of a passive scalar, presenting, for the first time in this flow configuration, the distributions of turbulent scalar fluxes, the mean scalar balance and budgets of scalar fluxes. The performance of turbulence models based on eddy viscosity and eddy diffusivity in predicting separating and reattaching flows is assessed in § 6. Finally, § 7 concludes the paper.

2. Numerical method and flow settings

The present study uses DNS datasets previously introduced by Cimarelli et al. (2024). For clarity and completeness, the numerical methodology and flow parameters used in those simulations are summarised here.

The governing equations solved are the continuity and Navier–Stokes equations for incompressible flows, and the scalar transport equation,

(2.1) \begin{equation} \left . \begin{array}{ll} \displaystyle \frac {\partial u_i}{\partial x_i} = 0, \\ \displaystyle \frac {\partial u_i}{\partial t} + u_j \frac {\partial u_i}{\partial x_j} = \displaystyle -\frac {\partial p}{\partial x_i} + \frac {1}{{\textit{Re}}} \frac {\partial ^2 u_i}{\partial x_j \partial x_j}, \\ \displaystyle \frac {\partial \theta }{\partial t} + u_i \frac {\partial \theta }{\partial x_i} = \displaystyle \frac {1}{{\textit{Re}}{\textit{Pr}}} \frac {\partial ^2 \theta }{\partial x_i \partial x_i}, \end{array}\right \} \end{equation}

where the subscripts $i$ , $j$ , $k$ take values $1$ , $2$ , $3$ to denote the streamwise, vertical and spanwise directions, with $(x_1, x_2, x_3) = (x, y, z)$ . Here, $u_i$ is the $i$ th velocity component, $(u_1, u_2, u_3) = (u, v, w)$ , $p$ is the pressure and $\theta$ represents the passive scalar field. All variables in (2.1) and throughout the paper are non-dimensionalised by the thickness of the rectangular cylinder $D$ , the free stream velocity $U_0$ , and the difference in scalar concentration between the body surface and the free stream, $\varTheta _w - \varTheta _0$ . The Reynolds number is defined as ${\textit{Re}} = U_0 D / \nu$ and the Schmidt number as ${\textit{Sc}} = \nu / \alpha$ , where $\nu$ is the kinematic viscosity and $\alpha$ the molecular diffusivity of the scalar. Note that for heat transfer with constant thermophysical properties, $\textit {Sc}$ corresponds to the Prandtl number, $\textit {Pr}$ . The velocity, pressure and scalar fields are decomposed into mean and fluctuating components as $u_i = U_i + u_{i}^{\prime }$ , $p = P + p^{\prime }$ and $\theta = \varTheta + \theta ^{\prime }$ .

DNS was performed using the code Nek5000 (Fischer, Lottes & Kerkemeier 2008), based on the high-order spectral element method for spatial discretisation. Within each spectral element, velocity and passive scalar are represented as $N$ th-order Lagrange interpolation polynomials on the Gauss–Lobatto–Legendre quadrature points, ensuring $C^0$ continuity across element interfaces. The pressure is represented by polynomials of order $N-2$ based upon the Gauss–Legendre quadrature points. In this study, seventh-order accuracy in space is achieved with $N = 7$ . Temporal discretisation uses an operator-splitting approach where viscous terms are treated implicitly via a third-order backward differentiation scheme, while nonlinear advective terms are treated explicitly using a third-order extrapolation scheme. In the DNS at the highest ${\textit{Re}}$ , the second-order accurate scheme is applied. Aliasing errors are prevented through over-integration of the advective terms by a $3/2$ factor. A low-pass explicit filter with a cut-off mode of $N-1$ and weight $0.02$ is employed for additional stabilisation (Fischer & Mullen 2001).

At the inlet, Dirichlet conditions are imposed, prescribing a uniform velocity field $(U_0, 0, 0)$ and a uniform scalar concentration $\varTheta _0 = 0$ . At the outlet and along the vertical directions, natural boundary conditions are applied, corresponding to $-p \hat {n}_j + 1/{\textit{Re}} (\partial u_i / \partial x_j) \hat {n}_j = 0$ and $(\partial \theta / \partial x_j) \hat {n}_j = 0$ , where $\hat {n}_j$ denotes the components of the unit vector normal to the boundary surface. On the body surfaces, no-slip and fixed scalar conditions, $\varTheta _w = 1$ , are enforced. Periodicity is imposed in the spanwise direction.

The computational domain size is $L_x \times L_y \times L_z = 80D \times 31D \times 5D$ . The rectangular 5 : 1 cylinder is located $20D$ downstream of the inflow boundary, is centred in the vertical direction and spans the domain width. The coordinate origin is positioned at the upper leading-edge corner, centred spanwise, as shown in figure 2. Simulations are performed for Reynolds numbers ${\textit{Re}} = 3000$ , $8000$ and $14\,000$ , with a constant Schmidt number of ${\textit{Sc}} = 0.71$ , representative of air when the passive scalar being transported is heat.

The grid is designed to resolve the smallest turbulent length scales across critical regions of the flow. For ${\textit{Re}} = 14\,000$ , the minimum grid spacing is $(\Delta x_{\textit{min}}, \Delta y_{\textit{min}}, \Delta z) = (0.0021, 0.0021, 0.007)$ at the leading edge. The maximum ratios of grid spacing to Kolmogorov scale $\eta$ in the shear layer are $(\Delta x/\eta , \Delta y/\eta , \Delta z/\eta )_{\textit{max}} = (4.2, 4.6, 6.3)$ . At the wall, the maximum grid spacings in viscous units are $(\Delta x^{+}, \Delta y^{+}_{w}, \Delta z^{+})_{\textit{max}} = (4.1, 0.66, 5.1)$ . Further details regarding the spatial resolution issue can be found from Corsini et al. (2022). The time step is fixed to satisfy the Courant–Friedrichs–Lewy condition ${\textrm{CFL}} \lt 0.5$ , with $\Delta t = 5.5 \times 10^{-4}$ , $3.0 \times 10^{-4}$ and $2.6 \times 10^{-4}$ for ${\textit{Re}} = 3000$ , $8000$ and $14\,000$ , respectively, ensuring $\Delta t$ is significantly smaller than the smallest Kolmogorov time scale, e.g. $\tau _\eta \geqslant 0.007$ at ${Re} = 14\,000$ .

For the statistical analysis, once a statistically stationary state is reached, the flow is sampled over at least $250 D/U_0$ time units, with data collected at intervals of $5 D/U_0$ . This yields a minimum of $51$ three-dimensional samples per DNS. Statistics are then computed by averaging over time, along spanwise homogeneous direction, and by exploiting the flow symmetry about the $xz$ plane at $y = -0.5$ . Consequently, the figures presented in the results only display the upper half of the rectangular cross-section, i.e. for $y \gt -0.5$ . The validations of the present datasets, along with an assessment of their quality in terms of both spatial discretization and statistical convergence, are provided in Corsini et al. (2022) for the case at $Re=3000$ and in Cimarelli et al. (2024) for $Re=14000$ .

For the sake of reproducibility and to limit the number of governing parameters, the present simulations do not include imposed inflow disturbances. However, it is worth noting that previous studies on rectangular cylinders with aspect ratio 5 : 1 and ${\textit{Re}} \gt 10^4$ (Ricci et al. 2017; Mannini et al. 2018) have shown that free stream turbulence can affect the flow. Increased turbulence intensity in the incoming flow typically leads to earlier onset of shear-layer instabilities, enhanced entrainment and reduced reattachment length. These effects also depend on the turbulence scale: small-scale disturbances, comparable to the shear-layer thickness, promote shear-layer growth and reattachment, whereas turbulence at scales of the order of the body size tends to weaken vortex shedding coherence (Nakamura & Ozono 1987). Comparable behaviours have been observed for separating and reattaching flows around rounded-edge bluff bodies and streamlined geometries, such as aerofoils and turbine blades (Bearman & Morel 1983; Öztürk & Schobeiri 2006; Wang et al. 2014). Although these configurations differ substantially from the present sharp-edged body in smooth inflow, they provide useful insight into how inflow conditions and surface curvature affect flow topology.

3. General flow characteristics

This section provides an overview of the flow around a rectangular cylinder with an aspect ratio of 5 : 1, shown in figure 1, and highlights the primary phenomena and their interrelations. The distinguishing time-periodic flow features are shortly described in § 3.1.

Figure 1. Instantaneous flow visualisations for the case at ${\textit{Re}}=14\,000$ . (a) Isosurfaces of $\lambda _2 = -5$ coloured by streamwise velocity. (b) Volume rendering of the passive scalar field. An animated version is also provided as supplementary material.

An attached boundary layer develops on the front face of the cylinder, starting from the stagnation point and progressing towards the sharp leading edges, where it separates, forming two laminar free-shear layers along the sides. Transition to turbulence occurs shortly after separation. The Kelvin–Helmholtz instability triggers the roll-up of the shear layers into a sequence of spanwise vortices. These vortices pair and undergo distortion under the influence of growing disturbances, eventually breaking down into small-scale turbulent motions. The evolution of the vortical structures is visualised in figure 1(a) using the $\lambda _2$ criterion by Jeong & Hussain (1995). The evolution of spanwise vortices into turbulence drives the rapid growth of the shear layer and exhibits similarities to a canonical mixing layer (Rogers & Moser 1992; Cimarelli et al. 2024). Reattaching shear layers, however, differ from plane mixing layers in several ways (Castro & Haque 1987). First, the low-speed side of the shear layer is highly turbulent, unlike the typically low turbulence levels in a plane mixing layer. Second, the reattaching shear layer exhibits a low-frequency flapping motion, which increases turbulence intensity by 5 %–10 % compared with plane mixing layers, albeit with minimal contribution to Reynolds shear stress (Eaton & Johnston 1982; Kiya & Sasaki 1983; Fang & Wang 2024). The instantaneous reattachment location fluctuates due to vortex shedding, leading to periodic reattachment–detachment cycles synchronised with large-scale vortex formation in the wake, see the discussion by Cimarelli et al. (2018b ).

Figure 2. Schematic of the flow configuration and definition of symbols. Shaded regions indicate areas with vorticity. The dashed line denotes the mean dividing streamline, $y_\psi (x)$ , the upper dotted line indicates the mean turbulent/non-turbulent interface, $y_\varOmega (x)$ , and the lower dotted line marks the boundary of the backflow region, $y_b(x)$ , on one side of the body. These lines are defined in the text. The arrows depict the conceptual flow model, with turbulent structures supplying the backflow near the wall.

The mean reattachment length of the turbulent shear layer, denoted as $x_r$ in figure 2, is governed by pressure recovery in the wake and the entrainment process of the shear layer (Smith, Pisetta & Viola 2021). The entrainment rate depends on the Reynolds number (Dimotakis & Brown 1976), which influences the onset of turbulence and the scale of turbulent structures. A higher entrainment rate causes a larger growth rate of the shear layer, thus reducing the mean reattachment length. The proximity of the shear layer to the wall further affects entrainment, limiting fluid transport across the lower interface (Cimarelli et al. 2024; Li & Wang 2024). However, the detailed mechanisms regulating momentum and energy transfer in the near-wall reattachment zone remain unsettled (Ota & Kon 1980; Gorin 2012).

The recirculating flow region, also known as the separation bubble, is bounded by the turbulent shear layer above and a near-wall backflow region below. The backflow region, characterised by relatively slow reverse flow, is sustained by the adverse pressure gradient, which deflects part of the shear-layer fluid upstream. As noted by Simpson, Chew & Shivaprasad (1981), mean backflow originates not only from downstream regions near reattachment but also intermittently from large-scale structures passing through the separated flow and supplying the near-wall flow, see again figure 2. Within the backflow region, a reverse boundary layer initially accelerates due to the pressure gradient, but it is prone to separate as it approaches the leading edge under the action of an adverse pressure gradient, forming a secondary separation bubble near the wall. The backflow contains weak but dynamically relevant motion, as evident in the blue-shaded regions of low streamwise velocity in figure 1(a). For instance, in backward-facing step flows, reverse flow velocities can exceed 20 % of the free stream velocity (Adams & Johnston 1988b ), with skin-friction coefficients as low as $C_f = -0.0033$ (Jović & Driver 1994). For the rectangular 5 : 1 cylinder, skin-friction values reach up to $C_f = -0.013$ (Corsini et al. 2022). This region lacks large-scale structures; instead, it consists of small-scale turbulent fluctuations generated by vortex impingements and driven upstream by the pressure gradient (Pronchick 1983; Cimarelli et al. 2018b ).

In the reattachment zone, the adverse pressure gradient and wall interaction significantly affect the shear layer. Reynolds normal and shear stresses decay rapidly in this region and continue to diminish downstream. This behaviour is observed for aspect ratios of the cylinder as high as 15 : 1 (Li et al. 2024b ). Concurrently, a boundary layer begins to form within the reattached shear layer. In its inner layer, mean velocity profiles follow a logarithmic law, similar to that of ordinary boundary layers. However, the outer part retains free-shear-layer characteristics for a considerable distance downstream of reattachment, with large eddies generated during separation persisting downstream. Transitioning to an ordinary boundary-layer structure is generally a slow process (Bradshaw & Wong 1972; Jović 1998), not accomplished in the 5 : 1 rectangle.

At the trailing edge, the turbulent boundary layer detaches, forming a wake characterised by free-shear layers and two counter-rotating recirculation bubbles in the base region. Vortex shedding occurs as large-scale oscillations due to alternating vortex formation on either side of the body. These oscillations dominate wake dynamics, contributing to aerodynamic drag, structural vibrations and scalar transport. Figure 1(b) shows how large-scale vortex shedding promotes downstream transport and mixing of passive scalar, with high-concentration scalar released from the body surface and recirculation region being entrained into the wake and wrapped into the alternating vortical structures. The shedding frequency, quantified by the Strouhal number $\textit {St}$ , depends on the Reynolds number ${\textit{Re}}$ and the aspect ratio of the cylinder (Nakamura et al. 1996; Chiarini et al. 2022b ; Li et al. 2024b ). For a 5 : 1 aspect ratio, $\textit {St}$ approaches a constant value of approximately 0.11 for ${\textit{Re}} \gt 1000$ (Schewe 2013; Moore et al. 2019). Cimarelli et al. (2024) attribute this asymptotic behaviour of $\textit {St}$ and other integral quantities, such as drag coefficients, to the flow approaching a high-Reynolds-number regime.

3.1. Phase-averaged statistics

Vortex shedding and the oscillatory motion of the separated flow region is a distinguishing, general feature of the flow. Decomposing the flow in a periodic and random part can aid in the characterisation of this oscillatory motion. Phase-averaged statistics are obtained based on the temporal evolution of the lift coefficient, which exhibits regular oscillations with a frequency corresponding to the vortex shedding. Around each lift extremum, a tolerance window of duration $\pm 1 D/U_0$ was defined. Snapshots near lift minima were mirrored about the horizontal symmetry plane $y = -0.5$ and averaged together with those near lift maxima to produce statistically coherent phase-averaged fields. In total, $35$ and $21$ snapshots were used for the cases at ${\textit{Re}}=3000$ and ${\textit{Re}}=14\,000$ , respectively. The two phases analysed, $\phi _{\textit{max}}$ and $\phi _{\textit{min}}$ , correspond to the averaged solution obtained on the upper and lower sides of the body, respectively. The case at ${\textit{Re}}=8000$ is omitted for brevity, as it presents intermediate behaviour.

Figure 3. Phase-averaged velocity streamlines and pressure field at (a) ${\textit{Re}}=3000$ and (b) ${\textit{Re}}=14\,000$ . The upper side of the rectangle corresponds to phase $\phi _{\textit{max}}$ , while the lower side corresponds to $\phi _{\textit{min}}$ . Cross symbols mark the mean reattachment points on the body sides.

Figure 4. Phase-averaged distributions of (a) skin friction coefficient, (b) Nusselt number, (c) pressure coefficient and (d) standard deviation of pressure coefficient. Line styles denote results at ${\textit{Re}} = 3000$ (dotted line) and ${\textit{Re}} = 14\,000$ (solid line). Curve colours indicate averages over all samples (black), and at two phases $\phi _{\textit{max}}$ (red) and $\phi _{\textit{min}}$ (blue).

The phase-averaged flow structure is displayed in figure 3, showing velocity streamlines superimposed to the pressure field for ${\textit{Re}} = 3000$ and $14\,000$ . In phase with the streamwise oscillation of the low-pressure region, the separation bubbles along the rectangle sides periodically elongate and contract, while its height remains nearly constant. A positive peak in the lift coefficient corresponds to a longer separation bubble on the upper side. The displacement of the reattachment point increases with ${\textit{Re}}$ . As displayed in figure 4(a), showing the phase-averaged friction coefficient $c_f$ , the reattachment coordinate $x_r$ is in the range $4.04 \lt x_r \lt 4.29$ at ${\textit{Re}}=3000$ and within $4.07 \lt x_r \lt 4.55$ at ${\textit{Re}}=14\,000$ . When the separation bubble shortens (at $\phi _{\textit{min}}$ ), the reverse flow strengthens, as apparent from the reduced skin friction in the interval $2\lt x\lt 4$ . However, separation of the reverse boundary layer is observed at both phases. This indicates that strong reverse flow does not penetrate into the upstream portion of the backflow.

The unsteadiness of the separation bubble affects scalar transfer at the wall, as shown by the local Nusselt number in figure 4(b). A longer separation bubble at $\phi _{\textit{max}}$ leads to enhanced scalar near reattachment and after it, while the reduced bubble at $\phi _{\textit{min}}$ corresponds to diminished transfer with respect to the mean. This difference is attributed to vortices impinging on the wall, which enhance the wall-normal scalar gradient more effectively than the reattached turbulent boundary layer. These impingements increase wall-normal velocity fluctuations and promote entrainment of fluid at low scalar concentration from the free stream. The phase difference in scalar transfer is amplified as ${\textit{Re}}$ increases.

Figure 4(c) shows the phase-averaged pressure coefficient $c_p$ . The pressure difference between the two phases generates the fluctuating lift force on the body and grows with ${\textit{Re}}$ along the entire lateral face, consistent with the increase in lift amplitude reported by Cimarelli et al. (2024) (see their figure 5d). While pressure minima maintain similar magnitudes across phases, their positions and pressure recovery rates change. At $\phi _{\textit{max}}$ (red curves), the pressure minimum shifts upstream with increasing ${\textit{Re}}$ , and the recovery becomes smoother, leading to lower pressure levels in the reattachment region. At $\phi _{\textit{min}}$ (blue curves), pressure is higher near the leading edge and recovers more steeply downstream, resulting in higher pressure in the reattached flow. The $c_p$ distribution motivates also the persistence of the secondary separation bubble, as the more energetic reverse flow is accompained by a stronger adverse pressure gradient.

Figure 4(d) shows the standard deviation of the phase-averaged pressure coefficient, $c_{p_{\textit{std}}}$ , which represents the amplitude of random pressure fluctuations at a fixed phase. According to phase-averaged decomposition (Cantwell & Coles 1983), total fluctuations from the global mean include contributions from periodic large-scale motion and from local random turbulence: $ p = P + \widetilde {P} + p^{\prime \prime }$ , where $p^\prime = \widetilde {P} + p^{\prime \prime }$ , $\widetilde {P}$ is the periodic mean component and $p^{\prime \prime }$ is the random component. Fluctuations at $\phi _{\textit{max}}$ and $\phi _{\textit{min}}$ are similar and decrease in intensity with increasing ${\textit{Re}}$ . The difference between the global mean (black curves) and phase-locked results (coloured curves) arises from the contribution of large-scale periodic fluctuations. The growing gap with ${\textit{Re}}$ suggests that periodic, coherent unsteadiness becomes increasingly dominant relative to random turbulence.

4. Momentum transport

4.1. Mean flow

The analysis begins with the mean momentum distribution around the 5 : 1 rectangular cylinder. Effects of a Reynolds number increase on the turbulent flow are considered. A variation in the Reynolds number changes the location of laminar to turbulent transition, affecting the pressure distribution and the momentum transfer along the shear layer and in the base region (Chapman, Kuehn & Larson 1958).

Figure 5. Mean (a,c) streamwise and (b,d) vertical velocity components at (a,b) ${\textit{Re}}=3000$ and (c,d) $14\,000$ . Dotted and dashed lines, from top to bottom, designate the turbulent/non-turbulent interface $y_\varOmega$ , the shear-layer centreline $y_{sl}$ and the boundary of the backflow region $y_b$ . Cross symbols indicate the locations of maximum and minimum values.

Figure 5 shows the contours of the streamwise ( $U$ ) and vertical ( $V$ ) mean velocities at ${\textit{Re}}=3000$ and $14\,000$ , where the crosses indicate velocity extrema. To provide a support for the interpretation of statistical quantities above the rectangle, figure 5, as well as most of those following, include indications of the turbulent/non-turbulent interface of the leading-edge shear layer $y_\varOmega (x)$ , the shear-layer centreline $y_{sl}(x)$ , the line encompassing the backflow region $y_b(x)$ and the dividing streamline $y_\psi (x)$ . The turbulent/non-turbulent interface $y_\varOmega$ can be interpreted as the upper boundary of the shear layer, as it lies between the free flow and the sheared region. It is defined here as the set of outermost points of the rotational region with absolute values of vorticity greater than $0.01|\varOmega _{z}|_{\textit{max}}(x)$ , where $|\varOmega _{z}|_{\textit{max}}(x)$ is the maximum absolute value of the mean spanwise vorticity at each $x$ -location. The shear-layer centreline $y_{sl}$ is identified by the position of $|\varOmega _{z}|_{\textit{max}}(x)$ . Concerning the boundary of the backflow region $y_b$ , it separates regions of the flow with different characteristics beneath the shear layer, as will be shown later, and is representative of the recirculation region. This boundary corresponds to the line where the mean streamwise velocity is zero, $U=0$ . Finally, the dividing streamline $y_\psi$ is the mean streamline originating from the leading edge, and thus connects the separation point of the shear layer to its mean reattachment point.

Two distinct regions of negative streamwise velocity are observed, one above the rectangle and the other in its wake. These regions correspond to two mean separation bubbles induced by the leading and trailing edges of the rectangle, hereafter referred to as primary and wake separation bubbles. Minima of the streamwise velocity component are found in the backflow region inside the primary separation bubble, reaching $U_{\textit{min}}=-0.34$ and $-0.37$ for the cases at ${\textit{Re}}=3000$ and $14\,000$ , respectively. In the case at high ${\textit{Re}}$ , low values of streamwise velocity extend more upstream toward the leading edge, due to the influence of a more extended favourable pressure gradient along the wall (see figure 7c of Cimarelli et al. 2024).

The shear layer separating from the leading edge bounds the regions of flow reversal ( $U \lt 0$ ) on the lower side and high streamwise mean velocity ( $U \gt U_0$ ) on the upper side. The maximum value of the streamwise velocity component is attained close to the irrotational flow region, where $U_{\textit{max}}=1.32$ is reached for both the Reynolds numbers considered. While the location of the leading-edge shear layer proves to be almost independent from ${\textit{Re}}$ in the range analysed, its spreading primarily depends on the transition process (and therefore on ${\textit{Re}}$ ) and on the distance from the wall. In figure 5, $y_{sl}$ remains constant with respect to ${\textit{Re}}$ , while $y_\varOmega$ begins to spread at varying distances from the leading edge. The proximity of the wall to one side of the free shear layer limits the entrainment of fluid on that side, inducing a lower pressure in that region. This pressure difference generates a downward curvature in the flow, which further constrains the inflow of fluid into the turbulent shear layer. The process intensifies until the shear layer impinges on the rectangle wall. The mean reattachment point, $x_r$ , nearly coincides with the intersection of $y_b$ with the rectangle wall. These locations are found to move slightly downstream for ${\textit{Re}}$ ranging between $3000$ and $14\,000$ (see figure 10 of Cimarelli et al. 2024), following a behaviour similar to the reattaching flow behind a backward-facing step (Adams & Johnston 1988a ).

The fluid entrained from the backflow region into the shear layer balances the flow reversed by the pressure rise through the reattachment zone. The flow rate per unit span entrained from the backflow region, $Q_e$ , increases from $0.0689 \, U_0 D$ at ${\textit{Re}} = 3000$ to $0.0754 \, U_0 D$ at ${\textit{Re}} = 14\,000$ . This quantity is computed as the line integral of the mean velocity field $\boldsymbol{U}$ across the curve $\varGamma$ , which represents the upstream boundary of the backflow region:

(4.1) \begin{equation} Q_e = \int _{\varGamma } \boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\hat {n}} \, {\textrm{d}}s. \end{equation}

In the equation, $\hat {n}$ is the unit vector normal to the curve and ${\textrm{d}}s$ is the arc length element. The curve $\varGamma$ is parametrised as $\varGamma = \lbrace (x, y_b(x)) \, | \, x \in [0, x_0]\rbrace$ , where $y_b(x)$ is the height at which the mean streamwise velocity vanishes and $x_0$ is the streamwise location where the volumetric flux changes sign (indicated by an asterisk in figure 6). Figure 6 shows the volume flux across $y_b$ as a function of the streamwise direction. Positive values of the volume flux indicate fluid flowing outward from the backflow region, while negative values correspond to incoming flow. The slight increase in $Q_e$ with ${\textit{Re}}$ is consistent with the enlargement of the primary separation bubble, as its vertical thickness remains approximately constant over the analysed ${\textit{Re}}$ range. As the Reynolds number increases and the transition of the laminar shear layer shifts upstream, a higher fraction of fluid is entrained from the backflow region into the shear layer close to the leading edge.

Figure 6. Volume flux across the line of zero streamwise velocity as a function of the streamwise distance normalised with the reattachment length: dotted line, ${\textit{Re}} = 3000$ ; dashed line, ${\textit{Re}} = 8000$ and solid line, ${\textit{Re}} = 14\,000$ . The asterisk symbol marks the zero-crossing point $x_0$ .

Maxima of the vertical mean velocity ( $V_{\textit{max}}=1.14$ and $1.20$ at ${\textit{Re}}=3000$ and $14\,000$ ) are observed at the leading edge, where the laminar boundary layer formed on the front face of the rectangle separates. The minimum value of the vertical velocity component is found above the reattachment region and corresponds to $V_{\textit{min}}=-0.172$ and $-0.156$ at ${\textit{Re}}=3000$ and $14\,000$ , respectively.

4.2. Reynolds stresses

Reynolds stresses reflect how turbulence redistributes momentum, affecting dominant phenomena such as shear-layer growth, reattachment and wake formation. In this section, we analyse their evolution across the Reynolds number range $3000 \lt {\textit{Re}} \lt 14\,000$ , also establishing connections to findings from prior works. Reynolds stresses in the flow around the rectangular 5 : 1 cylinder have been examined in previous numerical studies by Li et al. (2024a ), Cimarelli et al. (2018a ) and Chiarini & Quadrio (2021), based on DNS data ranging from ${\textit{Re}}=1000$ to $3000$ , as well as in the experimental study by Kumahor & Tachie (2022), based on PIV measurements at ${\textit{Re}}=16\,200$ .

Figure 7. Distributions of the Reynolds normal stresses (a,c) $\langle u^{\prime 2} \rangle$ and (b,d) $\langle v^{\prime 2} \rangle$ at (a, b) ${\textit{Re}}=3000$ and (c,d) $14\,000$ . The cross symbols mark the locations of the maximum values.

Figure 8. Distributions of the (a,c) Reynolds normal stress $\langle w^{\prime 2} \rangle$ and (b,d) shear stress $\langle u^{\prime } v^{\prime } \rangle$ at (a,b) ${\textit{Re}}=3000$ and (c,d) $14\,000$ . The cross symbols mark the locations of maximum value of $\langle w^{\prime 2} \rangle$ and the minimum value of $\langle u^{\prime } v^{\prime } \rangle$ .

In the statistically two-dimensional flow considered, the only non-zero components of the Reynolds-stress tensor are the normal stresses $\langle u^{\prime 2} \rangle$ , $\langle v^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ , and the shear stress $\langle u^{\prime } v^{\prime } \rangle$ . Throughout this section, we use the term Reynolds stress to refer to quantities $\langle u_i^{\prime } u^{\prime}_{\!j} \rangle$ , excluding the negative sign that is introduced when these terms are transferred from the left-hand side to the right-hand side of the mean momentum equation to represent turbulent stresses. The spatial evolution of the four terms in the present flow configuration is shown in figures 7 and 8 for the cases at ${\textit{Re}} = 3000$ and $14\,000$ . For sharp-edged bluff bodies with uniform incoming flow, the regions of maximum Reynolds stresses are strictly related to the process of laminar to turbulent transition and occur along the turbulent free shear layers. Turbulent fluctuations are initiated by the Kelvin–Helmholtz-like roll-up of the vortex sheet separating from the leading edge and by the pairing of vortices. By increasing the Reynolds number, the instability mechanisms are triggered at shorter distances from the leading edge and turbulence arises more rapidly under the action of mean shears of growing intensity. To illustrate the deformation experienced by the fluid in the shear layer, the principal mean strain rate $S_{\lambda }$ , i.e. the largest eigenvalue of the mean strain-rate tensor $ S_{\textit{ij}} = ({1}/{2}) ( \partial U_i / \partial x_j + \partial U_j / \partial x_i )$ , is evaluated along the shear-layer centreline $y_{sl}$ and shown in figure 9. Note that $S_{\textit{ij}}$ is a second-order symmetric and traceless tensor, thus its eigenvalues are real and opposite, indicating stretching ( $S_{\lambda }$ ) in one direction and compression ( $-S_{\lambda }$ ) in the other. The data in figure 9 are presented on a logarithmic scale to highlight the power-law behaviour of $S_{\lambda }$ . In the initial, laminar stage of the shear layer, $S_{\lambda }$ exhibits larger values with ${\textit{Re}}$ , but converges to similar values at the trailing edge over the range of ${\textit{Re}}$ examined. Two distinct power-law decay regions, marked by the dash-dotted lines in figure 9, are observed. The first region exhibits a weaker decay, associated with transitional turbulence mechanisms, while the second shows a steeper decay, corresponding to the fully turbulent shear layer. This variation indicates a change in the deformation of fluid elements due to turbulence. The decay rates are found to be almost invariant with respect to ${\textit{Re}}$ and can be associated with the different growth stages of the shear layer thickness, as described by Cimarelli et al. (2024).

Figure 9. Largest eigenvalue of the mean strain-rate tensor along the leading-edge shear layer centreline: dotted line, ${\textit{Re}} = 3000$ ; dashed line, ${\textit{Re}} = 8000$ ; and solid line, ${\textit{Re}} = 14\,000$ . The dash-dotted lines mark the two power-law decay regions.

The streamwise Reynolds normal stress $\langle u^{\prime 2} \rangle$ is the dominant component in the leading-edge shear layer for all the Reynolds numbers considered. At ${\textit{Re}}=3000$ , high streamwise fluctuations arise in the rear part of the primary separation bubble, reaching a peak of $0.116$ at $(x, y) = (3.37, 0.45)$ . As shown by Cimarelli et al. (2024), in the case at ${\textit{Re}} = 3000$ , the intense turbulent fluctuations originating in this region are due to the positive combination of the transition process and the shedding of large-scale vortices from the separated region, occurring in close proximity. As the Reynolds number increases, transition mechanisms shift upstream along the shear layer centreline and separate (at least in terms of distances) from the vortex shedding region. Hence, levels of $\langle u^{\prime 2} \rangle$ remain significant across the bulk of the shear layer and aligned with the streamwise direction. The maximum values of $\langle u^{\prime 2} \rangle$ are reached in the early turbulent region and increase with ${\textit{Re}}$ . In detail, the peak values of $0.114$ and $0.122$ are found at $x=1.38$ and $0.79$ for ${\textit{Re}} = 8000$ and $14\,000$ , respectively.

The vertical and spanwise Reynolds normal stresses $\langle v^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ are characterised by a lower intensity in the leading-edge shear layer with respect to the streamwise component. This difference becomes wider with the growth of the Reynolds number. In particular, while at ${\textit{Re}}=3000$ , the maximum values of $\langle v^{\prime 2} \rangle _{\textit{max}} = 0.090$ and $\langle w^{\prime 2} \rangle _{\textit{max}} = 0.090$ are observed respectively at $(2.96,0.47)$ and $(3.23,0.40)$ , and at ${\textit{Re}}=14\,000$ , $\langle v^{\prime 2} \rangle _{\textit{max}} = 0.067$ and $\langle w^{\prime 2} \rangle _{\textit{max}} = 0.080$ are found at $(0.71, 0.38)$ and $(0.89, 0.42)$ , respectively. Vertical fluctuations on the low-velocity side of the shear layer rapidly fall as the solid wall is approached. Due to the curvature of the mean streamline close to the reattachment region, turbulent motions are conveyed towards the wall and wall-normal fluctuations are suppressed because of the impermeability condition. On the other side, fluctuations parallel to the wall are generated close to the wall in the reattachment region and in the backflow region at decreasing distances from the leading edge with ${\textit{Re}}$ .

The streamwise Reynolds normal stress also prevails in the shear layer separating from the trailing edge of the rectangle. Here, turbulent motions developed in the reattached boundary layer, while transported downstream, are amplified by the mean velocity gradients in the wake region and acquire streamwise turbulent energy. The magnitude of $\langle u^{\prime 2} \rangle$ in this region increases with the Reynolds number. In the case at ${\textit{Re}} = 14\,000$ , ${\langle u^{\prime 2} \rangle }_{\textit{max}} = 0.084$ is obtained at $(6.03, 0.07)$ . Unlike the leading-edge shear layer, high streamwise fluctuations are not centred in the bulk of the trailing-edge shear layer, but are spread across the vertical direction. One of the reasons lies in the vertical oscillatory motion induced on the trailing-edge shear layer by the von Kármán vortex shedding motion occurring in the wake region. Other factors that contribute to the fast growth of the shear-layer thickness are the turbulent state of the boundary layer separating at the trailing edge and the supply of turbulent motions generated upstream in the reattached boundary layer.

The large scale oscillations in the vertical direction occurring in the wake are associated with high vertical fluctuations centred in the wake centreline, i.e. the symmetry plane $y=-0.5$ . More intense oscillations take place at greater Reynolds numbers, as shown by the significant rise of $\langle v^{\prime 2} \rangle$ for $x\gt 5$ in figure 7(b,d). The maximum value in the wake increase from $\langle v^{\prime 2} \rangle _{\textit{max}} = 0.070$ at ${\textit{Re}}=3000$ to $0.074$ and $0.105$ at ${\textit{Re}}=8000$ and $14\,000$ , respectively. Moreover, high levels of $\langle v^{\prime 2} \rangle$ last over longer distances in the wake at larger ${\textit{Re}}$ (not shown in figure).

The Reynolds shear stress $\langle u^{\prime } v^{\prime } \rangle$ provides a relevant contribution to momentum transfer across the flow, as well as to the production of turbulent kinetic energy. As shown in figure 8(b, d), $\langle u^{\prime } v^{\prime } \rangle$ is mostly negative around the rectangle and exhibits high magnitude levels in the core of the turbulent leading-edge shear layer and in the wake. The minimum value, located where streamwise and vertical fluctuations are best anti-correlated, reaches $-0.053$ and is located within the separation bubble at $(3.21, 0.47)$ for ${\textit{Re}}=3000$ , while equals $-0.045$ and is located in the wake region at $(6.89, 0.064)$ for ${\textit{Re}}=14\,000$ . Regions of positive values for the shear stress $\langle u^{\prime } v^{\prime } \rangle$ are limited to the most upstream part of the separation bubble below the shear layer and within the wake separation bubble.

To enable a quantitative comparison of the Reynolds stresses generated in each section $x$ at different Reynolds numbers, we compute the integrals of the Reynolds stress tensor components, $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ , across the turbulent flow, from the rectangle wall to the outer edge of the shear layer $y_{\varOmega }$ ,

(4.2) \begin{equation} \gamma _{\textit{ij}}(x) = \int _{0}^{y_{\varOmega }} \big|\big\langle u_{i}^{\prime } u_{j}^{\prime } \big\rangle (x,y)\big| \, {\textrm{d}} y, \end{equation}

where the absolute value ensures that $\gamma _{12}$ is positive. The resulting integrals for the three ${\textit{Re}}$ values investigated are reported in figure 10. Globally, a larger amount of $\langle uu\rangle$ rises along the shear layer at higher Reynolds numbers. However, cross-stream Reynolds normal stresses $\langle v^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ do not experience a similar increase. Although higher vertical and spanwise fluctuations occur near the leading edge by increasing the Reynolds number, lower values of $\gamma _{22}$ and $\gamma _{33}$ are observed in the downstream half of the rectangle. The same is true for the shear component $\gamma _{12}$ .

Figure 10. Integrals of the Reynolds stresses across the separated flow above the rectangle, (a) $\gamma _{11}$ , (b) $\gamma _{22}$ , (c) $\gamma _{33}$ and (d) $\gamma _{12}$ : dotted line, ${\textit{Re}} = 3000$ ; dashed line, ${\textit{Re}} = 8000$ ; and solid line, ${\textit{Re}} = 14\,000$ .

4.3. Mean momentum balance

Turbulence significantly affects the mean flow (see § 3) by redistributing momentum through Reynolds stresses. This section quantifies this effect by analysing the mean momentum balance. We highlight how the pressure gradient and the divergence of the Reynolds stress tensor act to bend the shear layer, promote reattachment and sustain the backflow. To isolate these effects, we reformulate the momentum equations in a streamline-aligned coordinate system, which allows direct interpretation of turbulent momentum transport across the mean dividing streamline.

Previous sections compared flow features across three different Reynolds numbers. From this point forward, the analysis focuses on the case at ${\textit{Re}} = 14\,000$ , where the flow approaches an asymptotic high-Reynolds-number regime.

For the statistically two-dimensional and steady flow under consideration, the conservation equations for mean momentum in the streamwise and vertical directions reduce to

(4.3) \begin{equation} \displaystyle \displaystyle U_j \frac {\partial U_i}{\partial x_j} = \displaystyle - \frac {\partial P}{\partial x_i} - \frac {\partial \big\langle u_{i}^{\prime } u_{j}^{\prime } \big\rangle }{\partial x_j}. \end{equation}

The left-hand side represents the mean advection of momentum that is balanced by the gain or loss in mean momentum resulting from the action of the mean pressure field and the Reynolds stresses. The viscous diffusion term is omitted from the mean momentum balance, as it is negligible across most of the domain for the Reynolds numbers considered, except in the initial part of the separated shear layer, where turbulence is still developing and Reynolds stresses remain low. This simplification was verified by comparing its magnitude to those of the pressure gradient and Reynolds stress divergence.

Figure 11. Field lines of (a) the mean pressure gradient $(-\partial P / \partial x, -\partial P / \partial y)$ and (b) the divergence of the Reynolds stress tensor $(-\partial \langle u^{\prime } u_j^{\prime }\rangle / \partial x_j, -\partial \langle v^{\prime } u_j^{\prime }\rangle / \partial x_j)$ , with $j=1,2$ . The colour represents the magnitude of the vector field. Results are evaluated at ${\textit{Re}} = 14\,000$ . The green line marks the dividing streamline originating from the leading edge $y_\psi$ .

The effect of the mean pressure gradient force on the flow is first examined. Figure 11(a) shows the vector field of the mean pressure gradient around the rectangle. The incoming flow about $y = -0.5$ loses streamwise momentum because of the increasing pressure levels near the stagnation point, $(x,y) = (0, -0.5)$ , and deviates vertically, gaining $y$ -momentum, under the influence of a favourable pressure gradient. Downstream the leading edge, curved mean streamlines which circumvent the front corner are accompanied by a transverse mean pressure gradient that is proportional to the centripetal acceleration. Peak pressure gradients are thus observed near the leading edge, $(x,y) = (0, 0)$ , in the region where high velocities occur together with small radii of curvature of the streamlines. In the shear layer, the pressure gradient field grows in intensity (indicated by a darker shade), starting from the initial stages of the turbulent shear layer. The pressure gradient peak is located at $x \approx 0.4$ at ${\textit{Re}} = 14\,000$ . Here, the pressure gradient force is directed towards the shear-layer centreline, drawing fluid from the outside region. Downstream the centre of the primary separation bubble, $(x,y) = (2.17, 0.33)$ , the mean pressure raises and slows down the streamwise component of the mean flow. The presence of the impermeable wall is felt by the mean flow moving downstream through a vertical adverse pressure gradient. This downwash is observed starting from $x \approx 2.2$ .

The divergence of the Reynolds stress tensor, shown in figure 11(b), represents the momentum changes due to turbulent transport within the flow. Turbulent motion significantly affects the momentum distribution along the separated shear layers and in the vicinity of the reattachment region. In the leading-edge shear layer, two peaks of Reynolds stress divergence are observed on the upper and lower sides. On the upper side, the divergence of Reynolds stresses removes streamwise ( $x$ ) momentum while supplying wall-normal ( $y$ ) momentum. Conversely, on the lower side, $x$ -momentum is gained while $y$ -momentum is drained. Due to the transport of mean momentum by turbulent motion, fluid is pushed from the shear layer centreline to the free stream and the backflow region. Turbulence is responsible for the spreading of the shear layer and for smoothing the longitudinal velocity gradient between its two sides, similar to what is observed in plane mixing layers (Castro & Haque 1987; Rogers & Moser 1994). The $x$ -momentum lost in the upper side is released in the lower side. As the shear layer spreads and approaches the wall, turbulent transport extracts $y$ -momentum from the region near the wall and promotes flow reattachment.

To clarify the mechanisms regulating shear layer bending and reattachment, we analyse the dominant terms in the mean momentum balance along the dividing streamline, $y_\psi$ , which originates from the leading edge (see the green solid line in figure 11). For this purpose, we examine the mean momentum balance in the streamline coordinate system $(\tau , \eta , z)$ (Finnigan 1983), where $\tau$ is locally aligned with the mean velocity vector, $z$ corresponds to the Cartesian $z$ -axis, and $\eta$ is orthogonal to both $\tau$ and $z$ by the right-hand rule.

Reformulating the mean momentum equations in the streamline coordinate system offers several advantages. First, the terms in the transformed equations distinctly describe transport along or normal to the streamlines, removing any ambiguity caused by misaligned coordinates relative to the mean flow. In the second instance, it facilitates the direct interpretation of the Reynolds shear stress as the turbulent momentum transport across streamlines. Furthermore, the transformation emphasises quantities related to mean flow acceleration and streamline curvature, which are relevant characteristics of distorted shear flows. Finnigan et al. (1990) were the first to apply the streamline flow equations developed by Finnigan (1983) to investigate flow over a hill. Recent studies employing streamline coordinate analysis to capture the properties of external flows include those by Morse & Mahesh (2021), Plasseraud, Kumar & Mahesh (2023) and Prakash et al. (2024).

The mean momentum equations are expressed in the streamline coordinate system as follows:

(4.4) \begin{equation} \left . \begin{array}{rr} \displaystyle U_{\tau } \frac {\partial U_{\tau }}{\partial \tau } &= \displaystyle - \frac {\partial P}{\partial \tau } - \frac {\partial \big\langle u_{\tau }^{\prime }u_{\tau }^{\prime } \big\rangle }{\partial \tau } - \frac {\partial \big\langle u_{\tau }^{\prime }u_{\eta }^{\prime } \big\rangle }{\partial \eta } + \frac {\big\langle u_{\tau }u_{\tau } \big\rangle - \big\langle u_{\eta }^{\prime }u_{\eta }^{\prime } \big\rangle }{L_{a}} + 2 \frac {\big\langle u_{\tau }^{\prime }u_{\eta }^{\prime } \big\rangle }{R},\\[10pt] \displaystyle \frac {U^{2}_{\tau }}{R} &= \displaystyle - \frac {\partial P}{\partial \eta } - \frac {\partial \big\langle u_{\tau }^{\prime }u_{\eta }^{\prime } \big\rangle }{\partial \tau } - \frac {\partial \big\langle u_{\eta }^{\prime }u_{\eta }^{\prime } \big\rangle }{\partial \eta } + 2 \frac {\big\langle u_{\tau }^{\prime }u_{\eta }^{\prime } \big\rangle }{L_{a}} + \frac {\big\langle u_{\tau }^{\prime }u_{\tau }^{\prime } \big\rangle - \big\langle u_{\eta }^{\prime }u_{\eta }^{\prime } \big\rangle }{R}, \end{array}\right \} \end{equation}

where viscous terms have been again neglected. In (4.4), the subscripts denote the directional components of the velocity vector (e.g. $U_{\tau }$ is the mean velocity component in the tangential $\tau$ direction). Here, $L_{a}$ and $R$ are the local radii of curvature of the $\eta$ and $\tau$ coordinate lines, respectively, and are defined by the following relations:

(4.5) \begin{align} \frac {1}{L_{a}} &= \frac {1}{U_{\tau }} \frac {\partial U_{\tau }}{\partial \tau }, \end{align}

(4.6) \begin{align} \frac {1}{R} &= \frac {1}{U_{\tau }}\left ( \varOmega + \frac {\partial U_{\tau }}{\partial \eta } \right )\!. \end{align}

In the streamline and streamline-normal momentum equations (4.4), the left-hand side represents advection along a streamline and the mean centripetal acceleration, respectively. On the right-hand side, additional terms appear in the expression for Reynolds stress divergence due to the curvilinear nature of the coordinate system adopted. For a detailed derivation of the streamline-coordinate momentum equations, the reader is referred to Finnigan (1983).

Figure 12. Mean momentum balance evaluated along the dividing streamline $y_\psi$ . The dotted vertical line indicates the location of mean reattachment $x_r$ . (a) Terms in the streamline momentum equation. (b) Terms in the streamline-normal momentum equation. (c) Vectors of the mean pressure gradient (blue arrows) and the divergence of the Reynolds stress tensor (red arrows) along the dividing streamline $y_\psi$ , outlined by the dashed line.

Figure 12 shows the terms from the streamline and streamline-normal momentum equations along the dividing streamline $y_\psi$ . The mean advection term is substantially driven by the pressure gradient in both the tangential ( $\tau$ ) and normal ( $\eta$ ) directions. Turbulence effects are confined to the early turbulent and reattachment regions, where stress divergence terms gain relevance due to fluctuations production process and wall impermeability. The following analysis focuses on the range $0 \lt x \lt x_r$ .

In the streamline balance (see figure 12 a) the tangential pressure gradient supplies $\tau$ -momentum in the early turbulent region ( $x\lt 1$ ), compensating the negative peak of transport due to Reynolds tangential stress, $-\partial \langle u_{\tau }^{\prime }u_{\tau }^{\prime } \rangle / \partial \tau$ , and inducing a net acceleration along the streamline. At greater streamwise distances, $\tau$ -momentum is removed ( $U_{\tau } \partial U_{\tau } / \partial \tau \lt 0$ ) under the influence of the adverse tangential pressure gradient. When the flow approaches the wall ( $x\gt 3$ ), turbulent transport terms deliver momentum, balancing the pressure gradient at the mean reattachment location. Here, the Reynolds stress gradient $-\partial \langle u_{\eta }^{\prime }u_{\eta }^{\prime } \rangle / \partial \eta$ is the leading term, peaking at $x \simeq 4$ . Explicit effects of the curvature are represented by $(\langle u_{\tau }^{\prime }u_{\tau }^{\prime } \rangle - \langle u_{\eta }^{\prime }u_{\eta }^{\prime } \rangle ) / L_{a}$ , which extracts tangential momentum in the reattachment region.

As the trajectory of the dividing streamline indicates, fluid at the leading edge initially experiences a negative centripetal acceleration, i.e. the local centre of curvature lies in the direction of decreasing $\eta$ . In examining the streamline-normal balance (figure 12 b), the contributions to the centripetal acceleration are dominated by two main terms, the normal pressure gradient and the transport due to normal Reynolds stress, $-\partial \langle u_{\eta }^{\prime }u_{\eta }^{\prime } \rangle / \partial \eta$ . The magnitudes of both terms peak in the early turbulent region, with turbulent transport providing a centrifugal contribution, $-\partial \langle u_{\eta }^{\prime }u_{\eta }^{\prime } \rangle / \partial \eta \gt 0$ , and the pressure gradient contributing centripetally, $-\partial P / \partial \eta \lt 0$ . As the turbulent shear layer evolves ( $x\gt 1$ ), the turbulent transport term becomes negligible and the centripetal acceleration along the curved streamline is dictated solely by the normal pressure gradient. The pressure gradient term decreases in magnitude with the streamwise distance and changes its sign to positive while approaching the wall, inducing a change in streamline curvature ( $U^{2}_{\tau } / R \gt 2$ ) at $x \simeq 4$ . From reattachment onwards, the dividing streamline flattens against the wall, exhibiting zero curvature $1/R \approx 0$ . Here, the normal pressure gradient balances the normal Reynolds stress gradient.

Figure 12(c) provides an overview of the acceleration terms from pressure and turbulent stresses in the mean momentum balance along the dividing streamline. Notice that these two contributions are invariant in Cartesian and streamline coordinate systems, provided all terms emerging in the expressions are combined into a single vector. In general, the Reynolds stress divergence counteracts the pressure gradient, working to increase the thickness of the separation bubble in the early turbulent region of the shear layer and promoting the flow reattachment in the downstream region of the rectangle.

4.4. Reynolds-stress balances

This section investigates the production, redistribution, transport and dissipation of the individual components of the Reynolds-stress tensor. While the previous study by Chiarini & Quadrio (2021) has provided a detailed description of the dominant mechanisms at ${\textit{Re}}=3000$ , particularly the role of production in the separated shear layer and intercomponent redistribution via the pressure–strain tensor, here we extend the analysis to the highest Reynolds number case. Attention is given to the scaling and spatial displacement of peak values, as well as changes in the relative weight between production, dissipation and pressure–strain terms as the Reynolds number increases up to ${\textit{Re}} = 14\,000$ .

The budgets for the Reynolds stresses can be written as

(4.7) \begin{equation} \frac {\partial T_{k,{\textit{ij}}}}{\partial x_{k}} = P_{\textit{ij}} - \epsilon _{\textit{ij}} + \varPi _{\textit{ij}}, \end{equation}

where the left-hand side term represents the mean, pressure, turbulent and diffusive transport terms, all of which are combined and expressed as the divergence of the Reynolds-stress flux $T_{k,ij}$ . The Reynolds-stress flux is defined as

(4.8) \begin{equation} T_{k,ij} = U_{k} \big\langle u_{i}^{\prime } u_{j}^{\prime } \big\rangle + \big [\big\langle p^{\prime } u_{i}^{\prime }\big\rangle \delta _{jk} + \big \langle p^{\prime } u_{j}^{\prime }\big\rangle \delta _{ik} \big ] + \big\langle u_{k}^{\prime }u_{i}^{\prime }u_{j}^{\prime }\big\rangle - \frac {1}{{\textit{Re}}} \frac {\partial \big\langle u_{i}^{\prime } u_{j}^{\prime } \big\rangle }{\partial x_{k}}, \end{equation}

where $\delta _{\textit{ij}}$ is Kronecker’s delta. The terms on the right-hand side of (4.7) correspond to the production tensor $P_{\textit{ij}}$ , the dissipation tensor $\epsilon _{\textit{ij}}$ and the pressure–strain tensor $\varPi _{\textit{ij}}$ . These are defined as follows:

(4.9) \begin{align} P_{\textit{ij}} &= -\left\langle u_i^{\prime } u_k^{\prime } \right\rangle \frac {\partial U_{j}}{\partial x_{k}} -\big\langle u_j^{\prime } u_k^{\prime } \big\rangle \frac {\partial U_{i}}{\partial x_{k}}, \end{align}

(4.10) \begin{align} \epsilon _{\textit{ij}} &= \frac {2}{{\textit{Re}}} \left\langle \frac {\partial u_i^{\prime }}{\partial x_k} \frac {\partial u_j^{\prime }}{\partial x_k} \right\rangle\!. \end{align}

(4.11) \begin{align} \varPi _{\textit{ij}} &= \left \langle p^{\prime }\, \frac {\partial u_i^{\prime }}{\partial x_j} \right \rangle + \left \langle p^{\prime }\, \frac {\partial u_j^{\prime }}{\partial x_i} \right \rangle\! . \end{align}

These terms can be combined to form the net source term $\varPsi _{\textit{ij}} = P_{\textit{ij}} - \epsilon _{\textit{ij}} + \varPi _{\textit{ij}}$ , representing the overall source or sink of Reynolds stresses. In the region where $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle \gt 0$ , a negative value of $\varPsi _{\textit{ij}}$ ( $\varPsi _{\textit{ij}} \lt 0$ ) indicates a loss, whereas a positive value ( $\varPsi _{\textit{ij}} \gt 0$ ) indicates a gain. For cases where $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle \lt 0$ , such as for $\langle u^{\prime } v^{\prime } \rangle$ , the opposite is true, i.e. $\varPsi _{\textit{ij}} \lt 0$ denotes a source of the shear stress, corresponding to an increase in the magnitude of the shear stress, and $\varPsi _{\textit{ij}} \gt 0$ indicates a sink. Each term in (4.7) is indexed according to the specific Reynolds stress component under consideration.

The non-zero components of the production tensor $P_{\textit{ij}}$ are shown in figure 13. Note that $P_{33}$ is zero, as the mean spanwise velocity and its spatial derivatives are zero. The production terms attain their highest values within the leading and trailing shear layers, where the work of the Reynolds stress against the mean velocity gradient is enhanced by intense turbulent fluctuations and high mean strain rates. In the leading-edge shear layer, $P_{11}$ acts as a source for $\langle u^{\prime 2} \rangle$ , drawing energy from the mean flow, while $P_{22}$ serves as a sink for $\langle v^{\prime 2} \rangle$ , transferring energy from the fluctuating motion to sustain the mean vertical flow. This reversed energy exchange, from turbulence to mean flow, intensifies with increasing Reynolds number. At ${\textit{Re}} = 3000$ , the negative peak of $P_{22}$ is approximately $9\,\%$ of the maximum of $P_{11}$ . By ${\textit{Re}} = 14\,000$ , the negative production level reaches $34\,\%$ . The increased energy transfer from turbulence to the mean vertical flow likely contributes to the weak Reynolds number sensitivity observed of the thickness of the main separation bubble. As shown by Cimarelli et al. (2024), the separated region preserves its vertical size across the analysed ${\textit{Re}}$ range. The large-scale oscillations in the lee of the body correspond to a region where the production rate $P_{22}$ is positive, driven by the vertical deceleration of the mean flow ( $\partial V / \partial y \lt 0$ ) and significant vertical fluctuations $\langle v^{\prime 2} \rangle$ . The shear-stress production term $P_{12}$ is generally negative around the rectangle, except for small positive values in the backflow region, where the adverse pressure gradient induces the separation of the backward boundary layer. The negative values of $P_{12}$ , combined with $\langle u^{\prime } v^{\prime } \rangle \lt 0$ (see figure 8 d), indicate areas of positive production of $\langle u^{\prime } v^{\prime } \rangle$ , and thus enhance the shear stress magnitude. Table 1 shows the location and magnitude of peaks of the main components of the production tensor. As the Reynolds number increases, all components exhibit a marked intensification, with $P_{11}$ growing from approximately 0.5 at ${\textit{Re}} = 3000$ to nearly 2.5 at ${\textit{Re}} = 14\,000$ , and the peaks shifting upstream and closer to the wall. The component $P_{22}$ becomes significantly more negative, while $P_{12}$ also increases in magnitude, indicating intensified shear production in the separated region.

Table 1. Peak values of production tensor terms at increasing Reynolds numbers.

Figure 13. Production tensor terms (a) $P_{11}$ , (b) $P_{22}$ , and (c) $P_{12}$ at ${\textit{Re}}=14\,000$ . The cross symbols mark the locations of the maximum or minimum values of $P_{\textit{ij}}$ .

Figure 14 displays the non-zero components of the dissipation tensor, $\epsilon _{\textit{ij}}$ , highlighting regions of significant viscous dissipation within the flow. The dissipation $\epsilon _{12}$ is not shown as negligible. High dissipation levels are concentrated within the leading-edge shear layer, reaching peaks in the same locations of turbulence production. As the turbulent shear layer develops downstream and approaches the wall near reattachment, dissipation rates decrease. In comparison to the production and pressure-rate-of-strain tensors, the dissipation terms exhibit relatively lower magnitudes. The peak values of dissipation components, summarised in table 2, increase with the Reynolds number, albeit less steeply than the production terms. The dissipation peaks for $\epsilon _{11}$ and $\epsilon _{33}$ nearly triple from ${\textit{Re}} = 3000$ to ${\textit{Re}} = 14\,000$ , and their locations move upstream and downward, following the earlier onset of transition in the shear layer.

Table 2. Peak values of dissipation tensor terms at increasing Reynolds numbers.

Figure 14. Dissipation tensor terms (a) $\epsilon _{11}$ , (b) $\epsilon _{22}$ and (c) $\epsilon _{33}$ at ${\textit{Re}}=14\,000$ . The cross symbols mark the locations of the maximum values of $\epsilon _{\textit{ij}}$ .

Figure 15. Pressure–strain tensor terms (a) $\varPi _{11}$ , (b) $\varPi _{22}$ , (c) $\varPi _{33}$ and (d) $\varPi _{12}$ at ${\textit{Re}}=14\,000$ . The cross symbols mark the locations of the maximum or minimum values of $\varPi _{\textit{ij}}$ .

Figure 16. Pressure–strain normal terms $\varPi _{ii}$ evaluated along (a) the leading-edge shear layer centreline and (b) in the near wall region, at $y=0.024$ : solid line, $\varPi _{11}$ ; dashed line, $\varPi _{22}$ ; and dotted line, $\varPi _{33}$ .

In the Reynolds-stress equations (4.7), the pressure–strain tensor accounts for the intercomponent energy transfer and the tendency towards isotropy of the fluctuating motion. Figure 15 shows the pressure–strain tensor terms $\varPi _{\textit{ij}}$ . In the separated shear layer, the pressure–strain correlation redistributes energy from the streamwise component $\langle u^{\prime 2} \rangle$ ( $\varPi _{11}\lt 0$ ) into the cross-flow component $\langle v^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ ( $\varPi _{22}\gt 0$ and $\varPi _{33}\gt 0$ ). This redistribution reaches a maximum in the early turbulent region of the shear layer. As displayed in figure 16(a), showing the normal terms $\varPi _{ii}$ along the centreline of the shear layer $y_{sl}$ , the energy transferred from $\langle u^{\prime 2} \rangle$ to $\langle v^{\prime 2} \rangle$ is larger than that transferred to $\langle w^{\prime 2} \rangle$ . Along $y_{sl}$ , the peak values $\varPi _{22,\,{max}}$ and $\varPi _{33,\,{max}}$ amount to $82\,\%$ and $26\,\%$ of $|\varPi _{11}|_{\textit{max}}$ , respectively. In the energy budget for $\langle w^{\prime 2} \rangle$ , there is no production term, making $\varPi _{33}$ the sole source of energy for the spanwise fluctuations. Thus, energy from the mean flow is primarily transferred to $\langle u^{\prime 2} \rangle$ , while $\langle v^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ gain energy from the streamwise fluctuations through the action of pressure fluctuations. Near the wall, the impermeability condition causes the pressure–strain normal terms to dampen the vertical velocity fluctuations and transfer energy from $\langle v^{\prime 2} \rangle$ to the wall-parallel components, $\langle u^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ . This process is associated with the strong vertical deceleration of vortical structures in the turbulent shear layer as they impinge upon the wall, driven by the reattaching mean flow. Figure 16(b) shows the profiles of $\varPi_{\textit{ii}}$ close to the wall at $y=0.024$ , the distance at which $\overline {y^+} = 10$ . Here, the overline denotes the streamwise average and the superscript $+$ indicates the normalisation with friction units. In this region, the pressure–strain correlations are significant in the reattachment area, though their magnitude is lower than in the separated shear layer. As the Reynolds number increases, the rate of energy drained from $\langle v^{\prime 2} \rangle$ near the wall and redistributed decreases. At ${\textit{Re}}=3000$ (not shown in the figure), a larger portion of energy is transferred to $\langle w^{\prime 2} \rangle$ , but this imbalance decreases with increasing Reynolds number, and at ${\textit{Re}} = 14\,000$ , the energy exchange between $\langle u^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ becomes nearly equal. For the off-diagonal term $\varPi _{12}$ , it is commonly understood that pressure fluctuations introduce a term in the $\langle u^{\prime } v^{\prime } \rangle$ budget proportional to $-\langle u^{\prime } v^{\prime } \rangle$ , reducing the absolute value of the Reynolds shear stress (Monin & Yaglom 1975). This reflects the natural tendency of turbulence to decay towards isotropy. Figure 15(d) shows $\varPi _{12}$ as positive within the core of the leading-edge shear layer, but exhibiting negative values on the upper side for $2\lt x\lt 4$ . Changes in pressure–strain peaks, reported in table 3, indicate a strengthened intercomponent energy redistribution with increasing Reynolds number. In particular, $\varPi _{11}$ exhibits increasingly intense negative values, while $\varPi _{22}$ and $\varPi _{33}$ show stronger positive peaks. For ${\textit{Re}} \lt 14\,000$ , approximately $75\,\%$ of the energy removed from the streamwise component is redirected to the vertical component, suggesting that the relative partitioning of turbulent kinetic energy within the evolving shear layer is preserved as ${\textit{Re}}$ increases.

Table 3. Peak values of pressure-strain tensor terms at increasing Reynolds numbers.

Figure 17. Distribution of the net source and sink terms, $\varPsi _{\textit{ij}} = P_{\textit{ij}} - \epsilon _{\textit{ij}} + \varPi _{\textit{ij}}$ overlaid with flux lines of $\boldsymbol{T_{\textit{ij}}}$ : (a) $\varPsi _{11}$ and $\boldsymbol{T_{11}}$ ; (b) $\varPsi _{22}$ and $\boldsymbol{T_{22}}$ ; (c) $\varPsi _{33}$ and $\boldsymbol{T_{33}}$ ; and (d) $\varPsi _{12}$ and $\boldsymbol{T_{12}}$ .

The combined sink and source terms for each of the Reynolds stress budgets (4.7),  $\varPsi _{\textit{ij}}$ , are analysed along with the lines of the Reynolds-stress flux vector $\boldsymbol{T_{\textit{ij}}}$ in figure 17. This representation summarises locations of net source, sink and paths of Reynolds stress transfer. In the $\langle u^{\prime 2} \rangle$ budget, energy is provided to the streamwise fluctuations directly from the mean flow along the free shear layers, and it is removed by $\epsilon _{11}$ and $\varPi _{11}$ , with $\varPi _{11}$ being the dominant sink away from the wall. Fluctuations $\langle u^{\prime 2} \rangle$ generated in the core of the leading-edge shear layer is primarily transported downstream towards the wake while also spreading towards the lower and upper sides of the shear layer. On the upper side, the upward transfer of $\langle u^{\prime 2} \rangle$ is mainly driven by the turbulent transport term, $\langle u^{\prime }u^{\prime }v^{\prime } \rangle$ , balancing the energy loss due to $\varPi _{11}$ in the outer region. On the lower side, the downward transfer of $\langle u^{\prime 2} \rangle$ is induced by both the turbulent transport and advection, ${\langle u^{\prime 2} \rangle } V$ . This downward transfer splits at $x \approx 2.8$ , feeding both upstream recirculation and the downstream vortex shedding regions. In the separation bubble, a singularity point where the $\langle u^{\prime 2} \rangle$ flux lines converge is identified at $(1.60, 0.33)$ , just along the boundary of the backflow region $y_b$ . A second region where $\varPsi _{11}\gt 0$ is observed in the flow impingement area ( $2.5\lt x\lt 4$ ) within the backflow region, resulting from the positive contribution of $\varPi _{11}$ . Along the reattached boundary layer ( $x \gt x_r$ ), $\langle u^{\prime 2} \rangle$ is fed back to the mean flow, $P_{11}\lt 0$ , or dissipated.

In the $\langle v^{\prime 2} \rangle$ budget, $\varPi _{22}$ acts as the primary source within the leading-edge shear layer, resulting in a net increase of $\langle v^{\prime 2} \rangle$ at the expense of $\langle u^{\prime 2} \rangle$ . In the early stages of the turbulent shear layer, negative production $P_{22}$ is outweighed by the positive contribution of $\varPi _{22}$ , except on the lower side of the shear layer, where $\varPsi _{22} \lt 0$ for $0.1\lt x\lt 0.6$ and $y_b \lt y \lt y_{sl}$ . In this area, $\langle v^{\prime 2} \rangle$ transported from the upper side of the shear layer is either returned to the mean flow kinetic energy $VV$ or dissipated by viscosity, $\epsilon _{22}$ . In the separated flow, flux lines curve downward to transfer $\langle v^{\prime 2} \rangle$ towards the wall, where $\varPsi _{22} \lt 0$ . Here, $\varPi _{22}$ is the dominant sink, affected by the wall impermeability constraint. A net rate of gain ( $\varPsi _{22}\gt 0$ ) is obtained along the trailing-edge shear layer and in the wake region where large-scale vertical oscillations occur. This excess of $\langle v^{\prime 2} \rangle$ is carried downstream by advection, $U \langle v^{\prime 2} \rangle$ , and across the wake by the turbulent transport term  $\langle v^{\prime }v^{\prime }v^{\prime } \rangle$ .

The lines of $\langle w^{\prime 2} \rangle$ flux closely follow the mean velocity streamlines, indicating that advection dominates the energy transfer process. The only other significant contributions to the transport of $\langle w^{\prime 2} \rangle$ come from the turbulent fluxes, $\langle w^{\prime }w^{\prime }u^{\prime }\rangle$ and $\langle w^{\prime }w^{\prime }v^{\prime }\rangle$ . The excess of $\langle w^{\prime 2} \rangle$ ( $\varPsi _{33}\gt 0$ ), generated in the core of the leading-edge shear layer through energy exchange with $\langle u^{\prime 2} \rangle$ ( $\varPi _{33}\gt 0$ ), is transported downstream and downward to the lower side of the shear layer, where it dissipates by viscous action, $\epsilon _{33}$ . Near the wall, the net energy gain provided by $\varPi _{33}$ in the impingement region is partly transported downstream along the reattached boundary layer ( $x\gt x_r$ ) and partly flows through the backflow region ( $x\lt x_r$ ), reaching both the leading edge and the upstream part of the shear layer.

In the $\langle u^{\prime } v^{\prime } \rangle$ budget, negative $\varPsi _{12}$ corresponds to a gain where $\langle u^{\prime } v^{\prime } \rangle \lt 0$ , thereby increasing the magnitude of the shear stress. Hence, lines of $\langle u^{\prime } v^{\prime } \rangle$ flux are viewed as carrying $\langle u^{\prime } v^{\prime } \rangle$ from regions where $\varPsi _{12} \lt 0$ to regions where $\varPsi _{12} \gt 0$ , which is opposite to the flux directions shown in figure 17(d). In the core of the leading-edge shear layer, production $P_{12}\lt 0$ exceeds the pressure–strain $\varPi _{12}\gt 0$ , resulting in a net production of $\langle u^{\prime } v^{\prime } \rangle$ . This produced $\langle u^{\prime } v^{\prime } \rangle$ is partially carried downstream by advection and partially outward by pressure and turbulent fluxes. The transport of $\langle u^{\prime } v^{\prime } \rangle$ in the upper side of the shear layer and at the wall is balanced by $\varPi _{12}$ , which serves as a sink. In the recirculating region, $\langle u^{\prime } v^{\prime } \rangle$ flux lines converges to a singularity point slightly downstream of the centre of rotation of the main separation bubble, located at $(1.95, 0.29)$ along $y_b$ .

5. Scalar transport

In addition to representing various physical phenomena, scalar fields provide a complementary perspective on turbulent mixing. This section addresses the turbulent transport of a passive scalar and presents the first detailed scalar transport budget for the rectangular 5 : 1 cylinder, enabling direct comparison with momentum transport mechanisms.

5.1. Turbulent scalar fluxes

At high Reynolds numbers, and thus at high Péclet numbers ( ${\textit{Pe}} = {\textit{Re}} {\textit{Sc}}$ ), the turbulent motion becomes more effective than molecular diffusion in mixing and transporting the scalar. In these conditions, turbulent scalar fluxes, $\langle u_i^{\prime } \theta ^{\prime } \rangle$ , emerge as the primary mechanism for scalar transport.

Figure 18. Scalar fluxes (a,c) $\langle u^{\prime } \theta ^{\prime } \rangle$ and (b,d) $\langle v^{\prime } \theta ^{\prime } \rangle$ at (a,b) ${\textit{Re}}=3000$ and (c,d) $14\,000$ . The cross symbols indicate the locations of maximum and minimum values.

Figure 18 shows the streamwise, $\langle u^{\prime } \theta ^{\prime } \rangle$ , and vertical, $\langle v^{\prime } \theta ^{\prime } \rangle$ , scalar fluxes around the rectangle at ${\textit{Re}}=3000$ and ${\textit{Re}}=14\,000$ . Negative streamwise scalar flux, $\langle u^{\prime } \theta ^{\prime } \rangle$ , is observed in the leading-edge shear layer and near thewall in the reattached boundary layer ( $x\gt x_r$ ). For both ${\textit{Re}}$ cases, ${\langle u^{\prime } \theta ^{\prime } \rangle }_{\textit{min}}$ occurs near the trailing-edge flow separation rather than in the leading-edge shear layer. As highlighted by Cimarelli et al. (2024) in the analysis of scalar variance $\langle \theta ^{\prime } \theta ^{\prime } \rangle$ , while streamwise velocity fluctuations reach their highest values downstream of the early turbulent region of the shear layer, maximum scalar fluctuations arise close to the wall. This is driven by large-scale vortices in the separated flow, which transport low-concentration scalar fluid and impinge on the wall, where $\theta = 1$ . Positive values of $\langle u^{\prime } \theta ^{\prime } \rangle$ appear in the backflow region and the wake separation bubble, with maximum values located near the separation point of the backward boundary layer.

Vertical scalar flux is positive ( ${\langle v^{\prime } \theta ^{\prime } \rangle } \gt 0$ ) around the rectangle, except within the wake separation bubble. High-intensity regions of $\langle v^{\prime } \theta ^{\prime } \rangle$ are found along the leading-edge shear layer, in the reattachment region near the wall, and along the trailing-edge shear layer. Due to the wall impermeability, vertical velocity fluctuations are dampened near the wall, leading ${\langle u^{\prime } \theta ^{\prime } \rangle }_{\textit{max}}$ to occur in the early turbulent region of the shear layer. As Reynolds number increases, the magnitude of scalar fluxes decreases and the gradients near the wall become less steep, indicating reduced turbulence-scalar interaction in this region.

5.2. Mean scalar balance

The role of turbulence in shaping the mean scalar distribution is encapsulated in the mean scalar balance. The conservation equation for the mean passive scalar, $\varTheta$ , in statistically steady and spanwise homogeneous flows is given by

(5.1) \begin{equation} \displaystyle U_i \frac {\partial \varTheta }{\partial x_i} = \frac {\partial }{\partial x_i} \left ( - \left\langle u_i^{\prime } \theta ^{\prime }\right \rangle + \frac {1}{{\textit{Pe}}} \frac {\partial \varTheta }{\partial x_i} \right )\!. \end{equation}

The term on the left-hand side represents the rate of change of $\varTheta$ due to advection, while terms on the right-hand side describe the transport of $\varTheta$ through turbulent and molecular diffusion. In this analysis, the Péclet number ( ${\textit{Pe}}=9940$ ) is sufficiently high to make molecular diffusion negligible in regions away from the wall and outside the laminar section of the separated shear layer. The rate of change in $\varTheta$ due to advection is balanced by turbulent transport.

Figure 19. Scalar transport term $-(\partial {\langle u^{\prime } \theta ^{\prime } \rangle } / \partial x +\partial {\langle v^{\prime } \theta ^{\prime } \rangle } / \partial y)$ with flux lines of $\langle \boldsymbol{u}^{\prime }\theta ^{\prime } \rangle$ at ${\textit{Re}} = 14\,000$ .

Figure 19 shows the turbulent transport term $-\partial \langle u_i^{\prime } \theta ^{\prime } \rangle / \partial x_i$ from (5.1), where positive values indicate a net addition of the mean scalar, while negative values denote scalar loss. The superimposed lines of the scalar flux vector, $\langle \boldsymbol{u}^{\prime }\theta ^{\prime } \rangle$ , display the spatial redistribution of $\varTheta$ , with scalar flux lines originating in regions of scalar removal and terminating in regions of scalar supply.

The distribution of the turbulent transport term varies significantly across different regions of the turbulent flow. In the leading-edge shear layer, turbulent motion carries  $\varTheta$ from the lower side to the upper side, with a more intense transfer observed in the early stages of turbulence development ( $x\lt 1$ ). Scalar removed from the wall region is redistributed into surrounding areas, specifically within the backflow region and along the reattaching boundary layer. In the reattachment region ( $3\lt x\lt 5$ ), the scalar flux is strong enough to carry $\varTheta$ toward the outer edge of the turbulent shear layer. In the upstream portion of the backflow region ( $0\lt x\lt 1$ ), turbulent transport drives $\varTheta$ upward through the turbulent shear layer. Downstream of the rectangle, the scalar flux draws $\varTheta$ from the trailing-edge shear layer and the core of the wake, redistributing it into the outer turbulent region of the flow. This analysis highlights the primary role of the scalar fluxes $\langle u_i^{\prime } \theta ^{\prime } \rangle$ in enhancing scalar mixing across the separated flow and the wake, promoting more uniform scalar distribution throughout these flow regions.

5.3. Scalar flux balances

To better understand how turbulent scalar fluxes are generated, redistributed, and removed, we examine the transport equations for $\langle u_i^{\prime } \theta ^{\prime } \rangle$ . These balances provide a detailed view of the mechanisms governing scalar-velocity correlations, which are responsible for turbulent scalar transport. This section analyses the individual budget terms to identify the dominant processes in separated and reattaching regions.

The balance equations of scalar flux $\langle u_i^{\prime } \theta ^{\prime } \rangle$ for steady conditions are expressed as

(5.2) \begin{equation} \frac {\partial T_{k,i\theta }}{\partial x_{k}} = P_{i\theta } - \epsilon _{i\theta } + \varPi _{i\theta }. \end{equation}

In analogy to (4.7), the left-hand side collects the spatial transport terms arising from advection, pressure, turbulence, and molecular diffusion, represented as the divergence of the flux $T_{k,i\theta }$ ,

(5.3) \begin{equation} T_{k,i\theta } = U_{k} \left\langle u_i^{\prime } \theta ^{\prime } \right\rangle + \langle p^{\prime } \theta ^{\prime } \rangle \delta _{ik} + \left \langle u_{k}^{\prime }u_{i}^{\prime }\theta ^{\prime }\right \rangle - \frac {1}{{\textit{Re}}{\textit{Sc}}} \left \langle u_i^{\prime } \frac {\partial \theta ^{\prime }}{\partial x_k} \right \rangle - \frac {1}{{\textit{Re}}} \left \langle \theta ^{\prime } \frac {\partial u_i^{\prime }}{\partial x_k} \right \rangle \!. \end{equation}

On the right-hand side of (5.2), $P_{i\theta }$ denotes the production term, representing the rate of generation of $\langle u_i^{\prime }\theta ^{\prime } \rangle$ through turbulent interactions with mean scalar and mean velocity gradients. The term $\epsilon _{i\theta }$ represents destruction, which quantifies the mean rate at which $\langle u_i^{\prime }\theta ^{\prime } \rangle$ is reduced by molecular diffusivity of momentum and scalar. This term is zero in isotropic turbulence and is usually neglected in non-isotropic turbulence at high Reynolds numbers or under the assumption of local isotropy (Launder 1976). Finally, $\varPi _{i\theta }$ represents the pressure-scalar gradient correlation and it may be viewed as the counterpart of the pressure-strain term in the $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ budget equations. However, since the sum of the $\varPi _{i\theta }$ terms may not be zero, the pressure-scalar gradient correlation does not necessarily redistribute scalar flux among components. These terms are defined as follows:

(5.4) \begin{align} P_{i\theta } &= -\left\langle u_i^{\prime } u_k^{\prime } \right\rangle \frac {\partial \varTheta }{\partial x_{k}} -\left\langle u_k^{\prime } \theta ^{\prime } \right\rangle \frac {\partial U_{i}}{\partial x_{k}}, \end{align}

(5.5) \begin{align} \epsilon _{i\theta } &= \left ( \frac {1}{{\textit{Re}}} + \frac {1}{{\textit{Re}}{\textit{Pr}}} \right ) \left \langle \frac {\partial u_i^{\prime }}{\partial x_k} \frac {\partial \theta ^{\prime }}{\partial x_k} \right \rangle\! , \end{align}

(5.6) \begin{align} \varPi _{i\theta } &= \left \langle p^{\prime }\, \frac {\partial \theta ^{\prime }}{\partial x_i} \right \rangle\! . \end{align}

Together, these terms form the net source term $\varPsi _{i\theta } = P_{i\theta }-\epsilon _{i\theta }+\varPi _{i\theta }$ representing the net scalar flux that is lost or gained. Similar to $\varPsi _{\textit{ij}}$ , the role of $\varPsi _{i\theta }$ as a source or sink depends on the sign of $\langle u_i^{\prime }\theta ^{\prime } \rangle$ . Each term in (5.2) is indexed to indicate the direction of the scalar flux under consideration ( $i=1$ for $\langle u^{\prime } \theta ^{\prime } \rangle$ or $i=2$ for $\langle v^{\prime } \theta ^{\prime } \rangle$ ).

Figure 20. Production terms of scalar flux (a) $P_{1\theta }$ and (b) $P_{2\theta }$ at ${\textit{Re}}=14\,000$ . The cross symbols mark the locations of the minimum value of $P_{1\theta }$ and maximum value of $P_{2\theta }$ . The green line denotes the dividing streamline originating from the leading edge, $y_\psi$ .

Figure 20 shows regions of positive and negative production rates of scalar flux, $P_{i\theta }$ , for the flow at ${\textit{Re}}=14\,000$ . The scalar flux $\langle u^{\prime } \theta ^{\prime } \rangle$ is primarily generated within the leading-edge shear layer. In this area, negative values of $P_{1\theta }$ act as a source since ${\langle u^{\prime } \theta ^{\prime } \rangle } \lt 0$ . Near the wall, in the backflow region, positive values of $P_{1\theta }$ correspond to zones where ${\langle u^{\prime } \theta ^{\prime } \rangle } \gt 0$ and hence are sources of $\langle u^{\prime } \theta ^{\prime } \rangle$ . The minimum value $P_{1\theta ,\,{min}}=-0.67$ is found at $(x, y) = (0.40, 0.30)$ . The term $P_{2\theta }$ acts as a source for $\langle v^{\prime } \theta ^{\prime } \rangle$ above the rectangle, with peak intensity along the shear-layer centreline, with $P_{2\theta ,\,{max}}=0.24$ at $(x, y) = (0.49, 0.34)$ , and in the reattachment region ( $2.2\lt x\lt 5$ ).

Figure 21. Production terms of (a) $\langle u^{\prime } \theta ^{\prime } \rangle$ and (b) $\langle v^{\prime } \theta ^{\prime } \rangle$ evaluated along the reattaching streamline $y_\psi$ : solid line, $P_{i\theta }$ ; dashed line, $P_{i\theta ,a}$ ; and dotted line, $P_{i\theta ,b}$ .

To elucidate the production mechanism of scalar fluxes, the individual contributions to $P_{1\theta }$ and $P_{2\theta }$ are examined. Each production term is expressed as the sum of two components, representing the effects of mean scalar gradients and mean velocity gradients interacting with turbulent motions, as follows:

(5.7) \begin{equation} P_{i\theta } = \underbrace {-\left\langle u_i^{\prime } u_k^{\prime } \right\rangle \frac {\partial \varTheta }{\partial x_{k}}}_{P_{i\theta ,a}} \underbrace {-\left\langle \theta ^{\prime } u_k^{\prime } \right\rangle \frac {\partial U_{i}}{\partial x_{k}}}_{P_{i\theta ,b}}. \end{equation}

Figure 21 illustrates the two terms $P_{i\theta ,a}$ and $P_{i\theta ,b}$ along with the total production rate $P_{i\theta }$ , evaluated along the dividing streamline $y_\psi$ . The production term $P_{1\theta }$ is evenly affected by the mean scalar and velocity gradients ( $P_{1\theta ,a} \approx P_{2\theta ,b}$ ) along the free shear layer ( $x\lt 3$ ) and in the reattached flow region ( $x\gt 4$ ). In contrast, the dominant contribution to $P_{2\theta }$ arises from $P_{2\theta ,a}$ , driven by the interaction of the wall-normal scalar gradient with vertical Reynolds stress, $-\langle v^{\prime 2} \rangle \partial \varTheta /\partial y$ . The contribution due to the mean velocity gradient term, $P_{2\theta ,b}$ , results in a loss of $\langle v^{\prime } \theta ^{\prime } \rangle$ within the shear layer and becomes negligible in the reattached region. This opposing effect reduces the magnitude of $P_{2\theta }$ .

Figure 22. Destruction terms of scalar flux (a) $\epsilon _{1\theta }$ and (b) $\epsilon _{2\theta }$ at ${\textit{Re}}=14\,000$ . The cross symbols mark the locations of the minimum values of $\epsilon _{i\theta }$ . (c) Detailed view of the distribution of $\epsilon _{2\theta }$ .

Assuming the Reynolds number is high enough, the small scales can be considered locally isotropic, making the destruction terms $\epsilon _{i\theta }$ negligible compared with $P_{i\theta }$ and  $\varPi _{i\theta }$ . However, this assumption does not fully account for the energy transfer to higher wavenumbers, a process that requires some level of destruction (Younis, Speziale & Clark 1995). Furthermore, DNS results from Rogers, Mansour & Reynolds (1989), obtained for a fully developed turbulent channel flow, indicate that $\epsilon _{i\theta }$ can act as a production term under certain conditions. Figure 22 shows the destruction terms $\epsilon _{i\theta }$ for the case under consideration. These terms are approximately an order of magnitude smaller than $P_{i\theta }$ and $\varPi _{i\theta }$ , and the colourbar limits have been adjusted accordingly for clarity. In the streamwise balance, $\epsilon _{1\theta }$ is a sink for $\langle u^{\prime } \theta ^{\prime } \rangle$ , since $-\epsilon _{1\theta }\gt 0$ (or $-\epsilon _{1\theta }\lt 0$ ) where ${\langle u^{\prime } \theta ^{\prime } \rangle }\lt 0$ (or ${\langle u^{\prime } \theta ^{\prime } \rangle }\gt 0$ ). The minimum value, $\epsilon _{1\theta ,\,{min}} = -0.015$ , occurs along the shear-layer centreline $y_{sl}$ at $(0.51, 0.35)$ . For the vertical balance, $\epsilon _{2\theta }$ predominantly leads to a loss of $\langle v^{\prime } \theta ^{\prime } \rangle$ around the rectangle, except for a localised zone along $y_{sl}$ ( $0.2\lt x\lt 0.6$ ) where $\epsilon _{2\theta }\lt 0$ is observed, see figure 22(c) for the detailed view. This negative region, with $\epsilon _{2\theta ,\,{\textit{min}}} = -0.0087$ at $(0.37, 0.28)$ , corresponds to positive values of $\langle v^{\prime } \theta ^{\prime } \rangle$ , indicating a production of scalar flux due to viscous effects. It is worth noting that this localised gain in scalar flux contributes significantly less than other sink and source terms in the same regions. For instance, at the location of $\epsilon _{2\theta ,\,{\textit{min}}}$ , $P_{2\theta }=0.103$ .

With the viscous destruction negligible, the pressure-scalar gradient correlations $\varPi _{i\theta }$ provide the mechanism that limits the growth of the scalar fluxes by counteracting the rate of production of $\langle u_i^{\prime }\theta ^{\prime } \rangle$ . Figure 23 shows the distribution of $\varPi _{i\theta }$ . The streamwise term $\varPi _{1\theta }$ is found to oppose to $P_{1\theta }$ discussed previously, serving as a sink of $\langle u^{\prime } \theta ^{\prime } \rangle$ along the free shear layer and in the reattachment zone. The positive peak $P_{1\theta ,\,{max}}$ is found at $(0.43, 0.31)$ . Regarding $\varPi _{2\theta }$ , the centreline $y_{sl}$ approximately divides a gain region ( $\varPi _{2\theta }\gt 0$ ) on the upper side of the shear layer from a loss region ( $\varPi _{2\theta }\lt 0$ ) on the lower side. In the latter, the negative peak occurs at $(0.43, 0.31)$ , with $P_{2\theta ,\,{\textit{min}}} = -0.34$ .

Figure 23. Pressure-scalar gradient correlation of scalar flux (a) $\varPi _{1\theta }$ and (b) $\varPi _{2\theta }$ at ${\textit{Re}}=14\,000$ . The cross symbols mark the locations of the maximum value of $\varPi _{1\theta }$ and minimum value of $\varPi _{2\theta }$ .

Figure 24. Distribution of the net source terms of the scalar flux, $\varPsi _{i\theta } = P_{i\theta } + \varPi _{i\theta } - \epsilon _{i\theta }$ overlaid with flux lines of $\boldsymbol{T_{i\theta }}$ : (a) $\varPsi _{1\theta }$ and (b) $\varPsi _{2\theta }$ .

The net sources of the streamwise and vertical scalar fluxes, $\varPsi _{1\theta }$ and $\varPsi _{2\theta }$ , are analysed along with lines of the spatial flux vector $\boldsymbol{T_{i\theta }}$ shown in figure 24. Regions where $\varPsi _{i\theta }\gt 0$ ( $\varPsi _{i\theta }\lt 0$ ) indicate sources for the scalar flux if $\langle u_i^{\prime } \theta ^{\prime } \rangle \gt 0$ ( $\langle u_i^{\prime } \theta ^{\prime } \rangle \lt 0$ ), corresponding to areas of net increase in scalar flux magnitude. In the shear layer and wake, where $\langle u^{\prime } \theta ^{\prime } \rangle$ is predominantly negative, the transport lines of scalar flux in figure 24(a) indicate flux $\langle u^{\prime } \theta ^{\prime } \rangle$ being carried from source regions ( $\varPsi _{1\theta }\lt 0$ ) on the upper side of the shear layer to the lower side and into the wake, moving opposite to the direction indicated by the arrows in the figure. Within the backflow region, the excess of $\langle u^{\prime } \theta ^{\prime } \rangle$ generated near the wall is transported upstream towards the leading edge of the rectangle and upward, beginning at $x \approx 1.5$ .

For the vertical flux budget, the flux $\langle v^{\prime } \theta ^{\prime } \rangle$ gained in the upper side of the shear layer $\varPsi _{2\theta } \gt 0$ is redistributed downward towards the the lower side and into the wake, where it is ultimately removed. In the near-wall region, in particular within $2\lt x\lt 5$ , a strong source of $\langle v^{\prime } \theta ^{\prime } \rangle$ is evident, associated with the impingement of turbulent vortices formed in the shear layer. The flux $\langle v^{\prime } \theta ^{\prime } \rangle$ generated near the wall is subsequently transported to the lower side of the shear layer, where it is drained ( $\varPsi _{2\theta } \lt 0$ ).

6. Assessment of eddy viscosity and eddy diffusivity models

Finally, this section evaluates common turbulence modelling assumptions by comparing DNS results with eddy-viscosity and eddy-diffusivity closures. This a priori assessment highlights and quantifies the limitations of common models in capturing complex features of separating and reattaching flows, by quantifying and localising their departures from the DNS reference.

In eddy-viscosity turbulence models, the Reynolds stresses $\langle u_i u_j \rangle$ and scalar fluxes $\langle u_i \theta \rangle$ , which appear as unknowns in the Reynolds-averaged Navier–Stokes equations (4.3) and mean passive scalar conservation equation (5.1), are modelled based on the eddy-viscosity and gradient-diffusion hypotheses (Boussinesq’s hypothesis). These approaches express the Reynolds stresses and scalar fluxes in terms of the mean velocity and scalar gradients, respectively, as follows:

(6.1) \begin{align} \big\langle u_{i}^{\prime } u_{j}^{\prime } \big\rangle - \frac {2}{3}k\delta _{\textit{ij}} &= - 2 \nu _{T} S_{\textit{ij}}, \end{align}

(6.2) \begin{align} \left\langle u_i^{\prime } \theta ^{\prime } \right\rangle &= - \alpha _{T} \frac {\partial \varTheta }{\partial x_{i}}, \end{align}

where $k = ({1}/{2}) \langle u_i^{\prime } u_i^{\prime } \rangle$ is the turbulent kinetic energy, and $\nu _{T}$ and $\alpha _{T}$ are scalar coefficients known as the eddy viscosity and eddy diffusivity, respectively. The eddy-viscosity hypothesis implies that the deviatoric Reynolds stress, $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle - ({2}/{3})k\delta _{\textit{ij}}$ , is proportional to the mean rate of strain, whereas the gradient-diffusion hypothesis postulates that the turbulent scalar flux is aligned with the gradient of the mean scalar field.

These hypotheses provide a practical closure for the governing equations, making them widely adopted in engineering applications due to their simplicity and ease of implementation. However, their limitations are well documented in the literature (Spalart 2000). In particular, their predictions often lack accuracy in complex flows, such as those involving separation and reattachment, where non-equilibrium effects and anisotropy play significant roles. For instance, in separated turbulent flows, eddy-viscosity models may predict the onset of separation adequately, but may fail in describing the flow downstream of the separation line. This inadequacy can have significant repercussions in the prediction of the aerodynamic performance of bluff bodies, as highlighted by Menter & Kuntz (2004) and Rumsey (2009). The main shortcomings include: under-prediction of eddy viscosity in the shear layer of separation bubbles, leading to insufficient mixing and delayed reattachment; inaccurate modelling of flow recovery after reattachment (Johnson, Menter & Rumsey 1994); poor representation of laminar-to-turbulent transition mechanisms (Crivellini, Ghidoni & Noventa 2023); difficulties in capturing backflow regions (Patrick 1987); and unreliable predictions of heat transfer distribution in the presence of significant adverse pressure gradients (Barnett & Carter 1986).

If the eddy-viscosity and gradient-diffusion hypotheses are accepted as adequate approximations, determining suitable formulations for $\nu _{T}$ and $\alpha _{T}$ becomes critical. Although DNS does not incorporate an explicit eddy viscosity, it provides high-fidelity data that can be used to compute an equivalent eddy viscosity for benchmarking turbulence models. In this study, we adopt the following definition, also used in the work by Spalart et al. (2015):

(6.3) \begin{equation} \nu _{T} = -\frac {\big\langle u_{i}^{\prime } u_{j}^{\prime }\big \rangle S_{\textit{ij}}}{2 S_{\textit{ij}} S_{\textit{ij}}}. \end{equation}

where $\nu _{T}$ can be physically interpreted as the optimal value required to capture the correct turbulent kinetic energy production, $-\langle u_{i}^{\prime } u_{j}^{\prime } \rangle S_{\textit{ij}}$ , which is critical for the consistency of the mean flow energy balance.

The eddy diffusivity, $\alpha _{T}$ , is commonly related to the eddy viscosity via the turbulent Schmidt number, $\sigma _{T}$ :

(6.4) \begin{equation} \alpha _{T} = \frac {\nu _{T}}{\sigma _{T}}, \end{equation}

where $\sigma _{T}$ quantifies the relative efficiency of turbulent transport of momentum and scalar quantities. Typical values for $\sigma _{T}$ range from $0.7$ to $0.9$ , although lower values (e.g. $\sigma _{T} \approx 0.4$ ) are occasionally used in specific configurations (Chambers, Antonia & Fulachier 1985). Experimental evidence suggests that assuming a constant $\sigma _{T}$ can be overly simplistic, particularly in flows with heat transfer and adverse pressure gradients (Barnett & Carter 1986).

To derive an effective eddy diffusivity from DNS data, the following definition is employed:

(6.5) \begin{equation} \alpha _{T} = -\frac {\displaystyle \big\langle u_i^{\prime } \theta ^{\prime } \big\rangle \frac {\displaystyle \partial \varTheta }{\partial x_{i}}}{\displaystyle \frac {\partial \varTheta }{\partial x_{i}} \, \frac {\partial \varTheta }{\partial x_{i}}}. \end{equation}

This formulation ensures the same production rate of scalar variance $\langle \theta ^{\prime } \theta ^{\prime } \rangle$ , given by $-2 \langle u_i^{\prime } \theta ^{\prime } \rangle \partial \varTheta /\partial x_{i}$ . Consequently, a spatially varying and flow dependent turbulent Schmidt number $\sigma _{T}$ can be computed based on the effective $\nu _{T}$ and $\alpha _{T}$ as

(6.6) \begin{equation} \sigma _{T} = \frac {\nu _{T}}{\alpha _{T}}. \end{equation}

Closures for vector quantities, such as turbulent scalar fluxes in (6.2), require matching of both magnitude and direction between modelled and actual fluxes. For tensor quantities, such as in (6.1), proportionality implies that the tensors act identically on their shared eigenvectors, differing only by a scaling factor. Specifically, if $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle - ({2}/{3})k\delta _{\textit{ij}}$ and $-S_{\textit{ij}}$ are proportional, their eigenvectors $\boldsymbol{e}_{i}^{R}$ and $\boldsymbol{e}_{i}^{S}$ should be equal, and their eigenvalues $\lambda _{i}^{R}$ and $\lambda _{i}^{S}$ should be related by a proportionality constant $C_{\lambda }$ , i.e.

(6.7) \begin{equation} \frac {\lambda _{1}^{R}}{\lambda _{1}^{S}} = \frac {\lambda _{2}^{R}}{\lambda _{2}^{S}} = C_{\lambda }. \end{equation}

Figure 25. Alignment between the Reynolds stress tensor $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ and the mean rate-of-strain tensor $S_{\textit{ij}}$ . (a) Distribution of the angle $\beta$ between the principal axes of $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ and $S_{\textit{ij}}$ . (b) Visualisation of principal axes of $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ (black arrows) and $S_{\textit{ij}}$ (red arrows). (c) Probability density function of $\beta$ in the turbulent flow region. The dashed line indicates the abscissa $\beta = \pi /6$ , corresponding to a validity threshold. (d) Ratio of the tensors eigenvalues $(\lambda _{1}^{R} \lambda _{2}^{S}) / (\lambda _{1}^{S} \lambda _{2}^{R})$ .

To assess the validity of the eddy-viscosity hypothesis, a check of the relative alignment between the principal axes of the Reynolds stress tensor $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ and those of the mean rate-of-strain tensor $S_{\textit{ij}}$ is displayed in figure 25(a–c). Note the deviatoric component of a tensor ( $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle - ({2}/{3})k\delta _{\textit{ij}}$ ) shares the same principal axes as the full tensor ( $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ ). Additionally, for two symmetric tensors, such as $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ and $S_{\textit{ij}}$ , their principal axes form orthogonal pairs, ensuring that the angle between the corresponding principal axes are the same. Since eigenvectors are defined up to a scalar multiple and their sign is therefore arbitrary, a consistent sign convention is adopted: each eigenvector is flipped, if needed, to produce a positive dot product with a fixed reference direction, i.e. the $x$ -axis. This ensures consistency and avoids ambiguities in alignment analysis across different computations. The angle $\beta$ is then defined as the angle between the principal axes of the tensors $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ and $S_{\textit{ij}}$ , given by

(6.8) \begin{equation} \beta = \arccos \frac {\big|\boldsymbol{e}_{1}^{R} \boldsymbol{\cdot }\boldsymbol{e}_{1}^{S}\big|}{\big|\boldsymbol{e}_{1}^{R}\big| \, \big|\boldsymbol{e}_{1}^{S}\big|}. \end{equation}

Figure 25(a) shows the spatial distribution of $\beta$ , where values near zero indicate alignment between principal directions of $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle$ and $S_{\textit{ij}}$ , consistent with the eddy-viscosity hypothesis. Conversely, larger $\beta$ values indicate misalignment and reduced validity of the hypothesis. As a rule-of-thumb, angles smaller than $\pi /6$ are used here to indicate approximate alignment, as suggested by Schmitt (2007). This corresponds to $\cos (\beta ) \gt 0.866$ , implying a limited deviation from perfect alignment ( $\cos (\beta ) = 1$ ). While this threshold has no strict physical meaning, it offers a practical reference for interpreting alignment in turbulence modelling frameworks. As shown in figure 25(b), small angles are observed in the core of the shear layer and in the wake. The alignment deteriorates significantly near the wall and in the backflow region. This supports the suggestion by Simpson et al. (1981) that the Reynolds shear stress in the backflow region should be modelled based on the turbulence structure rather than local mean velocity gradients. Indeed, the mean velocity profiles in the backflow result from time-averaging the large turbulent fluctuations and are not indicative of the underlying causes of the turbulence. To provide unbiased results, the probability density function (p.d.f.) of the angle $\beta$ is measured in the turbulent region ( $0 \lt x \lt 7.5$ and $-0.5 \lt y \lt y_{\varOmega }$ ), as shown in figure 25(c). The p.d.f. highlights a dominant peak at $\beta = 0$ , confirming strong alignment in a large portion of the subdomain considered, with additional peaks at $\beta \simeq 0.35$ and $\pi /4$ , reflecting misaligned regions in the near-wall and shear layer zones. In addition to evaluating eigenvector alignment, the proportional relationship between the two tensors is assessed by analysing the ratio of their eigenvalues. Figure 25(d) shows the ratio $(\lambda _{1}^{R} \lambda _{2}^{S}) / (\lambda _{1}^{S} \lambda _{2}^{R})$ , which should equal unity if the proportionality condition described by (6.7) holds. The observed deviations from unity across the flow domain indicate that $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle - ({2}/{3})k\delta _{\textit{ij}}$ and $-S_{\textit{ij}}$ are not strictly proportional, even in regions where eigenvector alignment is satisfactory. This highlights that while alignment trends partially support the hypothesis, proportionality breaks down in large portions of the separated flow region, suggesting limitations of the eddy-viscosity assumption in capturing complex flow dynamics.

Figure 26. Alignment between scalar flux $\langle u_i^{\prime } \theta ^{\prime } \rangle$ and mean scalar gradient $-\partial \varTheta /\partial x_{i}$ . (a) Distribution of angle $\gamma$ between $\langle u_i^{\prime } \theta ^{\prime } \rangle$ and $-\partial \varTheta /\partial x_{i}$ . (b) Orientation of $\langle u_i^{\prime } \theta ^{\prime } \rangle$ (black arrows) and $-\partial \varTheta /\partial x_{i}$ (red arrows). (c) Probability density function of $\gamma$ in the turbulent flow region. The dashed line indicates the abscissa $\beta = \pi /6$ , corresponding to a validity threshold.

The gradient-diffusion hypothesis is evaluated by comparing the orientation of the scalar flux $\langle u_i^{\prime } \theta ^{\prime } \rangle$ relative to the mean scalar gradient $-\partial \varTheta / \partial x_i$ , shown in figure 26. The alignment between these two vectors is quantified by the angle $\gamma$ , computed as

(6.9) \begin{equation} \gamma = \arccos \frac {-\langle \boldsymbol{u}^{\prime } \theta ^{\prime } \rangle \boldsymbol{\cdot }\boldsymbol{\nabla }\varTheta }{|\langle \boldsymbol{u}^{\prime } \theta ^{\prime } \rangle | \, |\boldsymbol{\nabla }\varTheta |}. \end{equation}

Smaller values of $\gamma$ indicate better agreement with the gradient-diffusion hypothesis. Figure 26(a) reveals that $\gamma$ rarely approaches zero, even in zones with moderate alignment. In the free shear layer, large values of $\gamma$ occur, while downstream of the rectangular region ( $x \gt 6$ ), $\langle u_i^{\prime } \theta ^{\prime } \rangle$ and $-\partial \varTheta / \partial x_i$ exhibit nearly opposite directions, indicating significant misalignment. This misalignment is also evident near the wall in figure 26(b) and indicates significant deviations from the gradient-diffusion hypothesis. The p.d.f. of $\gamma$ , shown in figure 26(c), exhibits a broader distribution compared with $\beta$ , with a lower peak at $\gamma \simeq \pi /3$ . Non-zero values of the p.d.f. persist over a wide range of angles higher than $\pi /6$ , suggesting the inadequacy of the hypothesis over extensive regions of the flow.

Figure 27. Contours of (a) eddy viscosity $\nu _{T}$ , (b) eddy diffusivity $\alpha _{T}$ and (c) turbulent Schmidt number $\sigma _{T}$ computed from DNS data at ${\textit{Re}}=14\,000$ .

Figure 27 presents the distribution of eddy viscosity $\nu _{T}$ , eddy diffusivity $\alpha _{T}$ and the turbulent Schmidt number $\sigma _{T}$ , derived a priori from DNS data using (6.3), (6.5) and (6.6), respectively. High values of eddy viscosity $\nu _{T}$ (figure 27 a) are observed in the downstream portion of the shear layer, with significant growth in the wake and maxima for $x\gt 6$ . Non-negligible positive values also appear in the backflow region. Conversely, $\nu _{T}$ assumes negative values in the near-wall region of the reattachment zone ( $2.3 \lt x \lt 5$ ) and within the separation bubble in the wake. These negative values indicate reversed momentum transport with respect to that aligned with $S_{\textit{ij}}$ . The eddy diffusivity $\alpha _{T}$ (figure 27 b) shows high values above the reattachment zone and in the wake for $y\gt 0$ . In the wake, a region of negative $\alpha _{T}$ appears, corresponding with the area of intense vertical velocity fluctuations induced by the shedding of large-scale vortices from the trailing edge. The turbulent Schmidt number $\sigma _{T}$ (figure 27 c) exhibits significant spatial variability around the rectangle, highlighting the limitations of assuming $\sigma _{T}$ as constant in the flow. Regions with high $\sigma _{T}$ non-uniformity are found in areas of complex flow, such as zones of flow recirculation and the wake. Simple models, such as the Reynolds analogy ( $\sigma _{T} = 1$ ), fail to accurately represent the detailed transport behaviour in this flow configuration.

In conclusion, flows involving separation and passive scalar transport, as in the present case, demand advanced turbulence models for accurate representation in RANS frameworks. Specifically, these models must independently account for turbulent momentum flux and scalar flux components to capture the intricate transport mechanisms and spatial variability characteristic of separated turbulent flows.

7. Conclusions

The turbulent flow and scalar transport around a 5 : 1 rectangular cylinder involve several complex phenomena, including fixed point laminar separation, shear layer formation and evolution, flow reattachment, recirculating bubble formation, boundary layer development post-reattachment, and wake formation. The general picture of these flow characteristics is gained through a detailed statistical analysis.

Observation of the principal mean strain-rate along the shear layer enables the identification of distinct $S_{\lambda }$ decay laws in the transitional and turbulent regimes. This allows for the location of turbulence onset in the shear layer at different Reynolds numbers. Also, Reynolds stress distribution is investigated in a Reynolds number effect perspective: it is the $x$ -distribution of Reynolds stress integrals, rather than their intensity, that varies with ${\textit{Re}}$ . The mean momentum balance is analysed in both the Cartesian frame and in a curvilinear coordinate system, where coordinates are aligned and orthogonal to the mean streamlines. The momentum balance along the streamline dividing the recirculating bubble from the outer flow highlights the influence of pressure and the divergence of Reynolds stress in shaping the streamlines. The two terms induce both the initial curvature and flow reattachment, thereby defining the boundary of the separation bubble. The Reynolds stress budgets provide a detailed description of the sources, transport paths, redistribution among different components and dissipation of each tensor component. The leading-edge shear layer is the location where these phenomena are most intense. There, the production of $\langle u^{\prime 2} \rangle$ and $\langle u^{\prime } v^{\prime } \rangle$ occurs at the expense of the mean flow. Turbulent energy is partially lost through a corresponding sink of $\langle v^{\prime 2} \rangle$ , feeding the mean vertical flow and partially dissipated in the same region. Additionally, pressure–strain interactions redistribute turbulent fluctuations, transferring energy from the dominant streamwise component $\langle u^{\prime 2} \rangle$ to the cross-flow components $\langle v^{\prime 2} \rangle$ and $\langle w^{\prime 2} \rangle$ . Paths of scalar transport indicate that $\varTheta$ is drawn by turbulent fluxes from the wall through the shear layer and up to the external flow region, alternating between areas of mean scalar removal and mean scalar release. Unlike velocity fluctuations, which reach large values only in the free-shear layer, high-magnitude scalar fluctuations also arise close to the wall. This is driven by large-scale vortices in the separated flow, which transport low-concentration scalar fluid and impinge on the wall, where $\varTheta = 1$ . Budgets of the individual turbulent scalar fluxes show that the net source terms of $\langle u^{\prime } \theta ^{\prime } \rangle$ and $\langle v^{\prime } \theta ^{\prime } \rangle$ are largest in the leading-edge shear layer ( $0 \lt x \lt 1$ ) and at the reattached boundary layer ( $2 \lt x \lt 5$ ), due again to the impingement of turbulent vortices formed in the shear layer.

Modelling flows of this complexity using eddy viscosity and eddy diffusivity approaches may exhibit inaccuracies. An a priori analysis is conducted to examine the misalignment between the mean strain-rate tensor and the deviatoric component of the Reynolds stress tensor, and between the mean scalar gradient and the turbulent scalar fluxes. Although the principal axes of $S_{\textit{ij}}$ and $\langle u_{i}^{\prime } u_{j}^{\prime } \rangle - ({2}/{3})k\delta _{\textit{ij}}$ exhibit good alignment, with the misalignment angle $\beta$ generally below $\pi /6$ above the cylinder, to be proportional the two tensors must have eigenvalues related by a proportionality constant. This condition is not satisfied in the shear layer or, more broadly, in the separated flow region. The alignment between $-\partial \varTheta /\partial x_{i}$ and $\langle u_i^{\prime } \theta ^{\prime } \rangle$ is generally poor around the body. The p.d.f. of the misalignment angle $\gamma$ shows that most values lie well beyond the threshold for validity, with a peak near $\pi /3$ . There are regions where the a priori turbulent stress and turbulent scalar flux point in the opposite direction compared with fully resolved, direct numerical simulation results. Due to the variability of the a priori $\nu _T$ and $\alpha _T$ distributions, the turbulent Schmidt number varies widely, ranging from $-2$ to $2$ . Therefore, turbulence models are required to go beyond the eddy diffusivity and eddy viscosity concepts of most present methods for accurately predicting flows involving separation and reattachment phenomena.