2.1. Macroscopic forcing method

In this section, we discuss details on how to use MFM to measure the eddy viscosity, starting from the GMT equations.

2.1.1. GMT equation

In our earlier work, which was mainly on the transport of passive scalars, we briefly introduced how one can apply MFM to analyse momentum transport (Mani & Park Reference Mani and Park2021). To quantitatively determine the eddy viscosity operator, one first needs the detailed velocity field of the specific flow of interest. One method of obtaining such velocity fields is to perform a DNS simulation, which we call the donor simulation, as it donates a velocity field whose eddy viscosity is to be determined.

To analyse momentum transport by a given flow, we will now consider GMT, which can be derived from the Reynolds transport theorem for a fluid system with a Fickian model for molecular viscosity:

(2.5)$$\begin{gather} \frac{\partial v_i}{\partial t} + \frac{\partial u_j v_i}{\partial x_j} ={-}\frac{\partial q}{\partial x_i} + \frac{1}{{Re}} \frac{\partial^2 v_i}{\partial x_j\partial x_j} + s_i, \end{gather}$$

(2.6)$$\begin{gather}\frac{\partial v_j}{\partial x_j}=0, \end{gather}$$

where $v_i$ represents momentum per unit mass, and is considered to be different from $u_i$, the donor velocity field. Here $s_i$ is the macroscopic forcing. In addition, $q$ is the generalised pressure to ensure the incompressibility of the momentum field $v_i$.

Equations (2.5) and (2.6) then describe a passive solenoidal vector field that is transported by the background velocity field $u_j$ governed by (2.1). An advantage of working with GMT, as opposed to NS, is its linearity with respect to the transported quantity, $v_i$. Under such conditions, expressing the generalised eddy viscosity in the format given by (1.1) becomes meaningful. As discussed by Mani & Park (Reference Mani and Park2021), GMT spans a larger solution space than NS; NS is a special subset of the GMT space where $v_i=u_i$.

An important question that naturally follows is whether the computed RANS operator of GMT is the same as that of the NS equation. In our earlier work, we already showed analytically and numerically that the macroscopic operators of the GMT and NS equations are identical (Mani & Park Reference Mani and Park2021). In brief, we showed that the solutions of GMT and NS equations become microscopically the same after sufficient time regardless of the initial conditions when we apply the same boundary conditions to both equations. The time scale at which the solutions become identical was found to be $\tau _{mix} = 16.6 \delta / u_\tau$ for a turbulent channel flow. Therefore, it is justified that the macroscopic operator of the GMT equation obtained by MFM is the same as the RANS operator of the NS equations. In sum, GMT works as an auxiliary set of equations that probes RANS operator of NS and therefore we can obtain eddy viscosity of the RANS equations by investigating that of the GMT equations.

It is important to note that Hamba (Reference Hamba2005) wrote an equation very similar to GMT equations in spite of taking a conceptually different derivation path. His passive vector equation is indeed GMT subtracted by the mean of GMT. The main difference lies in the explicit inclusion of forcing in the equations, allowing for a general macroscopic field. In contrast, Hamba (Reference Hamba2005) implicitly applies forcing by specifically considering Dirac delta function mean fields.

2.1.2. Analysis strategy

We aim to study two aspects of the eddy viscosity kernel in a turbulent channel flow: the anisotropy and the non-locality. To fully investigate such non-Boussinesq effects, it is ideal to compute every value of the full eddy viscosity kernel ${\mathsf{D}}_{ijkl}$ in (1.1). Since the channel flow is homogeneous in $x_1$ and $x_3$ directions and statistically stationary, we integrate the mixing effect in these directions. The simplified Reynolds stress for GMT variables can be expressed as

(2.7)\begin{equation} -\overline{u_i'v_j'}(x_2)=\int {\mathsf{D}}_{ijkl}({x_2}, {y_2})\left.\frac{\partial V_l}{\partial x_k} \right\vert_{{y_2}} \,\mathrm{d}y_2. \end{equation}

Equation (2.7) incorporates anisotropy via tensorial representation and non-locality via the integration form. MFM has the capability to compute all the elements in the eddy viscosity kernel ${\mathsf{D}}_{ijkl}(x_2,y_2)$ by tracking the influence of each entry of ${\rm d}V_l/{\rm d}\kern0.7pt x_k$ on the entire Reynolds stress field. It has been demonstrated by Liu, Williams & Mani (Reference Liu, Williams and Mani2023) that such a brute force approach is theoretically equivalent to Hamba's Green's function approach (Hamba Reference Hamba2005).

However, one caveat is that the cost of each simulation is significant and consequently it is not desirable to conduct a full non-local MFM analysis. To conduct computation for ${\mathsf{D}}_{ijkl}$ for given $k$ and $l$, one requires as many DNS simulations as the number of degree of freedom of the RANS space. Therefore, to reduce the cost of the analysis, we conduct two separate analyses for the anisotropy and non-locality, both using MFM.

First, we focus on studying the anisotropic nature of the eddy viscosity. However, to focus exclusively on anisotropy, we systematically construct a local approximation of the eddy viscosity operator using the Kramers–Moyal expansion (Van Kampen Reference Van Kampen1992), as investigated by Mani & Park (Reference Mani and Park2021). For instance, in a parallel flow where ${\rm d}V_1/{\rm d}\kern0.7pt x_2$ is the only active component of the velocity gradient, the Reynolds stress $\overline {u_2^\prime v_1^\prime }$ in (2.7) can be written as the integral of only ${\mathsf{D}}_{2121}$ component of the eddy viscosity. By considering a Taylor series expansion of ${\rm d}V_1/{\rm d}\kern0.7pt x_2$ around $y_2=x_2$, one can re-express the eddy viscosity operator in terms of the following expansion:

(2.8)\begin{align} -\overline{u_2'v_1'}(x_2) &= \int {\mathsf{D}}_{2121}({x_2}, {y_2}) \left. \frac{\partial V_1}{\partial x_2} \right|_{y_2} \,{\rm d}{y_2} \end{align}

(2.9)\begin{align} &= \int {\mathsf{D}}_{2121}({x_2}, {y_2}) \left( \left. \frac{\partial V_1}{\partial x_2} \right|_{x_2} + (y_2-x_2) \left. \frac{\partial^2 V_1}{\partial x_2} \right|_{x_2} + \cdots \right) {\rm d}{y_2} \end{align}

(2.10)\begin{align} & = \sum_{n=0}^{\infty} {\mathsf{D}}^n_{2121}({x_2}) \frac{\partial^{n+1} V_1}{\partial x_2^{n+1}}, \end{align}

where ${\mathsf{D}}^n_{2121}=\int {{\mathsf{D}}_{2121}(x_2,y_2)(y_2-x_2)^n/n!\,{\rm d}y_2}$ represents the $n$th spatial moment of the eddy viscosity kernel.

As discussed by Mani & Park (Reference Mani and Park2021), the leading term in this expansion encapsulates the local limit eddy viscosity while the subsequent terms characterise finite moments associated with the non-local effects. The general form of this leading-order approximation for all components of the Reynolds stress and mean velocity gradient is as follows:

