The turbulent flow over the BARC rectangular cylinder: a DNS study

A Direct Numerical Simulation (DNS) of the incompressible flow around a rectangular cylinder with chord-to-thickness ratio 5:1 (also known as the BARC benchmark) is presented. The work replicates the first DNS of this kind recently presented by Cimarelli et al (2018), and intends to contribute to a solid numerical benchmark, albeit at a relatively low value of the Reynolds number. The study differentiates from previous work by using an in-house finite-differences solver instead of the finite-volumes toolbox OpenFOAM, and by employing finer spatial discretization and longer temporal average. The main features of the flow are described, and quantitative differences with the existing results are highlighted. The complete set of terms appearing in the budget equation for the components of the Reynolds stress tensor is provided for the first time. The different regions of the flow where production, redistribution and dissipation of each component take place are identified, and the anisotropic and inhomogeneous nature of the flow is discussed. Such information is valuable for the verification and fine-tuning of turbulence models in this complex separating and reattaching flow.


Introduction
The flow around bluff bodies with sharp corners is interesting for both fundamental research and industrial applications, particularly in the field of vortex-induced oscillations [41]. In fact, several types of structures (buildings, bridges, pylons etc) often present such cross-sectional shapes [38].
The simplest prototype of such bodies is the rectangular cylinder. It sports a simple geometry, yet the flow around it is rich with features that are found in flows around bodies of more complex shape: a corner-induced separation, a shear layer which at a certain point becomes unstable, a detached boundary layer that may reattach, several recirculating regions and a large wake. By varying the aspect ratio A, the rectangular cylinder spans the entire set of blunt bodies, from a zerothickness flat plate normal to the flow for A = 0, to a square cylinder for A = 1, and finally to a flat plate parallel to the flow for A → ∞.
Already at low values of the Reynolds number Re, the main flow features depend on A . For small aspect ratios, say A ≤ 2, the body is not long enough for the flow to reattach after the leading-edge separation. For 2 ≤ A ≤ 3 the flow does reattach along the top and bottom sides, but reattachment is intermittent, and vortex shedding still occurs from the leading-edge corners only. For even larger A, reattachment becomes permanent, creating a large recirculating region sometimes referred to as primary vortex, and vortex shedding occurs from both the leadingand trailing-edge corners, while the main flow features keep changing with the aspect ratio. For large enough Re, transition to turbulence complicates the matter further: the large scales associated with the flow instabilities coexist and nonlinearly interact with the smaller scales associated to the turbulent motions. This variability is accompanied by a significant scatter of experimental and numerical data. This is the main reason why a benchmark on the aerodynamics of a rectangular cylinder with A = 5 (Benchmark on the Aerodynamics of a Rectangular 5:1 Cylinder, or BARC) has been established [2]. Goal of the BARC benchmark is to set the standards for both numerical simulations and experiments, and to arrive at a quantitatively accurate description of the main patterns of the flow, e.g. the size of the primary vortex, the shedding frequency, the root-meansquare value of fluctuations of the vertical force. Unfortunately, as discussed by Bruno et al. in [5], even with A fixed, a large variability remains, which still affects the prediction of the main features of the mean flow. Besides the effect of the Reynolds number, this is due to the strong sensitivity of the BARC flow to several aspects of the experimental as well as the numerical studies. For experiments, these range from test conditions and body shape to measurement inaccuracies; for numerical simulations, they include RANS and LES turbulence modelling, discretisation choices, the numerical method itself and the computational procedure. The dispersion of BARC experimental data was considered in Ref. [24], where uncertainties on the angle of incidence, the free-stream longitudinal turbulence intensity and the free-stream turbulence lengthscale were studied via a probabilistic method and two-dimensional URANS simulations; such uncertainties were found to not fully explain the scatter of results, with the turbulence model playing a major role. Similarly Ref. [25] presented a stochastic analysis of the sensitivity of LES results to grid resolution and details of the model filter, finding that a finer grid translates into a significantly smaller primary vortex. More recently, Ref. [34] performed a sensitivity analysis of LES simulations of the BARC flow to the rounding of the leading-edge corners, to the aim of explaining the discrepancy between well-resolved numerical simulations and experiments. They observed that introducing a very small radius of curvature is sufficient to enhance the agreement between numerical and experimental data.
As an example, we summarise below the information available for some quantities of interest, namely the frequency (expressed by the non-dimensional Strouhal number St) of the vortex shedding and the length L 1 (expressed in terms of the body height D) of the main recirculating region created by the leading-edge shear layer. Further information can be found in Ref. [5]. In [26] Matsumoto and cowork-ers studied the BARC flow experimentally at Re ≈ 10 5 , with both smooth and turbulent incoming flows. The dominant shedding frequency was St = 0.132 for the smooth inflow and St = 0.197 for the turbulent inflow, with the primary vortex length varying from L 1 = 1.875 to L 1 = 4.375. Bruno et al. in [4] performed a LES simulation at Re = 4 · 10 4 finding St ≈ 0.11 and L 1 ≈ 4.68. Mannini et al. [21] performed URANS simulations at Re = 1 · 10 5 using a slightly modified Spalart-Allmaras turbulence model and a Linearised Explicit Algebraic model coupled with the standard k − ω model. Using the first model they found St ≈ 0.098 with no flow reattachment; with the second model, instead, they found St ≈ 0.105 and L 1 ≈ 4.65. Mannini et al. [22] performed a DES study at Re = 26400 focusing on the effect of the spanwise dimension of the domain on the main properties of the flow; they report St ≈ 0.1 and L 1 ≈ 4.75. Bruno et al. [3] studied the combined effect of spanwise discretisation and spanwise domain size with LES. For an increasing spanwise resolution L 1 decreases whereas St remains practically constant; increasing the spanwise domain size, instead, affects neither St nor L 1 . Patruno et al. [32] studied with LES and URANS the effect of small angles of incidence at Re ≈ 10 4 . At zero incidence they found a non-symmetric mean flow with LES, with L 1 = 4.01 and L 1 = 4.1 in the upper and lower cylinder sides and St = 0.132. Their URANS simulation, instead, using the k − ω − SST turbulence model yielded a symmetric mean flow with L 1 = 4.26 and St = 0.121. Ricci et al. [33] performed LES simulations at Re ≈ 5.5 · 10 4 to study the turbulent inflow conditions. They found that a higher level of incoming turbulence corresponds to a higher curvature of the shear layer and therefore to a shorter primary vortex, with a related upstream shift of the secondary vortex. A similar study was also carried out by Mannini et al. [20] experimentally, by varying the inflow conditions and the angle of incidence in the range 10 4 < Re < 10 5 . For zero incidence they found St ≈ 0.115 which weakly increases with Re and with the intensity of the free stream turbulence. Moore et al. [28] experimentally found for Re = 1.34 · 10 4 St ≈ 0.1114 and L 1 ≈ 4.4; moreover, they found that St does not change with the Reynolds number in the range 1.34 · 10 4 ≤ Re ≤ 1.18 · 10 5 unlike L 1 which has a decreasing trend. Lastly, Cimarelli et al. [7] investigated via high-order implicit LES the influence of geometrical characteristics of the body, as the sharp leading-edge corners, the presence of separation at the trailing edge and the coupling between the two sides of the plate. At Re = 3000 they report St = 0.14 and L 1 ≈ 4.05.
In such scenario of highly scattered experimental and numerical information, a remarkable achievement was the recent first Direct Numerical Simulation (DNS) of the BARC flow in the turbulent regime. It was presented by Cimarelli, Leonforte and Angeli [10], who employed the finite-volumes OpenFOAM toolbox [40], and was then used in the derivative works [8,9]. Albeit at a relatively low value of Re, these data represent a key step towards building a robust information set for the BARC flow, since DNS has the ability to remove from the picture the uncertainties related to turbulence modeling. It remains desirable, however, to assess the robustness of their results with respect to discretization, numerical method and computational procedures. This is one of the goals of the present work, which replicates this DNS study but employs a different solver (an in-house finite-differences code), and uses different discretization choices. As a second objective of the present contribution, we also intend to advance the statistical characterization and understanding of the BARC flow, by presenting for the first time a detailed discussion of the complete set of terms involved in the budget equations for the components 5D U∞ x y z Fig. 1 Sketch of the BARC geometry, with the reference system and the computational domain employed in the present work.
of the tensor of the Reynolds stresses. The availability of the complete budget for the Reynolds stresses is essential for improving LES and RANS closure models.

