A priori tests of eddy viscosity models in square duct flow

We carry out a priori tests of linear and nonlinear eddy viscosity models using direct numerical simulation (DNS) data of square duct flow up to friction Reynolds number Reτ=1055\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {Re}}_\tau =1055$$\end{document}. We focus on the ability of eddy viscosity models to reproduce the anisotropic Reynolds stress tensor components aij\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{ij}$$\end{document} responsible for turbulent secondary flows, namely the normal stress a22\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{22}$$\end{document} and the secondary shear stress a23\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{23}$$\end{document}. A priori tests on constitutive relations for aij\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{ij}$$\end{document} are performed using the tensor polynomial expansion of Pope (J Fluid Mech 72:331–340, 1975), whereby one tensor base corresponds to the linear eddy viscosity hypothesis and five bases return exact representation of aij\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{ij}$$\end{document}. We show that the bases subset has an important effect on the accuracy of the stresses and the best results are obtained when using tensor bases which contain both the strain rate and the rotation rate. Models performance is quantified using the mean correlation coefficient with respect to DNS data C~ij\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{C}}_{ij}$$\end{document}, which shows that the linear eddy viscosity hypothesis always returns very accurate values of the primary shear stress a12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{12}$$\end{document} (C~12>0.99\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{C}}_{12}>0.99$$\end{document}), whereas two bases are sufficient to achieve good accuracy of the normal stress and secondary shear stress (C~22=0.911\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{C}}_{22}=0.911$$\end{document}, C~23=0.743\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{C}}_{23}=0.743$$\end{document}). Unfortunately, RANS models rely on additional assumptions and a priori analysis carried out on popular models, including k–ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} and v2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^2$$\end{document}–f, reveals that none of them achieves ideal accuracy. The only model based on Pope’s expansion which approaches ideal performance is the quadratic correction of Spalart (Int J Heat Fluid Flow 21:252–263, 2000), which has similar accuracy to models using four or more tensor bases. Nevertheless, the best results are obtained when using the linear correction to the v2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^2$$\end{document}–f model developed by Pecnik and Iaccarino (AIAA Paper 2008-3852, 2008), although this is not built on the canonical tensor polynomial as the other models.