(2.11)\begin{equation} -\overline{u_i'v_j'} (x_2) = {\mathsf{D}}^0_{ijkl} (x_2) \frac{\partial V_l}{\partial x_k}, \end{equation}

where ${\mathsf{D}}^0_{ijkl}(x_2)$ is called the leading-order eddy viscosity tensor,

(2.12)\begin{equation} {\mathsf{D}}^0_{ijkl}(x_2)=\int{{\mathsf{D}}_{ijkl} \,{\rm d}y_2}. \end{equation}

Equation (2.11) would be exact only when ${\mathsf{D}}_{ijkl}(x_2, y_2)$ is local, i.e. ${\mathsf{D}}_{ijkl}(x_2, y_2) = {\mathsf{D}}^0_{ijkl}(x_2) \delta (y_2 - x_2)$ where $\delta (x)$ is a Dirac delta function.

The local eddy viscosity in (2.11) is no longer a scalar value varying in space; it is a fourth-order tensor with 81 coefficients. The tensor representation was suggested by previous researchers including Batchelor (Reference Batchelor1949), but the full quantification has not been conducted to the best of the authors’ knowledge. As presented in Appendix B, using only 9 MFM simulations, we computed all 81 coefficients of the eddy viscosity tensor. The resulting tensor elements are provided in Appendix C.

The next investigation focuses on the non-locality of the eddy viscosity. As conducting MFM to measure the full kernel can be costly for complex turbulent flow systems, we focus on calculating a subset of tensorial kernel components, specifically the kernel components that are multiplied to $\partial V_1/\partial x_2$ in (2.7). The computed tensorial kernel components are ${\mathsf{D}}_{ij21}(x_2,y_2)$ and they are associated with the Reynolds stresses which correspond to the velocity gradient $\partial U_1/\partial x_2$, the only velocity gradient appearing in the RANS closure for a channel flow. The detailed steps on how to measure eddy viscosity kernel using MFM is discussed in Appendix E.

2.1.3. Application of MFM

The next step involves how we actually compute the leading-order eddy viscosity tensor and the eddy viscosity kernel. Figure 2 illustrates how we conducted our MFM analysis. To apply MFM, we start with two sets of solvers: one for the NS equations and the other for GMT. At each time step, we solve the NS equation to obtain the velocity field $u_i$ and feed it as the advecting velocity to the GMT solver. For the GMT equations, we force the Reynolds-averaged GMT variable $V_i$ to be a specific value in order to acquire certain information about the eddy viscosity. A forcing field $s_i$ that results in $V_1=x_2$ and $V_2=V_3=0$ generates GMT data from which we can extract the leading-order eddy viscosity ${\mathsf{D}}^0_{ij21}$. More specifically, ${\mathsf{D}}^0_{2121}$ can be obtained by post-processing $\overline {u'_2 v'_1}$ from this GMT simulation, and re-evaluating (2.8)–(2.10) to observe

(2.13)\begin{align} -\overline{u'_2 v'_1}(x_2) = \int {\mathsf{D}}_{2121}({x_2}, {y_2}) \left. \frac{\partial V_1}{\partial x_2} \right|_{y_2} \,{\rm d}{y_2} \stackrel{V_1 = x_2}{=} \int {\mathsf{D}}_{2121}({x_2}, {y_2}) \,{\rm d}{y_2} = {\mathsf{D}}^0_{2121}(x_2). \end{align}

Figure 2. Schematics of the MFM analysis.

As discussed, the macroscopic operator of the Reynolds-averaged GMT and the RANS operator of NS are identical. Therefore, ${\mathsf{D}}^0_{2121}$ corresponds to the standard eddy viscosity $\nu _T$ used in the Boussinesq approximation in RANS models. Likewise, we can compute other components of the leading-order eddy viscosity tensor using different selections of the macroscopic forcing field, $s_i$, such that other components of the mean velocity gradient are activated.

In addition, the same set-up shown in figure 2, can be used to compute the full kernel of eddy viscosity. The main difference is to apply macroscopic forcings that would generate mean gradient fields in the form of Dirac delta functions. For example, a macroscopic field, $s_i$, that sustains $\partial V_1/\partial x_2=\delta (x_2-y_2^*)$, would result in GMT data from which we can extract ${\mathsf{D}}_{ij21}(x_2,y_2^*)$, by merely post-processing the $\overline {u'_i v'_j}$ data. This specific choice of forcing would result in data similar to those obtained by Hamba (Reference Hamba2005), with the difference that Hamba used only the symmetric portion of the momentum flux tensor in order to ensure symmetry of the Reynolds stresses. As we shall see, GMT does not produce symmetric eddy viscosity kernels and, thus, ${\mathsf{D}}_{ijkl} \ne {\mathsf{D}}_{jikl}$. This is intuitively understandable noting that ${\mathsf{D}}_{ijkl}$ quantifies the rate of mixing of the mean $j$-momentum by the $i$-component of the velocity fluctuations while ${\mathsf{D}}_{jikl}$ quantifies the rate of mixing of the mean $i$-momentum by the $j$-component of the velocity fluctuations. Since in this framework, momentum and velocity fields can be quantitatively different, the symmetry does not hold. Likewise, this asymmetry propagates to the Kramers–Moyal expansion of the eddy viscosity operator, and as we shall see, even the leading-order eddy viscosities are not symmetric.

Lastly, we note that the macroscopic forcing procedure used in this work is an inverse forcing method as discussed by Mani & Park (Reference Mani and Park2021), since we explicitly set the desired mean momentum field $V_i$ for each GMT simulation, as opposed to setting the macroscopic forcing field.

2.2. Simulation set-up

We adapt MFM solver to a three-dimensional incompressible NS solver originally developed by Bose, Moin & You (Reference Bose, Moin and You2010) and modified by Seo, García-Mayoral & Mani (Reference Seo, García-Mayoral and Mani2015). The present DNS uses the fractional step method with semi-implicit time advancement (Kim & Moin Reference Kim and Moin1985). For the temporal difference scheme, we use second order Crank–Nicholson for the wall-normal diffusion and Adams–Bashforth for the rest of the terms. The solver uses a second-order finite spatial discretisation on a staggered mesh (Morinishi et al. Reference Morinishi, Lund, Vasilyev and Moin1998). In addition, we use a uniform grid in the streamwise and spanwise directions and grid-stretching in the wall-normal direction. The domain is periodic both in the spanwise and the streamwise directions, and the no-slip boundary condition is applied at the two walls.

The numerical set-up for the GMT solver is almost identical to that of DNS, except for two differences. The first is that GMT obtains the background velocity from the NS solver at every time step. The other difference is that GMT utilises macroscopic forcing, in order to maintain a desired macroscopic momentum field $V_i(x_2)$. To be most rigorous, the selected macroscopic forcing, $s_i(x_2)$, must be independent of time. Likewise, the resulting mean velocity field needs to match the pre-set $V_i(x_2)$ only after time averaging. However, constraining the simulations in this fashion, would require expensive iterations over which the entire simulation must be repeated after each adjustment of $s_i(x_2)$. To avoid this cost, in our implementation, we computed ensemble averages by averaging fields only in the $x_1$ and $x_3$ directions, and we constrained $s_i(x_2)$ at each time step such that $V_i(x_2)$ is matched to the pre-set $V_i(x_2)$. In other words, one can set $s_i(x_2)$ such that the homogeneous plane average of the temporal term ${\partial v_i}/{\partial t}$ in (2.5) becomes zero at each timestep.