Numerical method and computational procedures
The BARC test case is a two-dimensional rectangular cylinder with streamwise length L and cross-stream dimension D, with aspect ratio A = L/D set at A = 5. Figure 1 sketches the geometry and the reference frame employed in the present work. A Cartesian coordinate system is placed at the center of the cylinder, with the x axis aligned with the flow direction, the y axis denoting the cross-stream direction and the z axis the spanwise direction. The incoming velocity is uniform upstream, aligned with the x axis and set at U∞. The Reynolds number Re = U∞D/ν is based on U∞, D and the kinematic viscosity ν of the fluid. As in the reference DNS work [10], the value of the Reynolds number considered here is Re = 3000. The incompressible Navier-Stokes equations are: where u, v and w denote the streamwise, cross-stream and spanwise components of the velocity, whereas p is the pressure. The mean field is indicated with capital letters, i.e. U = (U, V, W ), and the fluctuations with an apex, i.e. u = (u , v , w ). The boundary conditions of the problem are the no-slip and no-penetration conditions on the surface of the cylinder together with an unperturbed velocity (U∞, 0, 0) enforced at the far field, periodic conditions at the spanwise boundaries on account of spanwise homogeneity, and a convective outlet condition for the velocity, i.e. ∂u/∂t = −U∞∂u/∂x.
The DNS code, introduced by Luchini in [19], solves the incompressible Navier-Stokes equations on a staggered Cartesian grid. Second-order finite differences are used in every direction. The momentum equations are advanced in time by a fractional time-stepping method that employs a third-order Runge-Kutta scheme. The Poisson equation for the pressure is solved by an iterative SOR algorithm. The presence of the cylinder is dealt with via an implicit second-order accurate immersed-boundary method, implemented in staggered variables to be continuous with respect to boundary crossing and numerically stable at all distances from the boundary [18,19]. Hereinafter, all variables are in dimensionless form, with D as length scale, U∞ as velocity scale and D/U∞ as time scale.
The computational domain extends for −22.5 ≤ x ≤ 40 in the streamwise direction, −21 ≤ y ≤ 21 in the vertical direction and −2.5 ≤ z ≤ 2.5 in the spanwise direction, with the cylinder placed at −2.5 ≤ x ≤ 2.5 and −0.5 ≤ y ≤ 0.5. The computational domain is discretised with Nx = 1776, Ny = 942 and Nz = 150 points in the three directions, for a total of more than 250 millions grid points. The distribution of points is uniform in the spanwise direction, whereas geometrically varying grid spacing is employed in the streamwise and vertical directions to yield a higher resolution near the leading-and trailing-edge corners. There, the finest spacing occurs with ∆x = ∆y ≈ 0.0015. The ratio between neighboring streamwise cells is 1.005 over the body from the corners towards x = 0, and 1.02 before and after the body; the ratio between neighboring vertical cells is 1.01 for |y| ≤ 1.75 and 1.02 otherwise. This distributions leads to 892 points placed along the streamwise edge of the cylinder and 217 along the cross-stream edge.
When compared to the computational domain used in [10], i.e. −37.5 ≤ x ≤ 74.5 and −25 ≤ y ≤ 25, the present domain is slightly smaller, but discretized with a significantly finer mesh. Although known to be marginal for an accurate calculation of some high-order statistics [3], the spanwise dimension of the cylinder is identical to that of the reference study. In comparison, almost 15 times more points are placed over the body: about 8 times more points along the streamwise edge, and twice the number of points along the cross-stream edge.
To accumulate well-converged statistics, we have exploited both temporal averaging and ensemble averaging, by running several independent simulations. The first simulation is started from the solution corresponding to a steady potential flow, with its symmetry broken by injecting just after the leading-edge corners localized random noise of small amplitude, i.e. < U∞/10, during the first 5 time units. The simulation is then advanced for a considerably long time, approximately 400D/U∞, to allow for the flow to fully loose memory of the initial transition and reach a truly statistically-stationary state. Indeed, this flow takes a long time to develop as confirmed by [10], where a comparably long initial transient of 250D/U∞ was discarded too. Once the flow reaches the desired statistically steady state, the main simulation enters the production stage and accumulates statistics for further 341 time units. At the same time, from the latest stages of the preliminary run, 10 independent flow fields are saved and later used to produce independent initial conditions for 10 additional simulations. At the time of writing, they have accumulated useful statistics for additional 2004 time units, leading to a total averaging time of 2345D/U∞. The convergence of the statistics has been assessed by testing different sample sizes.
The creation of independent initial conditions is based on the idea of removing the large-scale coherent structures populating the flow, but without excessively distorting its spatial structure, so that the required initial transient can be kept to a minimum. The process is based on Fourier transforming the velocity values along each spanwise line, followed by a random phase change, and a final inverse transform. This procedure preserves the structure of the mean flow but effectively removes the large-scale structures, which are then regenerated quickly but independently on each sample. Although by visual observation the fields become fully independent after a very small simulation time, a conservative approach is adopted, and 20 further time units are discarded before starting the accumulation of statistics.
The temporal discretization uses a varying time step, to ensure that the Courant-Frederic-Levy CFL number remains at CF L ≤ 1; this condition produces an average value of the time step of ∆t ≈ 0.0013. The simulations have been run on the GALILEO supercomputer at CINECA. Each simulation uses 32 cores of a single computing node, and subdivides the spatial domain in smaller subdomains along the streamwise and spanwise directions. Each simulation needs approximately 14 hours to advance the flow by one time unit.