simplest form, eddy viscosity or explicit algebraic stress models (EASMs) rely on the assumption that the anisotropic stress tensor a i j ≡ τ i j − 2/3kδ i j (where τ i j ≡ u i u j is the Reynolds stress tensor and k ≡ 1/2τ ii the turbulent kinetic energy) is linearly proportional to the local mean rate-of-strain: where ν t is the eddy viscosity, S i j is the mean rate-of-strain: and u i are the components of the velocity vector. The overline symbol is used to indicate ensemble averages, the prime symbol indicates turbulent fluctuations and x i are the spatial coordinates in the three directions i = 1, 2, 3. Variables with tilde indicate modelled quantities in which the only modelling assumption is the constitutive relation (e.g. the anisotropic Reynolds stress tensor a i j , modelled with the linear eddy viscosity hypothesis in Eq. (1)). Despite their popularity, linear eddy viscosity models are inaccurate in predicting flows with large anisotropies and inequalities of the normal stresses [23]. A notable case in which linear eddy viscosity models fail to reproduce the correct flow physics is represented by turbulent secondary flows, as occurring in ducts of non-circular cross section [24], with the simplest prototype represented by the turbulent square duct flow [21]. The flow mechanism generating these secondary flows can be traced back to nonzero mean streamwise vorticity ω x [7,36], which is related to the cross-stream velocity components through the stream function ψ: The analysis of the mean streamwise vorticity equation, shows that the only sources of ω x are the derivatives of the secondary shear stress v w and the anisotropy of the normal stresses v and w and they both have the same order of magnitude [21]. The linear eddy viscosity hypothesis τ i j = 2/3kδ i j − 2ν t S i j inevitably returns τ 23 = 0 and τ 22 = τ 33 = 2/3k, leaving a homogeneous equation for ω x , with homogeneous boundary conditions; therefore, the only possible solution is ω x = 0. Hence, the linear eddy viscosity (1) cannot predict secondary flows [37]. These limitations of the linear eddy viscosity hypothesis led [22] to derive a more general form for the constitutive relation (1), in which the anisotropic stress tensor of a three-dimensional flow is exactly represented as a nonlinear polynomial expansion of ten independent tensor bases and ten coefficients. Moreover, Pope [22] provided analytical expressions for the tensor polynomial coefficients in the case of two-dimensional mean flow (two mean velocity components), for which exact representation of a i j only requires three tensor bases. Using the same set of integrity bases, Gatski and Speziale [8] found the analytical expressions of the ten tensor polynomial coefficients for a general three-dimensional flow. Different authors [9,12] showed that, except for some degenerate cases, five tensor bases are sufficient to represent the anisotropic stress tensor, as a i j is traceless and symmetric; therefore, it has only five independent components. Additional bases are only necessary if S i j has multiple eigenvalues (i.e. axisymmetric strain) or if the vorticity vector is aligned with one of the principal directions of strain [12]. After these seminal studies, many nonlinear RANS models have been proposed, the most common being quadratic [17,19,27,31,32,37] and cubic models [2]. Nonlinear eddy viscosity models have been successfully used to reproduce turbulent secondary flows in different geometries such as non-circular ducts [13,32,33], wing-body junction [3,28,41], turbine-hub flow [38], and to recreate the spanwise rollers in Couette flow [34].
Nevertheless, only few studies used DNS data to assess the accuracy of eddy viscosity models for predicting turbulent secondary flows. Mompean et al. [16] carried out a priori tests of nonlinear k-ε models using direct numerical simulation (DNS) data at friction Reynolds number, Re τ = u τ h/ν ≈ 150 (where u τ = √ τ w /ρ is the mean friction velocity, τ w the mean wall shear stress, h the duct half side, ν and ρ the kinematic viscosity and density of the fluid, respectively). Mompean [15] also carried out additional tests on the nonlinear model proposed by Speziale [37] showing that it underpredicts the intensity of the secondary flows. Rice et al. [25] carried out particle image velocimetry of the flow along supersonic streamwise corners and reported good agreement with RANS simulations using the quadratic correction of Spalart [32].
In this study, we aim at assessing the accuracy of popular RANS models using the square duct flow as benchmark. The accuracy of RANS models can be quantified using two tools, namely a priori and a posteriori tests, both with advantages and disadvantages. In a priori tests, modelled quantities such as eddy viscosity, Reynolds stresses and model coefficients are compared to reference high-fidelity data, without advancing the RANS equations or additional transport equations. The main drawback of this approach is that it only provides information on the accuracy of the modelled Reynolds stresses and not on the mean flow, which is the final goal of RANS. An advantage of this approach is that it gives us information on the model accuracy in the best-case scenario. Hence, if a priori tests show inaccurate prediction of the Reynolds stress tensor, solving additional transport equations can only worsen the results. On the other hand, a posteriori tests have the advantage of providing mean flow data, but solving the equations we introduce the effect of the numerical scheme. Hence, to carry out fully consistent a posteriori tests one should use the same solver and the same mesh to generate reference high-fidelity data and RANS results. In this study, we limit ourself to a priori tests; therefore, we are unable to inspect the mean flow velocity but we can only assess the accuracy of modelling assumptions on the Reynolds stress tensor.
We carry out a priori tests using recent DNS data of square duct flow at friction Reynolds number Re τ = 150-1055 [14,21], thus extending the Reynolds number range of previous studies. We show that results become independent from Reynolds number as Re τ increases and therefore we mainly focus on the highest Reynolds number data. We carry out two types of a priori tests, the first on nonlinear constitutive relations for the anisotropic stress tensor a i j to determine the maximum accuracy that one should expect for a given number of tensor bases, and the second on popular nonlinear RANS models to determine the validity of their assumptions and whether they approach the maximum accuracy.
where S is the mean rate-of-strain tensor (2) and the mean rate-of-rotation tensor: A i j indicates the matrix associated with the tensor A and the notation AB = A ik B k j is used for tensor multiplication, whereas {AB} = A ik B ki indicates the trace operation. An exact representation of the anisotropic stress tensor ( a = a) can be recovered for N = 5, apart from the degenerate cases of axisymmetric flow and vorticity vector aligned with an eigenvector of S, for which additional bases are required. Truncated versions of the tensor polynomial (5) constitute the theoretical framework for the development of nonlinear eddy viscosity models [2,32], machine learning algorithms for the optimization of EASMs [10,39] and uncertainty quantification of EASMs [4,40].
In order to carry out a priori tests on the constitutive relation (5), we compute the coefficients G (n) , using the anisotropic stress tensor from DNS data of square duct flow [21]. G (n) can be obtained by multiplying equation (5) by the basis T (m) and taking the trace operation [9], thus obtaining a linear system of N equations for the unknown coefficients G (n) : For N = 1, this approach leads to the classical definition of linear eddy viscosity: In the more general case, we solve Eq. (8) to find the coefficients G (n) for different 1 ≤ N ≤ 5 and then evaluate the modelled anisotropic stress tensor a from (5). By solving the algebraic system (8), we let the coefficients G (n) vary in space to approximate the "exact" anisotropic stress tensor from DNS a. Hence, the coefficients G (n) are the ones that bring the most accurate representation of the anisotropic stress tensor (5), for a given number of bases. This kind of analysis was proposed by Jongen and Gatski [9], who used experimental data to assess the validity of the eddy viscosity hypothesis. Schmitt [29,30] also carried out the same analysis for linear and nonlinear constitutive relations highlighting the limitations of the eddy viscosity hypothesis. Nevertheless, this approach has never been applied systematically for different truncations of Pope's polynomial expansion. In the following analysis, only three components of the Reynolds stress tensor are reported, namely a 12 , a 22 and a 23 , as in square duct flow a 12 (y, z) = −a 13 (z, y), a 22 (y, z) = a 33 (z, y) (where y and z are the cross-stream coordinates) and a 11 does not appear in Eq. (3) (i.e. it does not influence the secondary flows). In the following, we use the term primary shear stress to denote a 12 , as it appears in the mean streamwise momentum equation and it directly influences the skin friction coefficient, whereas we use secondary shear stress for a 23 , as it appears in the cross-stream mean momentum equation.
Moreover, we only report in detail the anisotropic stresses for selected N , corresponding to popular choices for the development of eddy viscosity models. In particular, N = 1 corresponds to the linear eddy viscosity hypothesis, N = 2 the nonlinear correction proposed by Spalart [32], N = 4 has been used to develop quadratic eddy viscosity models, and for N = 5 we recover the exact representation a = a. Note that the exact representation can only be recovered when the ideal coefficients obtained from (8) are used; therefore, using five bases alone, with arbitrary coefficients, is not sufficient to recover a = a. Figure 1 shows the scatter plots of the modelled anisotropic stress tensor components, a 12 , a 22 , a 23 as a function of the ones from DNS a 12 , a 22 , a 23 at Re τ = 1055, normalized with the friction velocity from DNS (a + = a/u τ ). Figure 1a-c shows that the linear eddy viscosity hypothesis (N = 1) returns accurate representation of a 12 , the stress component directly affecting the streamwise velocity. The success of the linear eddy viscosity hypothesis in predicting the primary shear stress can be traced back to the strong similarity between square duct and canonical wall-bounded flows, namely to the fact that this flow can be seen as superposition of four concurrent walls [21]. If we assume that S 12 and S 13 are the most relevant components of S i j , the eddy viscosity (9) becomes which reduces to −u v /(2∂u/∂ y) of canonical wall-bounded flows, away from the side wall. On the other hand, a 22 ≈ 0 and a 23 are uncorrelated with respect to reference data, in agreement with the fact that the linear eddy viscosity hypothesis cannot predict the occurrence of secondary flows (Fig. 1b, c).
Figure 1d-f shows that using N = 2 improves the prediction of both the normal stress and secondary shear stress. Nevertheless, a 22 is consistently overpredicted, as well as the most intense values of a 23 . For N = 4, exact representation of a 22 is recovered, together with a very good prediction of a 23 (Fig. 1g-i). As predicted from theory, five bases allow us to recover the stresses from DNS ( a ≡ a, Fig. 1j-l).
The scatter plots allow us to easily compare the modelled stresses to DNS data, but they do not provide information on the spatial distribution of the stress components and therefore it is not possible to understand how the error is spatially distributed. For this reason, we also report a 12 , a 22 , a 23 at Re τ = 1055 in a duct quadrant, for N = 1, 2, 4 ( Fig. 2a- 12 , a 22 and a 23 , (Fig. 2j-l). The spatial distribution of the primary shear stress a 12 is well predicted by the linear model, and only minor differences are observed between the linear and nonlinear constitutive relations. As observed from the scatter plots, the linear eddy viscosity hypothesis returns a 22 ≈ 0 (Fig. 2b). The spatial distribution of a 23 (Fig. 2c) is similar to DNS only in the core region of the duct, where a large patch of negative stress approximately matches DNS values, whereas an incorrect stress sign is predicted towards the walls, where the flow is more anisotropic. Using N = 2 leads to slightly less accurate a 12 with a small tail close to the side wall, which is absent in the linear case and in DNS (Fig. 2d).
A significant improvement in the prediction of both a 22 and a 23 is observed for N = 2. We find that the largest error in a 22 is localized towards the side wall, where the sharp variation of a 22 is not captured (Fig. 2e). The spatial organization of the secondary shear stress also improves compared to the linear model, but the values of a 23 are larger than DNS (Fig. 2f). Similarly to the linear case, the sign and intensity of a 23 are correct Table 1 Correlation coefficient C i j (11) between anisotropic stress tensor of DNS a i j and modelled anisotropic stress tensor a i j (Eq. (5)) at friction Reynolds number Re τ = 1055, for different numbers of tensor bases N in the duct core, where the flow is more isotropic, whereas the error with respect to DNS becomes larger close to the corner and the walls, where very high stress values are predicted. Using four tensor bases allows us to obtain a spatial distribution of a 22 which is undistinguishable from DNS data, whereas differences are still observed for a 23 , although the main stress topology is captured (Fig. 2h, i). Also in this case the greatest error on a 23 is localized in proximity of the corner and at the walls, where the stress is overpredicted. In order to quantify the error with respect to DNS, we also define the averaged correlation coefficient between a i j and a i j : where the angle brackets denote average over the duct cross section: By construction, C ∈ [−1, 1], where C = 1 indicates the perfect correlation C = −1 the negative correlation.
In Table 1 and Fig. 3a, we report the mean correlation coefficient C i j at Re τ = 1055 for increasing number of tensor bases. As first gauged from qualitative inspection of the scatter plots and mean flow fields, we find that correlation with respect to DNS increases with N . The primary shear stress has a correlation coefficient C 12 ≈ 1 already for N = 1; therefore, the linear eddy viscosity hypothesis can accurately predict a 12 , whereas a 22 and a 23 are essentially uncorrelated with respect to DNS. For N = 2, we find a good correlation for the normal stress component C 22 ≈ 0.9 and relatively good also for C 23 ≈ 0.74. On the other hand, passing from N = 2 to N = 4 leads to perfect correlation on the normal stress C 22 = 1, but it does not increase the accuracy on a 23 . Figure 3b, c shows C i j as a function of the friction Reynolds number Re τ = [150, 227, 519, 1055] for fixed number of bases N = 2 and N = 4, respectively. We observe a minor Reynolds number effect, which tends to saturate for increasing Reynolds number, suggesting that Reynolds number independence is achieved for Re τ 500.
The analysis has been carried out using sets of tensor bases ordered as they appear in the polynomial expansion (5). This is motivated by the fact that nonlinear models are typically built following this approach, but for a given number of tensor bases different combinations are possible and results can be affected by the choice of the subset. For this reason, we explore also this possibility and in Table 2 we report the mean correlation coefficient C i j for all possible subsets with N = 2, 3, 4. The subsets are limited to the first five bases, as only using these ones it is possible recover exact representation of the anisotropic stress tensor (i.e. using T (1) , T (2) , T (3) , T (4) , T (5) , a = a, whereas using T (1) , T (2) , T (3) , T (4) , T (6) , a = a). Table 2 shows that a 12 is barely affected by the choice of the subset, whereas the accuracy on a 22 and a 23 can have a large variation. The best subset for a given number of bases is highlighted in boldface, showing that T (2) and T (5) are essential to describe both a 22 and a 23 , whereas bases T (3) and T (4) alone do not bring any improvement with respect to the linear case. The reason for this is that T (2) and T (5) have the same structure and contain both S and , contributing more to the flow anisotropy.
Moreover, in Fig. 4 we report the coefficient G (1) in a duct quadrant, for different truncations of Pope's polynomial expansion. Ideally, for a smooth flow, the tensor bases are smooth functions and so are the coefficients. Nevertheless, increasing N can give rise to sharp spatial variations of T (n) . This causes system (8) to be ill-conditioned and difficult to solve numerically, a problem which we overcome using a truncated singular value decomposition to invert the system matrix. The figure shows that for N > 2, G (1) exhibits sharp spatial variations, and a similar result holds for the other coefficients (not shown). More aggressive matrix  Table 2 Correlation coefficient C i j (11) between anisotropic stress tensor of DNS a i j and modelled anisotropic stress tensor a i j (Eq. (5)) at friction Reynolds number Re τ = 1055, for different combinations of tensor bases. The optimal subset for each N is highlighted in boldface  regularization techniques are possible and lead to smoother coefficients, but in this case we do not recover the exact representation of a i j for N = 5. We recall that even though coefficients G (n) are spiky, a i j is always smooth (Fig. 2). The analysis reveals that, even for a smooth flow solution as incompressible duct flow, the exact coefficients can exhibit unphysical spatial variations and therefore modelling, or even using, the exact coefficients G (n) for N > 2 might be impossible or at least very difficult in practical calculations.

A priori tests of RANS models
The idea of generalized eddy viscosity hypothesis of Pope [22] has inspired several nonlinear eddy viscosity models using truncated forms of the tensor polynomial (5). Nevertheless, eddy viscosity models inevitably rely on additional assumptions to represent the coefficients G (n) needed by the polynomial expansion and therefore they can only approach the accuracy of the nonlinear constitutive relations using the same number of tensor bases. In reality, it is not possible to tune the coefficients G (n) to achieve the maximum accuracy for all variety of flows, and therefore some models try to account for this by using more than five bases. We carry out a priori analysis on different turbulence models and compare their accuracy to the one of the constitutive relations using exact coefficients G (n) . First, we focus on the popular k-ε turbulence model [11] and then on nonlinear [32] and linear [20] corrections specifically developed to improve the prediction of turbulent secondary flows. In the following, we denote the anisotropic stress tensor of the RANS models with a hat symbolâ i j to distinguish it from a i j , obtained using the exact coefficients G (n) . We also introduce the correlation coefficient between a i j and DNS data a i j : where the angle brackets denote averaging over the duct cross section (12). For each RANS model, we report the correlation coefficient with respect to DNS C i j , Table 3. The value of C i j allows us to assess how the model performs with respect to DNS data, but this is an unfair comparison as RANS models are not supposed to match DNS data and can at most reach the ideal accuracy of the constitutive relation they are built on. Hence, in Fig. 5 we report C i j as a function of C i j , with the same number of tensor bases (e.g. C i j of the linear k-ε model as a function of C i j for N = 1); therefore, the axes bisector represents the maximum mean accuracy.
The k-ε turbulence model [11] is a two-equation model, requiring transport equations for the turbulent kinetic energy k and for the dissipation rate ε = 2ν S i j S i j (where S i j is the fluctuating rate-of-strain tensor). The class of k-ε models became very popular over the years, and it is among the most common choices in many commercial fluid dynamics solvers [23]. Besides the transport equations for k and ε, the model requires the specification of the eddy viscosity, which is defined using the assumption: where c μ is a constant and f d a damping function to avoid high values of eddy viscosity close to the wall. With this assumption, a generalized constitutive relation for the k-ε model can be derived based on the tensor polynomial expansion (5) [2]: where c 1 -c 7 are the model coefficients (Table 4). Equation (15) uses eight tensor bases to approximate a, whereas in Sect. 2 we showed that five bases allow us to recover the exact anisotropic stress tensor. The reason for using N > 5 is that the exact representation can only be recovered using the exact coefficients G (n) , which are not available in general. The spirit behind using additional tensor bases is that increasing the degrees of freedom could partially account for the lower accuracy on the coefficients G (n) , here modelled as functions of ν t , k and ε, simply based on a dimensional argument. Several k-ε models based on Eq. (15) are available in the literature, the most common being quadratic models using four tensor bases [19,27,31] and a cubic model with eight tensor bases proposed by Craft et al. [2].
Together with the transport equations for k and ε, the definition of the eddy viscosity (14) is one of the main assumptions of k-ε models; therefore, we analyse the validity of Eq. (14) by reporting the ratio ν t /(c μ k 2 /ε) with c μ = 0.075 (Fig. 6a). This quantity allows us to assess the trend of the modelled eddy viscosity close to the wall and to define a suitable damping function f d . The eddy viscosity ν t is evaluated from DNS data using the linear assumption (9), k and ε are the turbulent kinetic energy and dissipation of DNS, and no additional transport equations are solved. Profiles are plotted as a function of the viscous scaled wall-normal coordinate y + = y/δ v (where δ v = ν/u τ is the viscous length scale), for all spanwise locations z, up to the corner bisector (inset in Fig. 6a). Duct flow data are compared to channel flow DNS data at Re τ = 1000 [1], also with c μ = 0.075. Figure 6a shows that the ratio ν t /(c μ k 2 /ε) for square duct is similar to channel flow, when duct flow data are reported up to the corner bisector. This result is consistent with Pirozzoli et al. [21], who showed that duct flow statistics are in good agreement with canonical wall-bounded flows when profiles are reported in this fashion. For this reason, we modify the damping function proposed by Rodi and Mansour [26] for plane channel flow, substituting the wall-normal distance y (in their formulation) with the distance from the closest wall d: Nisizima [18] and Mompean et al. [16] proposed a different damping function for square duct flow to account for the corner: Figure 6a shows that the damping function in Eq. (16) follows the same trend of DNS data, whereas (17) (14) is further assessed by reporting in Fig. 6b the scatter plot of the modelled eddy viscosity ν t as a function of ν t , the linear eddy viscosity from DNS (Eq. (9)). The agreement is generally satisfactory, although there are regions in which ν t is underpredicted or overpredicted with respect to DNS. These regions correspond to the corner bisector, the duct centre and the wall bisector (z/ h ≈ 0) (Fig. 7). Along the corner bisector and at the duct centre, ν t is underpredicted with respect to DNS and the modelled eddy viscosity does not predict the bulging of the isolines towards the corner. Along the wall bisector, the k-ε model predicts slightly larger values of ν t with respect to DNS, but the spatial distribution is captured. The limited accuracy along the corner bisector can be traced back to the fact that Eq. (14) has been designed for canonical wall-bounded flows (i.e. a single wall), whereas in this region the flow experiences the simultaneous influence of two walls. On the other hand, the k-ε model is more accurate along the wall bisector where the influence of the side wall is minor and the flow is more similar to canonical wall-bounded flows.
We also investigate the accuracy of different nonlinear k-ε models in Table 4 using the generalized constitutive relation (15). Figure 8 shows the scatter plot of the modelled Reynolds stresses as a function of the DNS ones for several k-ε models. The shear stress a 12 predicted by the linear k-ε model (Fig. 8a) follows DNS data with good accuracy, although larger deviations are observed compared to maximum accuracy achieved by the linear eddy viscosity hypothesis (Fig. 1a). This can be attributed to the limited validity of Eq. (14) or to the form of the damping function (16). Between these two possible causes, we slightly favour the former as the damping function (16) seems rather accurate (Fig. 6a). On the other hand, the normal stress a 22 and secondary shear stress a 23 (Fig. 8b, c) are largely different from DNS, consistent with the fact that linear models cannot predict secondary flows.
All nonlinear k-ε models return a rather accurate prediction of a 12 , without substantial improvement with respect the linear model (Fig. 8d, g, l, o). As for the prediction of a 22 and a 23 , we find large differences depending on the model. The quadratic model of Nisizima and Yoshizawa [19] (Fig. 8d-f) improves the prediction of both a 22 and a 23 , but the agreement with DNS data is much lower compared to the ideal case with four tensor bases (Fig. 1h, i) and more similar to the maximum accuracy achieved with two bases (Fig. 1e,  f). On the other hand, the models of Rubinstein and Barton [27] (Fig. 8h, i) and Shih et al. [31] (Fig. 8m, n) do not improve the accuracy of the normal stress a 22 and the secondary shear stress a 23 , despite the fact that they are built on four tensor bases, as the one by Nisizima and Yoshizawa [19]. Both the models of Shih et al. [31] and Rubinstein and Barton [27] return a 22 ≈ 0 everywhere, whereas a 23 is similar to the linear model for Shih et al. [31] and negatively correlated with DNS data for Rubinstein and Barton [27]. The cubic model of Craft et al. [2] (Fig. 8p, q) only leads to minor improvement on a 22 and a 23 .
A more quantitative insight can be gauged from Fig. 5, which allows us to compare the correlation coefficient of the models C i j with respect to the ideal one C i j . The figure shows that both the linear and nonlinear k-ε models accurately predict a 12 , with a correlation coefficient C 12 0.98. Instead, very different performance between the models is visible for the normal and secondary shear stress components. The stresses predicted by the models of Rubinstein and Barton [27] and Shih et al. [31] are poorly correlated or even negatively  (Table 3), and far from the axes bisector. A slightly higher correlation is achieved by the cubic model of Craft et al. [2] using eight tensor bases, C 22 ≈ 0.6, C 23 ≈ 0.5, but we find the highest correlation with respect to DNS for the model of Nisizima and Yoshizawa [19], C 22 ≈ 0.7 C 23 ≈ 0.6 which is also the closest to the maximum mean accuracy with four tensor bases. Figure 5 shows that the accuracy of Reynolds stresses does not only depend on the number of tensor bases in the constitutive relation (15), but is also largely affected by the model coefficients, which show a large variation between different authors ( Table 4). The analysis also shows that using a number of tensor bases larger than five does not necessarily lead to an improved accuracy. Although a larger number of coefficients increases the "tuning power" of the tensor polynomial, Eq. (15) still relies on the assumption that the coefficients G (n) ∝ ν t k/ε. In order to verify this assumption, we report G (n) /(4ν t k/ε) for n = 2, 3, 4 as a function of the distance from the closest wall d. Figure 9 shows that the hypothesis of proportionality between G (n) and ν t k/ε is partially valid for G (3) , whereas it fails for G (2) and G (4) . The model of Shih et al. [31] uses nonuniform coefficients c 1 , c 2 , c 3 , but their spatial variation is limited, and it does not follow the large nonuniformities of G (n) /(4ν t k/ε) and therefore does not result in improved performance with respect to the constant coefficients models. We also note that the model of Nisizima and Yoshizawa [19] is the only one with a value of c 2 approximately matching DNS data, thus explaining the improved accuracy of this model with respect to the other k-ε variants. We recall that the proportionality G (3) ∝ ν t k/ε should not lead turbulence modellers to develop constitutive relations using T (1) and T (3) only, as the analysis carried out using exact coefficients shows that a good subset of bases always contains T (2) and T (5) ( Table 2).
The anisotropic Reynolds stresses in a duct quadrant are also reported in Fig. 10, for the linear k-ε model and the nonlinear model of Nisizima and Yoshizawa [19]. This model improves the spatial distribution of a 22 and a 23 , but these stresses are far from the ones obtained using four tensor bases and exact coefficients (Fig. 2h,  i). The spatial organization of a 22 and a 23 obtained with the model of Nisizima and Yoshizawa [19] is much closer to the ideal distributions for N = 2 (Fig. 2e, f), as the spatial nonuniformity of a 22 close to the side wall is not captured and values of a 23 are consistently overpredicted towards the walls.

Quadratic constitutive relation
Spalart [32] proposed a quadratic constitutive relation (QCR) to correct linear eddy viscosity models and showed that it improves the prediction of turbulence secondary flows for ducts of complex cross section [33]. The idea is to add to the linear stress a nonlinear correction based on the basis T (2) in (6):  Figure 11 shows the scatter plot of the anisotropic stress from Eq. (18) as a function of the DNS stresses. The primary shear stress a 12 is in good agreement with DNS results, whereas we observe only a minor improvement in the normal stress a 22 compared to the linear model (Fig. 1b). The secondary shear stress a 23 is in good agreement with the theoretical prediction obtained using the exact coefficients G (n) and two tensor bases (Fig. 1f). Figure 5 shows that correlation coefficients C 12 = 0.997 is the same as the ideal one with two bases and C 23 = 0.691 is also close to the maximum accuracy. Instead, for the normal stress component the QCR has a mean correlation coefficient C 22 ≈ 0.6, which is lower than the ideal correlation with two tensor bases C 22 ≈ 0.9, Table 3. Additional information on the spatial distribution of these modelled stresses can be  Fig. 12. The nonlinear correction slightly improves the distribution a 22 over the linear case, but its intensity remains everywhere underpredicted with respect to DNS. The secondary shear stress a 23 has a spatial organization which is in very good agreement with the one obtained using the exact model coefficients G (n) and two tensor bases (Fig. 2f).
The analysis shows that the model proposed by Spalart [32] brings a rather accurate prediction of all stress components, in line with maximum accuracy achievable with two tensor bases, although prediction of a 22 could be improved. The model also reveals to be more accurate than most k-ε models using four bases.

v 2f model
The v 2f model was originally developed by Durbin [5] as an improvement to the classical k-ε model for flows with large inhomogeneities. The model requires transport equations for the turbulent kinetic energy k, for the dissipation rate ε and for v 2 , namely the normal stress with respect to the closest wall. Even though the model requires a transport equation for the Reynolds stress component v 2 , it is typically considered an eddy viscosity models, with eddy viscosity where T , and c v μ = 0.2 and c T = 6 are model constants. Pecnik and Iaccarino [20] proposed a variation of the linear v 2f model to take into account the anisotropy of the normal stresses in (4). The idea is to modify the linear eddy viscosity hypothesis introducing a nonisotropic contribution N i j , without relying on additional nonlinear terms, and N i j must be traceless, where  (25) and φ is the solution of the Poisson problem ∇ 2 φ = −1, with homogeneous boundary conditions. Moreover, the damping function f d in Eq. (23) is Note that this correction is referred to as "linear" by its authors, as Eq. (22) is linear with respect to S i j . Nevertheless, the accuracy of this correction cannot be compared to the maximum accuracy of the linear eddy viscosity hypothesis in Fig. 1a-c, as Eq. (22) is not based on the tensor polynomial (5). The constitutive relation (22) does not strictly belong to the same class of models addressed in the rest of the manuscript, but it has specifically been developed for duct flow and therefore assessing its accuracy can be instructive.
We first assess the accuracy of the eddy viscosity hypothesis (20) by reporting in Fig. 13a the scatter plot of ν t as a function of ν t , the linear eddy viscosity evaluated from DNS (Eq. (9)). The eddy viscosity of the v 2f model is more accurate than the one of the k-ε model, and it follows very well DNS data. Minor differences with DNS are only visible along the corner bisector (Fig. 13b, c), where the effect of secondary flows is stronger.
Despite this minor discrepancy, the primary shear stress a 12 matches DNS data almost exactly, whereas a 22 and a 23 show large scatter with respect to DNS data (Fig. 14). Large improvement in the prediction of these stress components is achieved with the linear correction of Pecnik and Iaccarino [20], which returns almost exact predictions of a 22 and a more accurate distribution of a 23 . Figure 5 confirms the qualitative impression gauged from the scatter plots. The correlation with DNS data is excellent, C i j > 0.8; nevertheless, we observe that the linear v 2f model is slightly more accurate in predicting the primary shear stress. ( C 12 = 0.997 for the linear v 2f and C 12 = 0.988 for the correction.) Visualization of the stresses in the duct quadrant confirms that this linear correction predicts spatial distributions very similar to the DNS ones (Fig. 15). The model is the only one here tested which is able to predict the spatial nonuniformity of a 22 towards the side wall (Fig. 15b) and to correctly predict the sign and intensity of a 23 in the whole field.
The success of this model can be traced back to the ability of explicitly accounting for directions of anisotropy in the flow. The model has been first designed using the plane channel flow prototype, in which the term v 2 /k accounts for the anisotropy in the wall-normal direction, together with the vector n i , which is just [0, 1, 0] for plane channel flow. Extension to multiple walls is rather trivial for the term v 2 /k, as v becomes the normal velocity component with respect to the closest wall, whereas evaluating n i requires solution of a Poisson problem.
The model is the most accurate among the ones here tested; nevertheless, the correction has only been applied to simple geometries such as square and rectangular duct; therefore, its applicability to more complex

Effect of the Reynolds number
We also consider the effect of Reynolds number on the correlation coefficient C i j , for different RANS models (Fig. 16). Contrarily to what observed for the exact constitutive relation in Fig. 3b, c, we find a more prominent Reynolds number dependence for some RANS models. The correlation coefficient of the shear stress component C 12 is practically independent from Reynolds number for all models. On the other hand, an important Reynolds number effect is observed on C 22 and C 23 for the model of Rubinstein and Barton [27], for which the anisotropic stresses decorrelate from DNS as Reynolds number increases (Fig. 16c). This is probably an indication that the model coefficients have been tuned on low Reynolds number data. We also find a milder effect of the Reynolds number on C 22 and C 23 for the models of Craft et al. [2] and Spalart [32], but in these cases the correlation coefficients tend to become independent from Re τ , as Reynolds number increases.

Conclusions
We use DNS data of square duct flow up to friction Reynolds number Re τ = 1055 from Pirozzoli et al. [21] to carry out a priori tests of algebraic Reynolds stress models and assess their ability to predict turbulent secondary flows. We perform a priori analysis of constitutive relations for eddy viscosity models which gives us useful information on the maximum accuracy achievable for a given set of tensor bases.
According to this analysis, the linear eddy viscosity hypothesis brings accurate prediction of the shear stress component a 12 , with a mean correlation coefficient with respect to DNS C 12 = 0.998, whereas the normal and secondary shear components are essentially decorrelated from reference data, C 22 = −0.044 C 23 = 0.313. Moreover, we test all possible nonlinear bases combinations and show that optimal subsets always include T (2) and T (5) , in addition to the linear term T (1) . Two tensor bases bring a fairly accurate representation of the normal and secondary shear stresses, C 22 = 0.911, C 23 = 0.743. With three bases ( C 22 = 0.910, C 23 = 0.871) and four bases ( C 22 = 1.000, C 23 = 0.938), we achieve exact or almost exact representation of the normal stress and secondary shear stress.
Unfortunately, RANS models rely on additional assumptions on the coefficients G (n) , leading to a reduced accuracy compared to the one achieved with exact coefficients. The analysis of popular k-ε models shows that only few of them are able to improve the prediction of a 22 and a 23 compared to the linear model, namely the model of Craft et al. [2] ( C 22 = 0.643, C 23 = 0.532) and of Nisizima and Yoshizawa [19] ( C 22 = 0.696, C 23 = 0.616).
Moreover, we assess the performance of the nonlinear correction of Spalart [32], based on two tensor bases. Interestingly, we find that the model is as accurate ( C 22 = 0.635, C 23 = 0.691) as the one by Nisizima and Yoshizawa [19] using four tensor bases, but the QCR is closer to its maximum accuracy. A priori analysis of nonlinear models based on Pope's polynomial expansion shows that using three or four tensor bases could potentially lead to very accurate prediction of the turbulent stresses, but this is difficult to achieve in practice as the gain in tuning power is partially spoiled by the uncertainty on the additional model coefficients. Hence, using a quadratic constitutive relation as the one proposed by Spalart [32] leads de facto to the same accuracy of higher-order models.
We also carry out a priori analysis of the v 2f model, originally developed to overcome some of the limitations of the k-ε model. The assumption ν t = c v μ v 2 T is more accurate than the one used by the k-ε model, and it returns a slightly better prediction of the primary shear stress ( C 12 = 0.997 versus C 12 = 0.988 for k-ε).
Particularly interesting is the linear correction of Pecnik and Iaccarino [20] to the v 2f model, as it shows that representation of the normal stress and secondary shear stress does not necessarily require nonlinear constitutive relations. The model is the most accurate in predicting the normal and secondary shear stress ( C 22 = 0.984 and C 23 = 0.844), but the correction does not use Pope's theoretical framework and therefore direct comparison with other nonlinear models is not possible.