However, in this implementation the resulting $s_i(x_2)$ is not perfectly time independent. Due to finite number of samples per time step, fluctuations in time are observed. One remedy to reduce these fluctuations is to increase the number of samples by selecting a longer domain in the $x_1$ and $x_3$ directions. We have performed such domain convergence studies in Appendix A indicating the adequacy of the selected domain size in our MFM analysis.

There are two sets of forcings for MFM presented in this paper, each corresponding to the analysis of anisotropy and the non-locality of eddy viscosity (table 1). Within each set, multiple simulations are performed where the macroscopic forcings are varied to reveal different components of the eddy viscosity. The first set uses GMT simulations under different macroscopic forcings to reveal the leading-order eddy viscosity tensor ${\mathsf{D}}^0_{ijkl}$. We utilise these measurements to understand the anisotropy of the eddy viscosity. The second set probes a subset of the entire eddy viscosity kernel, ${\mathsf{D}}_{ij21}$, which quantifies the non-locality of the eddy viscosity in response to the most significant velocity gradient $\partial U_1 / \partial x_2$. In addition to the analysis method and the resulting eddy viscosity, table 1 presents the number of total DNSs in each set, the domain size, the spatial resolution and the sampling times. For the first set, only nine DNSs are needed corresponding to $k, l \in \{1, 2, 3\}$, and for the second set, MFM analyses require a set of simulations with the number of the macroscopic degrees of freedom. The results of each set are discussed in §§ 3 and 4, respectively, and the detailed simulation set-up of each set is discussed in Appendices B and C, respectively. In addition, the measured eddy viscosities, ${\mathsf{D}}^0_{ijkl}$ and ${\mathsf{D}}_{ij21}$, are provided as the supplementary data.

Table 1. Simulation set-up for anisotropy analysis and non-locality analysis of the eddy viscosity. Here $L_1 \times L_2 \times L_3$ is the domain size, $N_1 \times N_2 \times N_3$ is the number of grids and $T u_\tau /L_2$ is the total simulation time in wall units.

3.1. Standard eddy viscosity

In parallel flows, among all the components of the eddy viscosity tensor, by far the most important component is ${\mathsf{D}}^0_{2121}$ which represents the mixing effect by $\partial {U_1}/\partial {x_2}$. This component also corresponds to the standard eddy viscosity $\nu _T$. Figure 3 shows the MFM-measured ${\mathsf{D}}^0_{2121}$ across the wall-normal dimension $x_2$. An important observation here is that the MFM allows us to measure the eddy viscosity at the channel centre plane $x_2 = 0$, where the velocity gradient $\partial {U_1}/\partial {x_2}$ is zero due to the symmetry of the mean velocity profile. This value is hard to obtain in typical approaches: tuning $\nu _T$ to $\overline {u'_2 u'_1} / (\partial {U_1}/\partial {x_2})$. The eddy viscosity can be numerically determined by analysing the sampling points converging towards the centreline. However, in the vicinity of the centreline, both the numerator and denominator of the expression are significantly influenced by statistical noise. To achieve a reliable estimate, an extensive period of time integration is necessary to reduce the noise to acceptably low levels.

Figure 3. Eddy viscosity element ${\mathsf{D}}^0_{2121}$.

Figure 4 shows instantaneous field data for the normalised streamwise velocities $u_1'$ and $v_1'$ of the MFM simulation for evaluation of ${\mathsf{D}}^0_{2121}$ at the same instantaneous time. Figures (a,b) show the velocity profile over $(x_1, x_3)$ cross-section taken at $x_2 = -0.8492$ ($x_2^+ = 27$) and figures (c,d) show the velocity profile over $(x_3, x_2)$ cross-section taken at $x_1 = 3.08$. The key feature shown is that even though the forcings for the NS vector field $u_i$ and the GMT vector field $v_i$ are completely different macroscopically, MFM leads to similar features in the $u'_1$ and $v'_1$. The same qualitative observation holds across all three components of the $u'_i$ and $v'_i$ fields. Furthermore, while for $x_2<0$ we observe positive correlation between $u'_i$ and $v'_i$ fields, the sign of correlation flips for $x_2>0$. For this specific MFM analysis, the sole difference between $u'_i$ and $v'_i$ fields is in the enforced mean velocity profile. As shown by Mani & Park (Reference Mani and Park2021) without forcing, GMT would result in $v$-fields identical to $u$-fields after a few flow through times regardless of the choice of initial conditions. The case shown in figure 4 corresponds to a forced GMT in which the mean velocity gradient is kept constant $\partial V_1/\partial x_2=1$ in order to examine mixing by the leading-order (local limit) eddy viscosity. The observation in figure 4 suggests that mixing of the streamwise momentum in turbulent channel flow is substantially influenced by the leading-order effects.

Figure 4. Instantaneous velocity contours $u'_1$ and $v'_1$ normalised with each maximum value, $u'_1/(2\max (u'_1))$ and $v'_1/(2\max (v'_1))$: (a,b) correspond to the cross-section taken at $x_2 = -0.8492$ and (c,d) correspond to the cross-section taken at $x_1=3.08$. The shown vector field $v_i$ corresponds to a leading-order MFM in which the GMT equation is macroscopically forced to achieve $V_1=x_2$ and $V_2=V_3=0$. (a) Normalised $u'_1$ in $x_1$–$x_3$ plane. (b) Normalised $v'_1$ in $x_1$–$x_3$ plane. (c) Normalised $u'_1$ in $x_3$–$x_2$ plane. (d) Normalised $v'_1$ in $x_3$–$x_2$ plane.

To assess this conclusion quantitatively we next obtain the RANS solution using the measured ${\mathsf{D}}^0_{2121}$ to examine how accurate the leading-order eddy viscosity performs for the prediction of the mean velocity profile. Since the prediction of the mean channel flow only requires one component in the Reynolds stress, we conduct the RANS simulation using ${\mathsf{D}}^0_{2121}$ and compare the predicted solution with that of the DNS. As shown in figure 5, the MFM-based leading-order RANS solution predicts the DNS solution very accurately with an accuracy of 99 %. The accuracy is computed with max absolute error $\max (U_1^{{DNS}}-U_1^{{MFM}})/\max (U_1^{{DNS}})$, where $U_1^{{DNS}}$ is the streamwise velocity from DNS and $U_1^{{MFM}}$ is the streamwise velocity predicted from RANS using MFM-measured eddy viscosity ${\mathsf{D}}^0_{2121}$. The accuracy of this local RANS prediction indicates that the mean momentum mixing in the turbulent channel flow might be local. To assess this conclusion with more certainty we will later directly examine the non-locality of the eddy viscosity kernel.