Instantaneous flow
We start by a qualitative characterization of a typical snapshot of the flow, which provides the opportunity of describing its main structures. They are visualised with the λ 2 criterion [14], based on the second largest eigenvalue λ 2 of the tensor are the symmetric and antisymmetric part of the velocity gradient tensor ∂u i /∂x j . The spatial orientation of the structures is additionally described by means of isosurfaces of the vorticity components ωx and ωz. Figure 2 plots isosurfaces for λ 2 , ωx and ωz. The sharp leading-edge corner determines the detachment point. In the initial part for −2.5 < x < −2 the separated shear layer remains two-dimensional and laminar, as demonstrated by comparing the contours of λ 2 with those of ωz, and then transitions to turbulence. The separatrix originating from the corner defines the spatial extent of a large recirculating region, a.k.a. the primary vortex; a secondary smaller counter-rotating vortex is formed around x = −1.5, where no significant three-dimensional structures are observed. As already described in previous works [8,35,39], first at x ≈ −2 a Kelvin-Helmoltz like instability of the shear layer appears, leading to a breakdown into large-scale spanwise tubes. Further downstream, the spanwise tubes stretched by the mean flow roll up and originate hairpin-like vortices. A sudden transition to turbulence occurs at x ≈ −1.3. Then, further downstream (x ≥ 0), the hairpin-like vortices are stretched and break down to elongated streamwise vortices, identified by the contours of ωx, confirming previous findings [8,35]. At these streamwise positions the flow is now fully turbulent, and the coexistence between small-and large-scales structures is clearly observed. The flow then proceeds towards the trailing edge, where again the sharp trailing-edge corner fixes separation in space, and a large turbulent wake ensues.