Figure 5. Plot of $U_1$ reconstructed from ${\mathsf{D}}^0_{2121}$ and ${\mathsf{D}}^{0}_{2121}$ in wall units: (a) RANS prediction using ${\mathsf{D}}^0_{2121}$ where the dotted green line is the mean velocity prediction using ${\mathsf{D}}^0_{2121}$ and the blue solid line is its comparison to the DNS data; and (b) ${\mathsf{D}}^{0}_{2121}$ in wall-units where ${\mathsf{D}}^{0+}_{2121}={\mathsf{D}}^{0}_{2121}{Re}$ and $x_2^+=x_2{Re}/\delta$.

Hamba (Reference Hamba2005) also reported a small subset of the components of the leading-order (local limit) eddy viscosity, through a more expensive method of first computing the full eddy viscosity kernel for those components and then performing integration as in (2.12). Our result in figure 5 regarding accuracy of the leading-order eddy viscosity is in contrast to his result (see figure 4 in Hamba Reference Hamba2005). We attribute this difference to the fact that Hamba used the average of ${\mathsf{D}}^0_{2121}$ and ${\mathsf{D}}^0_{1221}$ as the representative local eddy viscosity. This averaging was motivated to enforce symmetric Reynolds stresses. However, conceptually these two eddy viscosities represent different mixing rates: the former represents mixing of the streamwise momentum in the wall-normal direction, whereas the latter represents mixing of the wall-normal momentum in the streamwise direction. As we shall see, while a full eddy viscosity kernel reproduces symmetric Reynolds stresses, the leading-order eddy viscosity causes errors not only in magnitude but also in symmetry of Reynolds stresses.

We next use MFM to quantify other components of ${\mathsf{D}}^0_{ijkl}$. Although these components do not affect prediction of the mean velocity profile in purely parallel flows, they provide an understanding of momentum mixing by this parallel flow, if hypothetical mean momentum gradients were imposed in other directions. One motivation to study these additional components of ${\mathsf{D}}^0_{ijkl}$ is to provide reference data of closure operators, as opposed to closure terms, for models that offer anisotropic eddy viscosity. Our analysis is additionally motivated by observation of spatially developing attached turbulent boundary layers, where weak momentum gradient could exist in both streamwise and spanwise directions. These mean gradients induce additional Reynolds stresses, due to components of ${\mathsf{D}}_{ijkl}$ other than ${\mathsf{D}}_{2121}$. In addition, it has been observed that turbulent boundary layers have similar hairpin structures in their velocity field as those seen in turbulent channel flows (Eitel-Amor et al. Reference Eitel-Amor, Örlü, Schlatter and Flores2015), and thus are expected to mix momentum in manners qualitatively similar to that of a turbulent channel flow. While quantitative differences are expected between the two flows, we expect anisotropy in eddy viscosity observed in turbulent channel flow be at least qualitatively representative of anisotropy encountered in wall-attached turbulent boundary layers in the absence of substantial wall curvature. Some of these qualitative similarities, such as components in ${\mathsf{D}}_{ijkl}$ with the highest magnitude, can already be confirmed from the study of Park, Liu & Mani (Reference Park, Liu and Mani2022) with a specific focus to their analysis of pre-separation zone of turbulent boundary layers. However, given the stringent statistical convergence requirements for MFM simulations, e.g. at least an order of magnitude longer simulations needed than commonly reported DNS, compared to turbulent boundary layers, turbulent channel flows have the advantage of cheaper runtime per time step and availability of additional homogeneous direction for statistical convergence.

3.2. Quantifying anisotropy

We computed all other components of the anisotropic eddy viscosity tensor ${\mathsf{D}}^0_{ijkl}$, a total of 81 coefficients as a function of the wall-normal coordinate. All the data are shown in Appendix C. Out of 81 components, 41 are non-zeros and 40 are inevitably zero due to the symmetry in spanwise direction.

Out of all the elements, the largest eddy viscosity component is ${\mathsf{D}}^0_{1111}$, with a maximum value of 1.318, and the smallest non-zero eddy viscosity component is ${\mathsf{D}}^0_{2331}$ with the maximum value of 0.00248. After comparing these values to a maximum value of the nominal eddy viscosity ${\mathsf{D}}^0_{2121}$, which is 0.0767, we determined that the largest coefficient in the eddy viscosity tensor is one order of magnitude larger than the nominal eddy viscosity and three orders of magnitude larger than the smallest coefficient, indicating a significant anisotropy. When we examine these ratios locally at each $x_2$, the differences are more drastic and may go up to a few orders of magnitude. After ${\mathsf{D}}^0_{1111}$, the largest eddy viscosity components are ${\mathsf{D}}^0_{1212}$ and ${\mathsf{D}}^0_{1313}$ with maximum values of 0.573 and 0.407, respectively. All three eddy viscosities have their first and third index represented by the streamwise direction. These indices represent the component of the velocity field that mixes momentum and the direction of the mean-momentum gradient, respectively. This observation coincides with the fact that $u'_1$ is the largest fluctuating velocity component in channel flow. Combining the two observations, we conclude that $u'_1$ is the strongest mixer of momentum and is most effective in mixing in the $x_1$ direction, as intuitively expected. Specifically, the rate of momentum mixing in the streamwise direction is substantially faster than the standard eddy viscosity which characterises the rate of (streamwise) momentum mixing in the wall-normal direction.

In addition, all three dominant eddy viscosity components have repeated second and fourth indices. These indices respectively represent the momentum component that is being mixed and the momentum component whose mean gradient is responsible for mixing. Based on this observation, we conclude that within ${\mathsf{D}}^0_{1j1l}$, mean gradient of component $l$ most effectively contributes to the generation of $\overline {u'_1 v'_j}$ when $j=l$. In other words, gradient of each momentum component most effectively generates fluxes of the same momentum component at least in the leading-order limit. This latter observation is extendable to ${\mathsf{D}}^0_{ijil}$ components, and is not a surprising outcome given that the production term in the transport equation for $\overline {u'_i v'_j}$ involves the mean gradient of $V_j$.

As we discussed previously, since flow structures and thus momentum mixing is similar between the channel flow and the attached boundary layers, we can use the measured eddy viscosity anisotropy in the former setting to identify important eddy viscosity components for the latter setting. To this end, we present in Appendix D a scaling analysis of various gradients contributing to the Reynolds stress tensor. Combining this analysis with the measured order of magnitude of each eddy viscosity component that acts as a pre-factor multiplying components of the velocity gradient tensor, we identify the key eddy viscosity components that contribute dominantly to the Reynolds stress tensor budget. Based on our analysis we identify ${\mathsf{D}}_{1111}$, ${\mathsf{D}}_{1121}$, ${\mathsf{D}}_{2121}$ and ${\mathsf{D}}_{2221}$ as the key four, out of 16, dominant eddy viscosity components for two-dimensional spatially developing turbulent boundary layers.

Motivated by this example, we next examine the identified anisotropy against the Boussinesq approximation. When we cast the Boussinesq approximation to our tensorial representation, the components in the eddy viscosity tensor are in ratio of 0, 1 or 2 to the standard eddy viscosity $\nu _T$. For instance, the four elements are prescribed with following ratios; ${\mathsf{D}}_{1111}=2\nu _T$, ${\mathsf{D}}_{1121}=0$, ${\mathsf{D}}_{2121}=\nu _T$ and ${\mathsf{D}}_{2221}=0$. Figure 6 shows a comparison of these eddy viscosity components to the Boussinesq approximation. In figure 6(a), we show the measured four elements using our MFM calculation. In figure 6(b) we set the standard eddy viscosity to the MFM-measured leading-order value, $\nu _T={\mathsf{D}}^0_{2121}$ and prescribe the other components with the ratio to $\nu _T$. As shown in the figure, a huge anisotropy is observed not only among all elements but specifically among these four critical elements, and the ratio of these plots can locally go up to hundreds. We conclude that whereas ${\mathsf{D}}_{2121}$ is the most important eddy viscosity component for parallel and semi-parallel flows, the presence of small non-parallel effects could lead to significant influence of anisotropy in momentum transport in wall bounded flows.

Figure 6. Comparison of the measured eddy viscosity elements ${\mathsf{D}}^0_{1111}$ (blue line), ${\mathsf{D}}^0_{1121}$ (orange line), ${\mathsf{D}}^0_{2121}$ (green line) and ${\mathsf{D}}^0_{2221}$ (red line) to the Boussinesq approximation.

Lastly, we point out that there have been attempts to include the anisotropy in RANS such as Spalart–Allmaras model with quadratic constitutive relation (SA-QCR) (Spalart Reference Spalart2000; Mani et al. Reference Mani, Babcock, Winkler and Spalart2013; Rumsey et al. Reference Rumsey, Carlson, Pulliam and Spalart2020). However, examining our results suggest that these models do not captured the level of the anisotropy that MFM measured. For instance, SA-QCR still prescribes ${\mathsf{D}}_{1111}=2\nu _T$ and the anisotropy is not yet introduced in needed directions.

3.3. Dependence on the rate of rotation

The dependence of Reynolds stress on rate of rotation is studied previously in the literature. Key main examples are (1) investigation of environmental flows that include the Coriolis term (Speziale, Raj & Gatski Reference Speziale, Raj, Gatski, Gatski, Speziale and Sarkar1992) and (2) corrections to the Boussinesq eddy viscosity based on Cayley–Hamilton theorem (Pope Reference Pope1975), which are also incorporated in quadratic constitutive relation (QCR) models (Spalart & Allmaras Reference Spalart and Allmaras1994). However, these models incorporate the mean rotation effects as higher-order nonlinear corrections to the Boussinesq eddy viscosity.

With our eddy viscosity tensor notation, insensitivity of Reynolds to the mean rotation implies that ${\mathsf{D}}^0_{ijkl}$ must be equal to ${\mathsf{D}}^0_{ijlk}$ because under this condition, each Reynolds stress component, $\overline {u'_i u'_j}$, would be equally sensitive to both $\partial U_l/\partial x_k$ and $\partial U_k/\partial x_l$, and thus is a function of the summation $\partial U_l/\partial x_k + \partial U_k/\partial x_l$, which is $2S_{kl}$. However, our measurement of the leading-order eddy viscosity tensor invalidates the relation ${\mathsf{D}}^0_{ijkl} {\mathsf{D}}^0_{ijlk}$. Figure 7 shows the comparison between ${\mathsf{D}}^0_{2121}$ and ${\mathsf{D}}^0_{2112}$. These two components have the same sign and their qualitative shape is similar, but the magnitudes are drastically different. This highlights an important conclusion: sensitivity of Reynolds stresses on mean rotation is not a secondary or higher-order effect and is present even at the leading-order term of the eddy viscosity expansion. Likewise, we reach the same conclusion with the case of ${\mathsf{D}}^0_{ij13}$ and ${\mathsf{D}}^0_{ij23}$.

Figure 7. Comparison of ${\mathsf{D}}^0_{2121}$ (blue line) and ${\mathsf{D}}^0_{2112}$ (green line).

3.4. Leading-order Reynolds stress

In § 3.1, we computed a RANS solution using ${\mathsf{D}}^0_{2121}$ and compare the solution with that of the DNS to assess appropriateness of the leading-order eddy viscosity for prediction of the mean velocity profile. Another way to make this assessment is to reconstruct Reynolds stress using the computed eddy viscosity tensor and compare it with the Reynolds stress of DNS. In this way, we can assess some other components in ${\mathsf{D}}^0_{ijkl}$. The Reynolds stress in the channel flow can be represented in the following way: $\overline {u'_i u'_j}^0=-{\mathsf{D}}^0_{ij21} {\partial U_1}/{\partial x_2}$, with the leading-order eddy viscosity tensor ${\mathsf{D}}^0_{ij21}$ computed using MFM and with the mean velocity gradient ${\partial U_1}/{\partial x_2}$ measured from the DNS data. Here, the superscript zero is added to the Reynolds stress to indicate that this is the leading-order reconstruction.

Figure 8 shows the five reconstructed Reynolds stresses associated with the RANS prediction of the channel flow, in comparison with the Reynolds stresses from the DNS data. There are three important observations with the Reynolds stresses that are reconstructed with the leading-order eddy viscosity tensor. The first finding is that while the Reynolds stresses reconstructed using the leading-order eddy viscosity show similar qualitative trends and orders of magnitudes to those from DNS, there is still a noticeable difference between the two. This difference is likely due to the leading-order truncation of the eddy viscosity operator. Among various components of the Reynolds stress tensor, only $\overline {u'_2 u'_1}^0$ is reasonably constructed. This observation is compatible with the previous observation of reasonable RANS solution for the mean flow, since only this Reynolds stress component is involved in mean momentum mixing for this flow.

Figure 8. Reynolds stresses constructed by the leading-order eddy viscosity tensor associated with (a) ${\mathsf{D}}^0_{1121}$, (b) ${\mathsf{D}}^0_{1221}$, (c) ${\mathsf{D}}^0_{2121}$, (d) ${\mathsf{D}}^0_{2221}$ and (e) ${\mathsf{D}}^0_{3321}$: green dotted line, the reconstructed Reynolds stress by the leading-order eddy viscosity tensor $\overline {u'_i u'_j}^0$; blue solid line, the DNS data.

The next important observation is that constructed Reynolds stresses from the leading-order eddy viscosity are not symmetric. When the eddy viscosity operator is modelled through a truncated Kramers–Moyal expansion, we introduce an approximation that inherently involves a loss of symmetry. This asymmetry arises from employing the GMT equation to examine momentum transport, leading to the observation that $\overline {u'_2 v'_1}^0$ is not equal to $\overline {u'_1 v'_2}^0$. However, considering that the NS equations serve as the attractor solution to the GMT framework, as established by Mani & Park (Reference Mani and Park2021), utilising exact closure operators ensures that $\overline {u'_i v'_j}$ aligns with $\overline {u'_i u'_j}$, thereby preserving the symmetry of the tensor. As we shall see, inclusion of the full non-local eddy viscosity will eliminate this error. However, the fact that the leading-order $\overline {u'_2 u'_1}^0$ match the DNS substantially better than $\overline {u'_1 u'_2}^0$ indicates that the former Reynolds stress is more local while the latter has substantial non-local sensitivity to the mean velocity gradient. Although we do not have an intuitive explanation for this observation, we note the coincidence that the former Reynolds stress represents flux of an active mean momentum component in the direction where its gradients are active. The only way that the latter Reynolds stress could be generated in this setting is through pressure coupling, whose fluctuations are known to non-locally depend on velocity fluctuations.