Temporal evolution of global and local quantities
A quantitative description of the flow begins with global quantities related to the mean flow field. Here and in the following, mean quantities, indicated by the operator · , are computed after time average and also by exploiting the homogeneity of the spanwise direction. The operators · z or · t will be used to indicate quantities averaged only along z, or in time respectively. For the BARC, the key quantities are the traditional dimensionless lift and drag coefficients, C and C d , i.e. the vertical and horizontal components of the aerodynamic force per unit width, normalised with 0.5ρU 2 ∞ D. Obviously for symmetry reasons C = 0 for an infinitely long simulation. Ref. [10] reports C d = 0.96 and does not mention the value of C ; our simulation provides a drag coefficient which is in agreement with theirs, at C d = 0.9437, and a pretty small residual value of C = −0.0155 which keeps decreasing with the integration time. Quantitative comparison of these and other quantities of interest with those of Ref. [10] are reported in Table 1, together with the range of values from other (experimental and numerical) studies collected and analysed in Ref. [5]. These studies are typically carried out at higher Re. Figure 3 shows the temporal evolution of the spanwise-averaged C z recorded in the primary (longest) simulation. Its excursions appear to be somewhat larger than the reference study, with the present signal reaching instantaneous values of up to +0.762 and −0.913, whereas in [10] they remain below 0.5 in absolute value. This is attributed to the much finer vertical grid spacing used in the present work, that allows to better capture the separation at the leading-edge corners; the same trend was observed also for the low-Reynolds laminar flow around a square cylinder in Ref. [37] (see table IV in their paper), and in Ref. [1] (see their figure  Table 1 Comparison between the present results and those from Ref. [10]. Additional columns report the data range (with their mean value in parentheses) for the experimental and numerical studies collected by Ref. [5].
The dominant time scales present in the flow are often [4,22,33,20] extracted by looking at localised peaks in the frequency spectrum S. The premultiplied periodogram of C z , shown in the top right panel and computed with data from the primary simulation only, shows a peak for the frequency f 1 = 0.1274. It is known from previous work [30,29,31,27,10] that this peak is associated to the vortex shedding in the wake. A more precise view of the dominant frequency can be obtained from the temporal autocovariance R C C (τ ) of the signal, shown in the bottom panel of figure 3. The first peak after the maximum at zero time separation is located at a time separation of τ = 1/f 1 . Overall, the identified frequency is quite similar to the time scale identified in [10], who measured f 1 ≈ 0.14. A frequency analysis of C d z (not shown) confirms the presence of a localised peak at a frequency of 2f 1 in the spectrum of the drag coefficient, induced by the alternate vortex shedding. Other studies, with different numerical approaches and at different Re, report for f 1 values of 0.098 or 0.105 [21], 0.1 [22], 0.11 [4], 0.132 [32]. Table 1 shows that often experiments and simulations tend to underestimate the shedding frequency. This is attributed to the sensitivity of the BARC flow to slight deviations from the ideal setting, or to simulation inaccuracies; the finite Re should not affect the computed value, because according to Ref. [29,36,28] at Re ≥ 3000 the Strouhal number has already reached its asymptotic value.
Several authors found that this flow is also characterised by a large-scale lowfrequency unsteadiness associated to a shrinkage and enlargement of the main recirculating region. Kiya and Sasaki in [15] and [16] report for this frequency a value of approximately 1/6 of the shedding frequency for a blunt flat plate at Re ≈ 10 4 . Similar results are also reported by [6] for the same flow and by [11] for a backward facing step. Ref. [10], instead, report f ≈ 0.042, that is approximately 1/3 of their shedding frequency. This low frequency is not detected clearly in the present data. However, the low Reynolds number and the amount of integration time preclude any firm statement. If such low frequency indeed exists at the present Re, it would require an extremely long integration to be reliably detected, as shown in [17] for the low-Re turbulent flow past a circular cylinder.
The As an example, we discuss below and show in figure 4 the frequency spectrum of the cross-stream velocity component v, although an equivalent picture is obtained when the streamwise velocity component is examined (not shown). The power spectrum at point a near the leading edge lacks a clearly dominating frequency, since the flow is essentially laminar. However, a small local peak is observed at the frequency f 1 that also emerges in the spectrum of C z , associated with the vortex shedding at the trailing edge. This suggests that at this Reynolds number the wake still influences the leading-edge region, as it happens at much lower Re [13]. The vortex-shedding frequency has been detected at these upstream sections also in Ref. [28] for Re = 13400, but not by Ref. [34] at Re = 40000; this suggests that the influence of the vortex shedding from the trailing edge on the first part of the shear layer gradually fades away as Re increases. Moving along the shear layer, a higher dominant frequency f ≈ 1.3 emerges, associated to the amplification of the velocity fluctuations in the shear layer [10] owing to the Kelvin-Helmotz instability [28]. The peak in the spectrum moves towards lower frequencies along the shear layer (f ≈ 1.44, 1.23, 1.21 for b, c, d respectively) and at the same time broadens. This frequency range is not entirely in agreement with Ref. [10], where the range is reportedly f ≈ 0.9 − 1.8. Similar results have been found also at larger Reynolds   Figure 5 shows the fields of the mean velocity (left) and pressure (right), obtained after spanwise and temporal average of the entire dataset (2345 time units). The mean flow separates at the sharp leading edge and reattaches downstream, before eventually separating again at the trailing-edge corner. Three recirculating regions are identified. The first is the large region identified by the shear layer separating at the leading edge. It starts from x s,1 = −2.5 and extends down to the reattachment point located at x e,1 = 1.455, with a length of L 1 ≡ x s,1 − x e,1 = 3.955. The centre of rotation of the primary vortex, defined as the stagnation point with U = V = 0 is found at (x c,1 , y c,1 ) = (−0.143, 0.83). In Ref. [10] it is noted that this region is associated with large negative values of pressure. Indeed, a pressure minimum is found close to the bubble centre, at (x p,1 , y p,1 ) = (−0.3, 0.99). Quantitative comparison of these and other quantities as measured here with those of the reference study [10], where L 1 = 3.65, are reported in Table 2. The larger size of the primary vortex may descend from the finer grid used here [25]. In terms of L 1 , literature values are 4.68 [4], 4.65 [21], 4.75 [22], 4.01 or 4.26 [32], 0.81 [22].

The mean flow
Within the large primary bubble, a second smaller counter-rotating recirculation bubble is observed. This secondary vortex is associated with the detachment of the reverse boundary layer caused by the adverse pressure gradient. Its characteristic length-scale is smaller: the secondary bubble extends between x s,2 = −1.87 and x e,2 = −0.91, with a shorter length of L 2 ≈ 0.96, with the centre of rotation placed at (x c,2 , y c,2 ) = (−1.3, 0.541). In Ref. [10] the length L 2 of this structure is L 2 = 1. Other literature values are 1.88 [4], 0.31 [21], 0.75 [7].
The third recirculating region is observed in the wake region, just after the trailing edge. This wake vortex extends from x s,3 = 2.5 up to x e,3 = 3.475, corresponding to a length of L 3 ≈ 0.975, and its centre of rotation is placed at (x c,3 , y c,3 ) = (2.915, 0.25). In Ref. [10] the length L 3 of this structure is L 3 = 1.2.
Other literature values are 0.76 [4], 1.4 or 0.7 [21], 0.94 [22], 0.9 or 0.81 [32], 1 [7]. Figure 6 plots the two components of the mean velocity vector at four streamwise stations, i.e. at x = −1.238 (corresponding to the small recirculation region), x = −0.1 (corresponding to the main recirculating region), x = 2 (after the reattachment of the primary bubble) and x = 2.85 (in correspondence of the recircula-  tion region in the wake). These velocity profiles are presented to provide a reference for experimental measurements, where hot-wire traverses could be employed, as for example in Refs. [20] and [33]. The evolution of the turbulent wake is described in figure 7, which assesses to what degree it obeys self-similarity. The left panel plots 1 − U (x, 0) where the abscissax = x − 2.5 measures the distance from the trailing edge. The centerline velocity defect follows a self-similar decay forx > 10, decaying proportionally to x −1/2 , and follows the curve where the two free parameters are fitted to A = 0.6 andx 0 = −0.5 (in Ref.    [10] of the normalised defect, with the cross-stream coordinate y 0.5 scaled with the characteristic wake thickness. The latter is defined such that U (x, y 0.5 ) = 0.5(1 + U (x, 0)). The profile is plotted at different streamwise stations downstream of x = 10. The profiles collapse reasonably well in the core of the wake, lending further support to the self-similarity of the wake after some distance from the body. The implied powerlaw spreading of y 0.5 ∼ x 1/2 is that of the plane turbulent wake. It should be noted, however, that the collapse is only marginal for |y/y 0.5 | > 1, indicating that at this distance from the body the self-similarity may not be complete. Moreover, the finite downstream domain length and the grid resolution in the far region might also be not enough to precisely capture the evolution of the turbulent wake. The mean flow field exerts its influence on the body, and determines the streamwise distribution of the coefficients c f (x) and cp (x), which express the longitudinal wall shear and the wall pressure made dimensionless with 0.5ρU 2 ∞ . These quantities are plotted in figure 8. The friction coefficient (left panel) has its largest negative value of −0.0166 at the leading edge, which is less than the value reported in Ref. [10], suggesting a different resolution of the interaction between the near-wall portion of the main shear layer and the secondary vortex. Friction then increases quickly to reach a plateau where it remains slightly positive for −1.87 < x < −0.91: this region is the trace at the wall of the small-scale recirculating bubble. Further downstream, c f becomes negative again, with a local minimum of value −0.0118 placed at x = 0.32, associated with the strong reverse flow observed in the large-scale recirculating region. When the trailing edge is approached, c f increases again to become positive for x ≥ x e,1 , i.e. downstream of the primary vortex. When compared to Ref. [10], the differences can be traced to the differences in the mean flow. The shift of the secondary vortex, for example, leads to a downstream shift of the first c f > 0 region, whereas the larger extension of the primary vortex results in a downstream shift of both the local minimum, which in Ref. [10] occurs at x ≈ −0.045 with a value of −0.013, and of the crossover point where c f is zero. The pressure coefficient cp (x) (right panel) is negative everywhere, in agreement with available information. A minimum of −1.205 is found in the region of the leading edge. After a sharp and localised increase, a mild decline starts at x ≥ −2.21 towards a local minimum of −0.855 at x = −0.470. This minimum is the footprint at the wall of the large negative pressure observed within the large-scale recirculating region. The pressure coefficient then increases again, and reaches its maximum of −0.124 at x ≈ 2.12, not long before the trailing edge separation. The computed cp (x) is in good general agreement with Ref. [10], although some differences are evident from figure 8. For example, for x < 0 our values are smaller (in absolute value), and the local minimum occurs later downstream. For x > 0, instead, the opposite is observed before the maximum, which takes place closer to the trailing edge and is less negative.

Single-point budget of the Reynolds stresses
The turbulent fluctuations in the flow are now described by studying the budget equation for the six independent components of the tensor for the Reynolds stresses. Terms for these budgets have been discussed e.g. in [23] for the channel flow, and have an obvious importance in turbulence modelling, but have not been fully documented yet for the BARC flow. As in the channel flow, two components, namely uw and vw , are zero because of statistical symmetry. Note that, in this Section, primes to indicate the fluctuating velocity components will be omitted for conciseness.
Very recently, Refs. [9,28,34] presented and discussed the production term in the equation for the turbulent kinetic energy, with the latter two studies focusing on the region close to the upstream corners. Hence, before presenting the Reynolds stresses, we start with figure 9, which shows the spatial map of the turbulent kinetic energy k = u i u i /2 (a repeated index imply summation), which is proportional to the tensor trace. The map is perfectly symmetrical on the two sides of the body, again supporting the adequacy of the statistical sample. Before x < −1.3, k is essentially zero, confirming the laminar state of the flow at the leading-edge separation and within the secondary vortex. For x > −1.3, however, k rapidly increases, signifying a quick transition to the turbulent state: the maximum is observed at (x, y) ≈ (0.2, 0.96). A further local maximum is observed in the wake at (x, y) ≈ (3.65, 0.36). However, this is only half the value of the maximum in the primary vortex. Overall, the map of k resembles the one reported in Ref. [8]: for example their global and local maxima are at (0.2, 0.9) and (3.7, 0.35). However, the intensity of k in Ref. [8] is larger than here: for example their maximum in the primary vortex is k ≈ 0.40 instead of the present maximum k ≈ 0.135.
We now proceed to describing the budget of the full Reynolds stress tensor. In tensorial notation, a compact form of the budget equation, stemming from the manipulation of the incompressible Navier-Stokes equations, and specialised to the present spanwise-homogeneous problem, can be written as: In Eq. (5), ψ x,ij and ψ y,ij are the fluxes in the x and y directions, defined as follows: whereas P ij , Π ij and ij denote the production, the pressure-strain and the pseudodissipation tensors, defined as: Figure 10 plots the four non-zero components of the Reynolds stress tensor in the x − y plane. Their general features have been already described in [10]. Once the flow becomes turbulent for x ≥ −2, the largest fluctuations are those of the streamwise component, in both the shear layer and the primary vortex, but the maxima of the three diagonal components are of comparable magnitude as also observed at larger Re [28]. The peak of uu is at (x, y) = (0.21, 0.97), i.e. above the centre of the primary vortex, and shows small values close to the cylinder's wall, confirming its dynamical association to the shear layer. The maximum of vv occurs slightly upstream than that of uu , i.e. at (x, y) = (0.05, 1.01), but it remains almost negligible near the wall, and is associated to the core of the primary vortex. Lastly, the position of the maximum of ww is almost identical to that of uu , at (x, y) = (0.24, 0.95). Interestingly, near the wall ww is the dominant component. This strikingly differs from e.g. plane Poiseuille flow, and is in agreement with the pattern seen in figure 2. In fact, large streamwise fluctuations are associated to the spanwise structures originated by the instability of the shear layer, whereas large spanwise and cross-stream fluctuations are linked to the streamwise-aligned structures populating the inner part of the large recirculating region. In the wake, uu is large in the shear layer separating from the trailing edge, with a peak at (x, y) = (3.46, 0.48), whereas smaller values are observed along the symmetry line y = 0. The other components, on the other hand, are large in core of the wake, with vv peaking at (x, y) = (3.59, 0). The ww component has quite a broad distribution, without a distinct peak. The general appearance is consistent with known results [10], but significant qualitative and quantitative differences emerge. Since the mean velocity field and in particular the extent of the primary vortex differ, different positions of the maxima are to be expected. Specifically, in [10] maxima at (0.07, 0.94), (0.46, 0.81) and (0.46, 0.78) are found for uu , vv and ww respectively. Hence the largest streamwise fluctuations occur before the other components, and the largest spanwise fluctuations occur closer to the cylinder side.
The off-diagonal component uv is also plotted in figure 10. The colour map is antisymmetric: here and in the following, in the text we refer to the top side of the Fig. 11 Contours of the production terms P 11 , P 22 , and P 12 .
cylinder. The most negative values are observed within the primary vortex, where the streamwise fluctuations correlate significantly to the vertical ones resulting in a minimum of uv located at (x, y) ≈ (0.19, 0.97). On the contrary, in the upstream portion for x < −1.5 below the separated shear layer, slightly positive uv are observed. Downstream the trailing edge uv is positive within the wake vortex, and negative elsewhere with a local minimum at (x, y) ≈ (3.5, 0.37). The shear stress uv plays a fundamental role in the production of the turbulent kinetic energy, which has been suggested [9] to become negative in a localised region in the BARC flow. This point will be discussed later. Figure 11 shows P 11 , P 22 and P 12 , i.e. the production terms for uu , vv and uv . We recall that P 33 = 0 because ∂W/∂x = ∂W/∂y = 0, and that the sum of the diagonal terms is half the production of turbulent kinetic energy. In the shear layer for x ≥ −2 P 11 is positive and P 22 is negative (albeit smaller). This implies that, once the shear layer instability takes place, energy is drained from the mean flow to feed streamwise fluctuations, whereas the opposite occurs for the cross-stream component v. Further downstream, in the core of the primary vortex P 11 and P 22 are both positive, and the mean flow feeds both uu and vv . Furthermore, close to the wall in the downstream portion of the body, i.e. for x ≥ 0, P 11 becomes mildly negative, hence a sink for uu . In this region flow reversal takes place, where ∂U/∂y < 0, and the streamwise fluctuations are weaker. Similar considerations can be put forward for the wake. Indeed, P 11 is mildly positive in the recirculating region behind the trailing edge, but becomes negative more downstream along the centerline y = 0, where indeed low values of uu take place. On the contrary, along the centerline line and outside the recirculation, P 22 is positive with relatively large values, denoting production of v fluctuations as shown also in figure 10. P 12 is shown in the bottom panel of figure 11. It is worth noting [12] that, since uv is not positive-definite, interpreting P 12 , Π 12 and 12 in terms of production or dissipation requires to account for the sign of uv . Negative values of P 12 are observed everywhere -except in a flat region very close to the cylinder side where it is slightly positive -with a global minimum within the primary vortex at (−0.59, 1.03) and a local minimum in correspondence of the shear layer separating from the trailing edge at (2.65, 0.5). Overall, since in the region with large negative P 12 uv is negative too, a production of −uv takes place. With analogous reasoning, the negative P 12 in the first part of the shear layer separating from the leading edge corners indicates a sink for the positive uv .
Let us now focus on the energy production, by examining the contributions to P 11 and P 22 . Each production term can be split in two, to isolate contributions related to the streamwise variation of the mean flow from those related to its cross-stream variation:  Figure 12 separately plots each of these four terms. P 11,a and P 22,b have opposite sign and similar appearance in the domain, owing to the incompressibility constraint ∂U/∂x = −∂V ∂y. The figure shows that the main contribution to P 11 changes depending on the region considered. On the shear layer, for x ≤ −1.5 positive production comes from P 11,a and is due to the negative values of ∂U/∂x Fig. 13 Left: contour of the production term for the turbulent kinetic energy P k = (P 11 + P 22 )/2; the black line denotes the 0 value. Right: contribution of to the production term for the turbulent kinetic energy from the interaction of uv and the mean shear ∂U/∂y + ∂V /∂x, i.e. (P 11,b + P 22,a )/2. associated with the shear layer. Here P 11,b is negative, since uv > 0 and ∂U/∂y is positive everywhere at this station. Further downstream, both P 11,a and P 11,b become positive, since uv < 0 (see figure 10). Production of uu in the core of the main recirculating region is instead given by P 11,b , and is determined by the positive values of ∂U/∂y: here ∂U/∂x < 0, leading to negative P 11,a . Near the body P 11 is dominated by P 11,a , as P 11,b → 0 faster for y → 0.
Unlike P 11 , the main contribution to P 22 comes from P 22,b , as the companion component P 22,a is significantly smaller similarly to what found experimentally by Ref. [28] in the first part of the shear layer at Re = 13400. Therefore, the main driver here is ∂V /∂y. The negative production in the shear layer and in the reverse boundary layer for x ≤ 0 is due to ∂V /∂y > 0, whereas the positive values of P 22 in the core of the main recirculating region to ∂V /∂y < 0.
Cimarelli et al. in [9] discuss the presence of a thin region with significant negative production rate of turbulent kinetic energy, localised in the shear layer close to the leading-edge corners (see figure 2 of their paper). Such a negative production region is not observed here. In fact, as shown in the left panel of figure  13 where the total production P k = (P 11 +P 22 )/2 is plotted, mildly negative values are indeed found, but they are close to the leading-edge corners below the shear layer, close to the cylinder wall for x > 0 and behind the trailing edge in the wake region. Overall, P 11 contributes to P k more than P 22 almost everywhere in the domain, as seen in figure 11; this is true also for x = −2.5, unlike what found at larger Re in Ref. [34], where a larger contribution from P 22,b is reported. In [9] negative production is linked to the interaction between uv and the mean streamwise and vertical shears ∂U/∂y and ∂V /∂x. This is confirmed in the right panel of figure 13 where the sum (P 11,b +P 22,a )/2 is plotted: a region with negative (P 11,b +P 22,a )/2 < 0 is indeed observed in the first portion of the shear layer (more precisely, this region is due to P 11,b alone, since P 22,a is negligible). However, in this region of the domain P 11,a is dominant, eventually resulting in a positive production term. Figure 14 shows the pressure-strain term, to elucidate how the pressure-mediated redistribution of energy affects the various components of the tensor of turbulent stresses. Π 11 is negative almost everywhere, implying that uu is redistributed to the other components. This takes place mainly in the shear layer and in the core of the primary vortex, where both Π 22 and Π 33 are positive, signaling that energy is received from uu by both vv and ww . Along the shear layer Π 22 > Π 33 , so that uu preferentially provides energy to vv ; the opposite occurs in the core of the recirculating region, where ww is the largest receiver. Following [12], the bottom left panel visualises this concept, by plotting in red (or blue) the areas where energy is preferentially transferred from uu to ww (or vv ). In fact, the incompressibility constraint mandates that: Hence, this panel plots the quantity Π 33 /|Π 11 | under the condition that Π 11 < 0, Π 22 > 0 and Π 33 > 0, and the colour scale is chosen to draw in red where the transfer from uu towards ww prevails over that towards vv , and in blue the opposite case. Near the cylinder wall, predominantly in the downstream part, Π 22 is negative, whereas both Π 11 and Π 33 are positive. This essentially visualises the same splatting phenomenon observed for example by [23] in the plane turbulent channel flow: owing to the non-penetration boundary condition v is redistributed among u and w. Interestingly, near the BARC wall Π 33 > Π 11 , so that here the redistribution mainly occurs towards the spanwise velocity component. This is visualised in the bottom right panel of figure 14 which plots the quantity Π 33 /|Π 22 | under the condition that Π 22 < 0, Π 11 > 0 and Π 33 > 0; the colour scale is chosen to highlight in red (blue) the transfer of vv towards ww ( uu ). This is in agreement with the observed prevalence of turbulent structures aligned with the x direction and with larger intensity of ww compared to uu and vv in the near-wall region (recall also that P 11 < 0 is a sink for uu in this region of the flow).
The middle right panel of figure 14 shows Π 12 , the pressure-strain term for uv . Π 12 is positive almost everywhere -except for the first portion of the shear layer separating from the leading edge -with the largest values at the cylinder wall for x > 0 and local maxima at (x, y) ≈ (−0.5, 1.03) within the primary vortex and at (x, y) ≈ (2.7, 0.53) in the shear layer separating from the trailing edge. Therefore, Π 12 is a sink for uv everywhere. Indeed, it dissipates positive uv in the shear layer for x < −1 and negative uv elsewhere with a peak of activity at the cylinder side.
The last term at the r.h.s. of equation (5) is the dissipation tensor. The relevant components are drawn in figure 15. For the diagonal components, large values of the dissipation occur in the shear layer for x ≥ −1.5 and in the core of the primary vortex. Interestingly, the largest values are observed for uu and ww . Moreover, for both the streamwise and spanwise components, unlike for vv , large dissipation is also seen close to the cylinder, where viscous effects are dominant. This qualitatively resembles the observations (see e.g. [23]) put forward for the channel flow, and confirms the larger degree of universality for the dissipative phenomena near a wall. The dissipation term of uv , 12 , is shown in the bottom right panel of figure 15. 12 is negative almost everywhere, except in the first portion of the shear layer separating form the leading edge where it is slightly positive. The global minimum occurs at the trailing-edge corner where the viscous phenomena are dominant, whereas a local minimum, of almost a lower order of magnitude, is seen in the primary vortex at (x, y) ≈ (−0.1, 1); similarly to 22 , 12 close to the cylinder side 12 is low. Overall, like for channel flow 12 is of a lower order of magnitude compared to P 12 and Π 12 almost everywhere, and therefore its contribution to the production/dissipation of uv is negligible.
The production, pressure-strain and dissipation tensors can be put together as a single source term at the r.h.s. of equation (5). Figure 16 plots contour of the source term ξ ij = P ij + Π ij − ij , together with the field lines of the flux vector with components (ψ x,ij , ψ y,ij ). These lines show how the excess of u i u j gets redistributed in the flow. For uu , a positive ξ 11 shows an excess of production of uu which is redistributed along the field lines originating from a singularity point placed in the shear layer at (x, y) ≈ (−1.7, 0.82), close to the maximum of ξ 11 (corresponding to the maximum of P 11 ). Following the field lines, some of them are observed to carry uu downstream towards the wake, whereas others enter the large recirculating region of the primary vortex. Much like the mean flow, these lines show a spiraling motion and vanish because of dissipation close to the rotation centre of the vortex, at (x, y) ≈ (−0.4, 0.8). Here, figure 15 confirms the importance of viscous effects. The field lines of vv are different. They originate in the shear layer at (x, y) ≈ (−2.12, 0.82), close to the maximum of ξ 22 (where, as seen in figures 11 and 14, Π 22 > |P 22 |). Some of them carry part of the excess of vv downstream towards the wake region, whereas others are attracted by the solid wall. In fact, near the wall ξ 22 < 0 because of the splatting effect, so that the wall acts as a sink. The field lines of ww qualitatively resemble those for uu . Indeed, ξ 33 is maximum in the shear layer, where the pressure strain produces spanwise fluctuations. As for uu , some lines advect part of the excess of ww downstream towards the wake, whereas others enter and remain in the recirculating bubble, spiralling inwards towards the centre of rotation of the primary vortex, and ending at (x, y) ≈ (−0.44, 0.83) due to viscous dissipation, i.e. slightly upstream the vanishing point of the lines of uu . Interestingly, these lines reach the most upstream part of the primary vortex, and highlight the zone where the reverse flow separates because of the adverse pressure gradient. The field lines of uv have a distinct shape. Some originate at the solid wall where ξ 12 is maximum, whereas others originate in the wake and carry uv upstream. The line set is attracted by the large sink region with ξ 12 < 0 placed in the shear layer, and approach the leading edge corner. Again, when discussing these field lines, fluxes of uv should not be interpreted in terms of energy transfer, as uv is not a positive-definite quantity [12].

Conclusions
We have studied with a Direct Numerical Simulation the BARC benchmark, i.e. a 5:1 rectangular cylinder immersed in a uniform flow, in the turbulent regime. Despite the simple and somewhat idealised shape of the cylinder, the BARC flow contains complex and fascinating features typical of separating and reattaching flows over complex geometries. The large scales associated with the flow instabilities coexist and interact with the smaller scales associated with turbulent motions to create a rich and intricate scenario. A fully reliable statistical characterisation of the BARC flow is still lacking: large scatter of data is observed already for first-order statistics [5], as the flow is very sensitive to geometry details and various types of external disturbances which are unavoidable in experiments, and to various modelling and discretisation choices in numerical simulations. Similarly, the agreement between experimental and numerical information is not entirely satisfying yet.
This study has replicated the first DNS of the BARC flow, recently carried out by Cimarelli et al. [10] at a value of the Reynolds number (based on the body thickness and the uniform incoming velocity) of Re = 3000. The numerical toolbox employed here is quite different though, as we have used a in-house finitedifferences solver instead of the finite-volumes OpenFOAM package. Moreover, we have used a finer spatial discretisation -the total number of point is more than one order of magnitude larger -and a longer averaging time. The goal of the study is to contribute to a solid and reliable DNS-based benchmark, that could be used as reference for RANS-or LES-based simulations of the BARC flow.
We have described and discussed the main features of the flow in terms of firstand second-order statistics, and the main differences with the available DNS results [10] have been reported. Our results present a near-perfect symmetry along the y = 0 centerline, which demonstrates the adequacy of the statistical sample. The longitudinal extent of the largest primary vortex has been found to be 10% larger, and the secondary vortex 10% smaller than the reference study. The position of the recirculating regions are different too. Also, we do not confirm the presence of a spot in the leading-edge shear layer where the production rate of turbulent kinetic energy becomes negative [9].
Moreover, the statistical understanding of the BARC flow has been furthered, and for the first time the whole set of terms involved in the budget equation for the tensor of the Reynolds stresses has been presented and discussed in detail. The analysis highlights the strongly anisotropic and inhomogeneous nature of the flow. The shear layer separating from the leading-edge corner is initially laminar, but its instability soon leads to the fluctuating field draining energy from the mean flow and feeding u u , which is the dominant contributor to the turbulent kinetic energy. The other components, v v and w w , are instead produced by a redistribution of u u driven by the pressure-strain term, and have a lower intensity. The cross-stream component of the mean flow in the shear layer is fed from the fluctuating field. The excess of turbulent kinetic energy is partially advected towards the wake, and partially transported within the large primary vortex. Here, u u and v v are fed by energy drained from the mean flow, and redistribution moves energy away from u u towards v v and w w . In that portion of the primary vortex where reverse flow takes place (close to the cylinder wall), a splatting effect is observed, with v v being redistributed to u u and (mainly) w w . In this region the production term for u u is negative, indicating that energy is drained from u u to feed the mean flow. Hence, in the reverse boundary layer the turbulent kinetic energy is mainly organised in spanwise fluctuations. Energy dissipation takes place mainly in the core of the primary vortex, where it is comparable for all the normal stresses, and in the vicinity of the cylinder wall, where the largest dissipation is observed for the wall-parallel components.
While the present state of affairs calls for further reduction of the uncertainty in the statistical description of this flow, the present dataset -which is made available to the research community at the following DOI: https://doi.org/10.5281/zenodo.4472682 -provides an useful addition to the existing knowledge on the BARC benchmark, an interesting flow for turbulence modelling.