Lastly, we observe that the leading-order eddy viscosity cannot reproduce the non-zero Reynolds stresses at the centreline, where the velocity gradient is zero due to the symmetry of the channel flow. Non-locality needs to be included in eddy viscosity to enable prediction of non-zero Reynolds stresses in regions of zero mean velocity gradient.

3.5. Positive definiteness

It is noted that the Reynolds stress is a positive semi-definite tensor (Du Vachat Reference Du Vachat1977; Schumann Reference Schumann1977). We often require eddy viscosity to satisfy the same condition as done in the Boussinesq approximation with $\nu _T \geqslant 0$ (Speziale, Abid & Durbin Reference Speziale, Abid and Durbin1994) for a well-posed closure model. In this section, we discuss whether this condition holds for our leading-order eddy viscosity tensor ${\mathsf{D}}^0_{ijkl}$ as well.

The positive definiteness of the eddy viscosity is closely related to the mean kinetic energy equation, which is the following:

(3.1)\begin{align} \frac{\partial}{\partial t}\left(\frac{U_i U_i}{2}\right) + U_j\frac{\partial}{\partial x_j}\left(\frac{U_i U_i}{2}\right) &= \frac{\partial}{\partial x_j}\left(-\frac{P}{\rho}U_j\right) + \nu \frac{\partial^2 U_i U_i / 2}{\partial x_j \partial x_j} \end{align}

(3.2)\begin{align} &\quad - \nu \frac{\partial U_i}{\partial x_j} \frac{\partial U_i}{\partial x_j} - \frac{\partial}{\partial x_j}\left( U_i \overline{u_j^\prime u_i^\prime} \right) + \overline{u_j^\prime u_i^\prime} \frac{\partial U_i}{\partial x_j}. \end{align}

The last term in (3.2) is the negative of turbulent kinetic energy production. It is well-known that this term drains the kinetic energy from the mean flow via interactions of the mean shear and the turbulent fluctuations, and provides energy to the turbulence production. We denote the turbulent kinetic energy production as $P_k = -\overline {u_j^\prime u_i^\prime } ({\partial U_i}/{\partial x_j})$. In all statistically stationary flows, the volumetric integral of $P_k$ must be non-negative, otherwise turbulent kinetic energy cannot be sustained. There are certain cases such as the separation of the shear layer where $P_k$ is locally negative, but even for those cases, the turbulent production is positive for most of the domain (Cimarelli et al. Reference Cimarelli, Leonforte, De Angelis, Crivellini and Angeli2019) rendering the total volume integral positive. The volumetric integral condition for $P_k$ can be expressed using our generalised eddy viscosity expression in (1.1):

(3.3)\begin{align} \int{P_k}\,{\rm d}^3{\boldsymbol{x}} &= \int -\overline{u_j^\prime u_i^\prime} \left. \frac{\partial U_i}{\partial x_j}\right\vert_{\boldsymbol{x}} \,{\rm d}^3{\boldsymbol{x}} , \end{align}

(3.4)\begin{align} &= \int\int {\mathsf{D}}_{ijkl}(\boldsymbol{x}, \boldsymbol{y}) \left. \frac{\partial U_l}{\partial x_k} \right\vert_{\boldsymbol{y}} \left. \frac{\partial U_i}{\partial x_j}\right\vert_{\boldsymbol{x}} \,{\rm d}^3{\boldsymbol{y}} \, {\rm d}^3{\boldsymbol{x}} \geqslant 0. \end{align}

The last relation is the same statement as conditioning the eddy viscosity operator to be positive semi-definite. For well-posedness of its RANS mathematical model, any given eddy viscosity field must satisfy this condition for all arbitrary admissible input mean velocity gradients $\partial U_i/\partial x_j$. Otherwise, there will be unstable modes of mean flow that can be energised by the turbulence model, leading to their time exponential blow-up.

Equation (3.4) can be further simplified to a single spatial integral when the eddy viscosity operator is local ${\mathsf{D}}_{ijkl}(\boldsymbol {x}, \boldsymbol {y})={\mathsf{D}}^0_{ijkl}(\boldsymbol {x})\delta (\boldsymbol {y}-\boldsymbol {x})$ such as in the case of the leading-order eddy viscosity. Substitution of a local model in (3.4) results in

(3.5)\begin{equation} \int{P_k}\,{\rm d}^3{\boldsymbol{x}} \simeq \int {\mathsf{D}}^0_{ijkl}(\boldsymbol{x}) \left. \frac{\partial U_l}{\partial x_k} \right\vert_{\boldsymbol{x}} \left. \frac{\partial U_i}{\partial x_j}\right\vert_{\boldsymbol{x}} \,{\rm d}^3{\boldsymbol{x}} . \end{equation}

Since the operator now involves the local interactions of the mean velocity gradient and since this condition must hold for all fields of $\partial U_i/\partial x_j$, the positive definiteness must be satisfied for each point. In other words, ${\mathsf{D}}^0_{ijkl} ({\partial U_j}/{\partial x_i} )({\partial U_l}/{\partial x_k}) \geqslant 0$ must also be satisfied pointwise for each local ${\mathsf{D}}^0_{ijkl}$ tensor and for all admissible values of mean velocity gradient. Therefore, the quadratic form of the eddy viscosity tensor must be non-negative implying that the eddy viscosity tensor must be positive semi-definite. This result is also similar to the condition considered by Milani (Reference Milani2020) in the context of local scalar transport.

Using our MFM measurement of the eddy viscosity, we can examine whether a local model from the truncated Kramer–Moyal expansion satisfies the positive semi-definite condition. Since the exact eddy viscosity must satisfy the positive semi-definite condition in (3.4), if this condition is not satisfied for the leading-order eddy viscosity, it is an indication that the local truncation is not valid and non-locality is needed for the positive definiteness condition.

To test the positive semi-definiteness of the leading-order eddy viscosity tensor, we flatten the eddy viscosity tensor and the velocity gradient. Then, the turbulent production becomes ${\mathsf{D}}^0_{ijkl} ({\partial U_j}/{\partial x_i})({\partial U_l}/{\partial x_k})=\boldsymbol {z}^T \boldsymbol{\mathsf{D}} \boldsymbol {z}$ where $\boldsymbol {z}$ is the flattened velocity gradient $[{\partial U_1} / {\partial x_1}\ {\partial U_2}/{\partial x_1}\ {\partial U_1}/{\partial x_2}\ {\partial U_2}/{\partial x_2}]^{\rm T}$ and where $\boldsymbol{\mathsf{D}}$ is the following matrix:

(3.6)\begin{equation} \boldsymbol{\mathsf{D}} =\begin{bmatrix} {\mathsf{D}}^0_{1111} & {\mathsf{D}}^0_{1112} & {\mathsf{D}}^0_{1121} & {\mathsf{D}}^0_{1122} \\ {\mathsf{D}}^0_{1211} & {\mathsf{D}}^0_{1212} & {\mathsf{D}}^0_{1221} & {\mathsf{D}}^0_{1222} \\ {\mathsf{D}}^0_{2111} & {\mathsf{D}}^0_{2112} & {\mathsf{D}}^0_{2121} & {\mathsf{D}}^0_{2122} \\ {\mathsf{D}}^0_{2211} & {\mathsf{D}}^0_{2212} & {\mathsf{D}}^0_{2221} & {\mathsf{D}}^0_{2222} \\ \end{bmatrix}. \end{equation}

It is well known that a symmetric $\boldsymbol{\mathsf{D}}$ is positive semi-definite if and only if all of its eigenvalues are non-negative However, for our case, $\boldsymbol{\mathsf{D}}$ is non-symmetric and $\boldsymbol {z}$ is limited to only certain value due to the incompressible condition. Therefore, we modified the quantity of interest. First, instead of the non-symmetric matrix $\boldsymbol{\mathsf{D}}$, we look at the positive definiteness of $\boldsymbol{\mathsf{D}}+\boldsymbol{\mathsf{D}}^T$. If $\boldsymbol{\mathsf{D}}+\boldsymbol{\mathsf{D}}^T$ is positive semi-definite, $\boldsymbol {z}^T \boldsymbol{\mathsf{D}} \boldsymbol {z} \geqslant 0$ also holds (Milani Reference Milani2020). Second, since the flow system is incompressible, $\boldsymbol {z}$ is limited to certain values satisfying ${\partial U_2} / {\partial x_2} = -{\partial U_1} / {\partial x_1}$. To expand the column vector multiplied to matrix $\boldsymbol{\mathsf{D}}$ to every non-zero real column vector, we must embed the incompressibility condition to the matrix $\boldsymbol{\mathsf{D}}$. We define $\boldsymbol {z} = \boldsymbol{\mathsf{C}} \boldsymbol {z^*}$ where $\boldsymbol {z^*}$ is the reduced flattened velocity gradient $[{\partial U_1} / {\partial x_1}\ {\partial U_2}/{\partial x_1}\ {\partial U_1}/{\partial x_2}]^{\rm T}$ and $\boldsymbol{\mathsf{C}}$ is the following matrix:

(3.7)\begin{equation} \boldsymbol{\mathsf{C}} =\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ -1 & 0 & 0 \end{bmatrix}. \end{equation}

Using this definition, $\boldsymbol {z}^T \boldsymbol{\mathsf{D}} \boldsymbol {z}$ becomes $\boldsymbol{z^{*T}} \boldsymbol{\mathsf{C}}^T\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{C}}\boldsymbol {z^*}$. Combining these two methods, we conclude that the eddy viscosity tensor is positive semi-definite when all the eigenvalues of the matrix $\boldsymbol{\mathsf{C}}(\boldsymbol{\mathsf{D}}+\boldsymbol{\mathsf{D}}^T)\boldsymbol{\mathsf{C}}^T$ is non-negative. We computed the smallest eigenvalue of this matrix at each $x_2$. The resulting plot is shown in figure 9. As shown, except for the thin zones near the wall, we see that the eigenvalues are positive, hence the eddy viscosity tensor is positive definite. Near the wall, however, the eigenvalues become negative, indicating that the leading term is not sufficient to capture the positive semi-definiteness of the eddy viscosity. In other words, this negativity occurs due to the local truncation of the eddy viscosity tensor and therefore non-locality should be incorporated to make the eddy viscosity positive semi-definite. The same analysis was conducted in the three dimension space context and the conclusion remains the same.

Figure 9. The minimum eigenvalue of the matrix $\boldsymbol{\mathsf{C}}(\boldsymbol{\mathsf{D}}+\boldsymbol{\mathsf{D}}^T)\boldsymbol{\mathsf{C}}^T$.

4.1. Non-locality

Figure 10 shows the full eddy viscosity kernel representation of the ${\mathsf{D}}_{2121}$ component. Each point in figure 10(a) represents the effect of the velocity gradient at the location $y_2$ to the Reynolds stress at the location $x_2$. The distribution of ${\mathsf{D}}_{2121}$ is confined to $x_2 \sim y_2$, indicating the approximate locality of this eddy viscosity component. At a given location $x_2$, we can visualise how much contribution the remote velocity gradient at a different location, $y_2$, makes to the Reynolds stress at the location $x_2$.

Figure 10. Distribution of ${\mathsf{D}}_{2121}$: (a) contour plot of ${\mathsf{D}}_{2121}$; (b) ${\mathsf{D}}_{2121}(x_2=x_2^*,y_2)$ with various $x_2^*$ where the blue line is at $x_2^*=0.502$ and the velocity gradient $\partial U_1 / \partial x_2$ is drawn as a red line; and (c) the blue solid line, $\int {\mathsf{D}}_{2121} \,\mathrm {d}y_2$, and the orange dashed line, ${\mathsf{D}}^0_{2121}$.

For instance, in figure 10(b), the thick blue line represents ${\mathsf{D}}_{2121}(x_2=0.502,y_2)$ and the distribution indicates the effects of the velocity gradient nearby. For a purely local eddy viscosity, a delta function around $y_2=0.502$ would b e expected. In the figure, even though the plotted profile is not a Dirac delta function, ${\mathsf{D}}_{2121}(x_2=0.502,y_2)$ shows concentrated behaviour around $y_2=0.502$. The blue line peak value is ${\mathsf{D}}_{2121}=0.53$ and the width that the curve drops to one third of its peak value is 0.12. Overall, our narrow banded results indicate that ${\mathsf{D}}_{2121}$ is relatively local throughout the domain. Such locality explains the earlier observation in § 3 where the leading-order eddy viscosity was shown to construct the shear component of the Reynolds stress within roughly 10 % error, and the mean velocity profile within 1 % error. However, a more quantitatively rigorous assessment would require examination of the mean velocity gradient across the kernel width. The leading-order eddy viscosity relegates the entire sensitivity of Reynolds stresses to the local pointwise value of the mean velocity gradient. If the mean velocity gradient happens to be relatively constant across the kernel width, the pointwise approximation, and hence the local model will be accurate. As shown in figure 10(b), the mean velocity gradient has a non-negligible variation across the sample kernel indicated by the blue curve. Therefore, strictly speaking, non-local effects in ${\mathsf{D}}_{2121}$ should not be negligible. We conclude that the accurate outcome of the local approximation for ${\mathsf{D}}_{2121}$ is partially owed to the error cancellation due to monotonic variation of the mean velocity gradient in the domain. In other words, the pointwise value of $\partial V_1/\partial x_2$ near the centroid of the kernel, reasonably represents the mean value given that errors from the left side and right side of the kernel partially cancel out each other.

Figure 10(c) compares the results of kernel integration $\int {{\mathsf{D}}_{2121} \,\mathrm {d}y_2}$ against the leading-order eddy viscosity for the same component ${\mathsf{D}}^0_{2121}$. The definition of ${\mathsf{D}}^0_{2121}$ is the leading-order moment of the eddy viscosity kernel ${\mathsf{D}}_{2121}$. In other word, with correct quantification the integration of the kernel $\int {{\mathsf{D}}_{2121} \,\mathrm {d}y_2}$ must match ${\mathsf{D}}^0_{2121}$. Figure 10(c) shows that the two results are collapsing verifying the consistency between our two different MFM measurements.

In § 3, we demonstrated that the leading-order eddy viscosity alone can predict an accurate RANS solution for the channel mean velocity with the prediction error around 1 %. This error can be further reduced by including the non-locality using the full kernel representation of the eddy viscosity. Figure 11 shows the two RANS results, one obtained using the leading-order eddy viscosity ${\mathsf{D}}^0_{2121}$ and the other obtained using eddy viscosity kernel ${\mathsf{D}}_{2121}$. Analytically, the full measurement of the kernel is expected to provide the RANS solution that is identical to the averaged DNS result. In our simulation, small errors are due to statistical noise that we expect to resolve with a larger data set. Still, the kernel result is significantly better than the leading-order result, indicating that the RANS solution to the leading-order eddy diffusivity model, which is highly local, can be improved using a non-local model.

Figure 11. Error in RANS prediction of the streamwise velocity using eddy viscosity tensor ${\mathsf{D}}^0_{2121}$ (denoted as MFM$^0$) and eddy viscosity kernel ${\mathsf{D}}_{2121}$ (denoted as MFM), $U_{1}-U_{1,{DNS}}$: dotted green line, $U_1$ predicted with ${\mathsf{D}}^0_{2121}$, and dashed orange line, $U_1$ predicted with ${\mathsf{D}}_{2121}$.

Next, we assess non-locality of other components of ${\mathsf{D}}_{ijkl}$ by examining the corresponding kernels. For example, the eddy viscosity kernel ${\mathsf{D}}_{1221}$ (figure 12) is widespread and shows significant non-locality, invalidating the intrinsic assumption in the Boussinesq approximation. The level of non-locality is drastic such that the velocity gradient at one half of the channel may affect the Reynolds stress at the other half of the domain. Moreover, the shape of ${\mathsf{D}}_{2121}$ differs from ${\mathsf{D}}_{1221}$, implying the non-universality of the kernel profile across different components of eddy diffusivity. Furthermore, the differences in the kernel shape between ${\mathsf{D}}_{2121}$ and ${\mathsf{D}}_{1221}$, clarifies why a truncated eddy viscosity operator based on the leading term of its Kramer–Moyal expansion can lead to asymmetric Reynolds stresses. Figure 13 shows the additional three non-zero eddy viscosity kernels ${\mathsf{D}}_{1121}$, ${\mathsf{D}}_{2221}$ and ${\mathsf{D}}_{3321}$. The rest of the components are zero due to channel symmetry. These three kernels correspond to the trace part of the eddy viscosity kernel and are also highly non-local.

Figure 12. Distribution of ${\mathsf{D}}_{1221}$: (a) contour plot of ${\mathsf{D}}_{1221}$; (b) ${\mathsf{D}}_{1221}(x_2=x_2^*,y_2)$ with various $x_2^*$ where the blue line is at $x_2^*=0.502$; and (c) blue solid line, $\int {\mathsf{D}}_{1221} \,\mathrm{d}y_2$, and orange dashed line, ${\mathsf{D}}^0_{1221}$.

Figure 13. Distribution of non-zero ${\mathsf{D}}_{ij21}$ where (a) ${\mathsf{D}}_{1121}$, (b) ${\mathsf{D}}_{2221}$ and (c) ${\mathsf{D}}_{3321}$.

4.2. Revisit of Reynolds stress

Lastly, we revisit the Reynolds stress reconstruction with the inclusion of the non-local effects. Figure 14 shows three different ways of constructing the Reynolds stresses. The first shown in the orange dashed line is the reconstructed Reynolds stress by the eddy viscosity kernel from MFM and the mean velocity gradient from DNS. The second shown in green dotted line is the reconstructed Reynolds stress by the leading-order eddy viscosity tensor from MFM and the mean velocity gradient from DNS. The last one is from the mean DNS data shown in blue solid line. The leading-order result and the mean DNS data are shown before in figure 8.

Figure 14. Reynolds stresses constructed by the leading-order eddy viscosity tensor and the eddy viscosity kernel associated with (a) ${\mathsf{D}}_{1121}$, (b) ${\mathsf{D}}_{1221}$, (c) ${\mathsf{D}}_{2121}$, (d) ${\mathsf{D}}_{2221}$ and (e) ${\mathsf{D}}_{3321}$: orange dashed line, the reconstructed Reynolds stress by the eddy viscosity kernel; green dotted line, the reconstructed Reynolds stress by the leading-order eddy viscosity tensor $\overline {u'_i u'_j}^0$; blue solid line, the DNS data.

Unlike the leading-order analysis, the results from the full kernel eddy viscosity match very well to the DNS data. These plots verify our computational method yielding two findings. First, with full kernels, the Reynolds stresses recover the symmetry that was lost in the leading-order approximation. Only ${\mathsf{D}}_{2121}$, which is relatively narrow banded, is applicable for the local approximation. Hence, this leads to the symmetry breakage after leading-order approximation. Second, now we can capture the non-zero Reynolds stress at the channel centreline. Thus, the measured non-local eddy viscosity allows prediction of non-zero Reynolds stresses near the centreline, whereas the leading-order approximation fails to do so.

4.3. Revisit of positive definiteness

In § 3.5, we discuss the positive definiteness of the local eddy viscosity tensor ${\mathsf{D}}^0_{ijkl}$. Due to the leading-order truncation of the eddy viscosity kernel, ${\mathsf{D}}_{ijkl}$, our result indicated that the local eddy viscosity tensor was not positive definite near the walls. In this section, we introduce the full kernel and examined whether including the non-locality restores the semi-positive definite condition.

Using the full eddy viscosity kernel expression and applying the fact that only one component of the velocity gradient tensor is non-zero, the turbulent production in (3.4) can be written as the following:

(4.1)\begin{equation} \int P_k \,\mathrm{d}\kern0.7pt x_2 = \int \int {\mathsf{D}}_{2121}({x_2}, {y_2})\frac{\partial U_1}{\partial y_2} \frac{\partial U_1}{\partial x_2} \,\mathrm{d}y_2 \,\mathrm{d}\kern0.7pt x_2 = \left[\frac{\partial U_1}{\partial x_2}\right]^{\rm T} \left[{\mathsf{D}}_{2121}\right] \left[\frac{\partial U_1}{\partial x_2}\right]. \end{equation}

The far right term represents the discrete form of the expression, where $[{\partial U_1}/{\partial x_2}]$ represents any velocity gradient vector at each point in $x_2$ and $[{\mathsf{D}}_{2121}]$ represents the discrete matrix value of ${\mathsf{D}}_{2121}({x_2}, {y_2})$. To make the turbulent production non-negative, the matrix $[{\mathsf{D}}_{2121}]$ needs to be semi-positive definite. Likewise in § 3.5, we computed eigenvalues of $[{\mathsf{D}}_{2121}] + [{\mathsf{D}}_{2121}]^{\rm T}$ to determine the positive definiteness. The computed eigenvalues range from 0.00 to 7.58, indicating that the eddy viscosity kernel ${\mathsf{D}}_{2121}$ is indeed semi-positive definite, recovering the stability condition that was lost by the leading-order truncation.