Direct calculation of the eddy viscosity operator in turbulent channel flow at Re_τ = 180

Abstract

This study aims to quantify how turbulence in a channel flow mixes momentum in the mean sense. We applied the macroscopic forcing method (Mani & Park, Phys. Rev. Fluids, 2021, 054607) to direct numerical simulation (DNS) of a turbulent channel flow at $Re_\tau =180$ using two different forcing strategies that are designed to separately assess the anisotropy and non-locality of momentum mixing. In the first strategy, the leading term of the Kramers–Moyal expansion of the eddy viscosity is quantified, revealing all 81 tensorial coefficients that essentially characterise the local-limit eddy viscosity. The results indicate the following: (1) the eddy viscosity has significant anisotropy, (2) Reynolds stresses are generated by both the mean strain rate and mean rotation rate tensors associated with the momentum field and (3) the local-limit eddy viscosity generates asymmetric Reynolds stress tensors. In the second strategy, the eddy viscosity is quantified as an integration kernel revealing the non-local influence of the mean momentum gradient at each wall-normal coordinate on all nine components of the Reynolds stresses over the channel width. Our results indicate that while the shear component of the Reynolds stress is reasonably reproduced by the local mean gradients, other components of the Reynolds stress are highly non-local. These results provide an understanding of anisotropy and non-locality requirements for closure modelling of momentum transport in attached wall-bounded turbulent flows.

1. Introduction

Many of the turbulence models in use today are based on the Boussinesq approximation (Boussinesq 1877) in which the Reynolds stresses are assumed to be a linear function of the local mean velocity gradients. This approximation furthermore assumes isotropy of the tensor representing the coefficients of this linear relation, which is commonly referred to as eddy viscosity. The two simplifications offered by the Boussinesq approximation reduce the job of turbulence modelling to a determination of a scalar eddy viscosity field from which local Reynolds stresses can be determined algebraically without the need to solve any additional equations. For cases in which a single component of the Reynolds stress plays the dominant role, such as in parallel flows, a scalar eddy viscosity can be tuned to yield acceptable Reynolds stress fields (Pope 2000). However, most turbulence models utilise this approximation even for multi-dimensional flows (Hanjalić & Launder 1972; Chien 1982; Durbin 1993; Menter 1994; Spalart & Allmaras 1994; Wilcox 2008). While some models allow anisotropic eddy viscosities (Spalart 2000; Mani et al. 2013; Rumsey et al. 2020), they still retain the locality of the Reynolds stress dependence on the mean velocity gradient.

Experimental measurements, as well as direct numerical simulation (DNS) data, suggest that the isotropy and locality assumptions of the Boussinesq approximation are not strictly valid. Several studies have shown significant misalignment between the principal axis of the Reynolds stress and strain rate tensors indicating non-negligible anisotropy of the eddy viscosity operator (Champagne, Harris & Corrsin 1970; Harris, Graham & Corrsin 1977; Rogallo 1981; Moin & Kim 1982; Rogers & Moin 1987; Coleman, Kim & Le 1996). Furthermore, the assumption of Reynolds stress locality is often not true because turbulent mixing may exist from the history of the straining in a given region of a turbulent flow. For instance, the experiment conducted by Warhaft (1980) showed that the Reynolds stress can arise from the history effects of straining, even with a locally zero mean strain rate. In this case, the Reynolds stress should incorporate temporal or spatial non-locality of the strain rate tensor.

Given these pieces of evidence, various modelling techniques have attempted to relax both locality and isotropy assumptions via development of second-order closure models (Launder, Reece & Rodi 1975; Wilcox 1998; Speziale, Sarkar & Gatski 1991; Gerolymos et al. 2012; Cécora et al. 2015) often using the Reynolds stress transport equation as a framework to identify the needed closures. Each of these models, provides a specific way in which Reynolds stresses could depend non-locally or anisotropically on the velocity gradient field. However, standard data of turbulent flows, either from DNS or experiments, do not provide sufficient information to allow proper discrimination between these models. While these data reveal anisotropy of the Reynolds stresses, they do not uniquely determine the anisotropy or non-locality of the closure operators that express their dependence on the mean velocity gradient. Closing this gap would require quantification of the eddy viscosity as an operator acting on the mean velocity gradient. With this goal in mind, this study presents a direct quantification of the eddy viscosity operator in a canonical turbulent flow via utilisation of the macroscopic forcing method (MFM), developed by Mani & Park (2021).

Prior to the description of our work, we start by reviewing generalised forms of the eddy diffusivity and eddy viscosity operators for scalar and momentum transport in turbulent flows. First, one way of generalising the Boussinesq approximation is to allow for the anisotropy of the eddy viscosity. Batchelor (1949) suggested using a second-order tensor replacing the diffusion coefficient in the Fickian model to describe the mean transport of a scalar quantity. Later, a similar concept was suggested by Rogers, Mansour & Reynolds (1989), where the mean turbulent flux of a passive scalar was approximated with an algebraic model expressed in a second-order tensor eddy diffusivity. This anisotropic eddy diffusivity model can be written as the following: $-\overline {u'_i c'}={\mathsf{D}}_{ij} \partial C/\partial x_j$ where $\overline {({\cdot })}$ represents ensemble-average, $u'_i$ represents the fluctuation of the velocity, $C$ and $c'$ represent the mean and the fluctuation of the scalar quantity being transported, $x_i$ represents the spatial Cartesian coordinate and ${\mathsf{D}}^0_{ij}$ represents the second-order eddy diffusivity tensor that is local.

Similarly, for the turbulent momentum flux, one method of generalising the Boussinesq approximation is to use a tensorial representation of the eddy viscosity. Hinze (1959) has suggested the use of the fourth-order tensor as the eddy viscosity. Later, Stanišić & Groves (1965) conducted a systematic investigation of the tensorial character of the eddy viscosity coefficient and revealed that the eddy viscosity tensor has to be at minimum fourth order. In parallel to the anisotropic eddy diffusivity model, the anisotropic eddy viscosity model for momentum transport can be written as $-\overline {u'_i u'_j} = {\mathsf{D}}^0_{ijkl} \partial {U_l}/\partial {x_k}$, where $U_l$ represents mean velocity field. Here, the Reynolds stress $\overline {u'_i u'_j}$ is locally closed in terms of the fourth-order tensorial eddy viscosity ${\mathsf{D}}^0_{ijkl}$ and the mean velocity gradient.

An even more general form of the eddy viscosity can be used to incorporate not only anisotropy but also non-locality. Hamba (2005, 2013) suggested writing the closure of the Reynolds stress in terms of the mean velocity gradient at remote times and locations. This form of eddy viscosity involves a fourth-order tensorial kernel, which we refer to as the eddy viscosity kernel. For statistically stationary flows, this relation can be expressed as

(1.1)\begin{equation} -\overline{u_i'u_j'} (\boldsymbol{x})= \int {\mathsf{D}}_{ijkl}(\boldsymbol{x}, \boldsymbol{y}) \left. \frac{\partial U_l}{\partial x_k} \right\vert_{\boldsymbol{y}} \,{\rm d}^{3}{\boldsymbol{y}}, \end{equation}

where ${\mathsf{D}}_{ijkl}(\boldsymbol {x}, \boldsymbol {y})$ is the eddy viscosity kernel indicating how mean gradients at location $\boldsymbol {y}$ result in Reynolds stresses at location $\boldsymbol {x}$. When written in dimensional form, the eddy viscosity kernel does not have the same dimension as the kinematic viscosity. Instead, the kernel represents increment of viscosity per unit volume of the non-locality dimension. In the case of (1.1) temporal non-locality is not considered due to the system's statistical stationarity, but full three-dimensional spatial non-locality is considered. As a result, the dimensional kernel would have unit of diffusivity per unit volume or ${\rm m}^{-1}\ {\rm s}^{-1}$.

Hamba (2005) reported the first quantification of the eddy viscosity kernel for a turbulent channel flow using a Green's function formulation approach based on an earlier work by Kraichnan (1987). However, their study focuses on a subset of the tensorial coefficients, i.e. ${\mathsf{D}}_{ij21}$. This choice is motivated since the mean velocity profile in the channel flow is insensitive to other components of the eddy viscosity kernel, given the mean velocity gradient, $\partial U_l / \partial x_k$ shown in (1.1), is non-zero only for $(k,l)=(2,1)$. Nevertheless, quantification of other components of eddy viscosity in this canonical setting would provide significant insights about momentum mixing in the broader context of wall-bounded shear flows. Aside from this shortcoming, Hamba (2005) chose to manually enforce the symmetry ${\mathsf{D}}_{ijkl}={\mathsf{D}}_{jikl}$ by performing arithmetic averaging of the respective components (i.e. $ij$ and $ji$) of the output data from their simulations. This choice was made given the expectation that the Reynolds stress tensor as the output of (1.1) must always be symmetric, while the raw kernels did not follow this symmetry.

Recently, Mani & Park (2021) presented an alternative interpretation of (1.1) in the context of the generalised momentum transport (GMT) equation. GMT can be derived by applying the Reynolds transport theorem to momentum transport without constraining the momentum field to be identical to the velocity field. In this context, the Reynolds stress, expressed as $\overline {u_i'v_j'}$, is interpreted as the mean product of two conceptually different fields, with $u_i$ representing the kinematic displacement of volume acting as a transporter of momentum, and $v_j$ representing momentum per unit mass, the quantity of interest that results in friction and pressure. Navier–Stokes (NS) is rendered as a special solution to GMT in which the two fields are constrained to be equal. Specifically, when GMT is supplied with the same boundary conditions and forcing conditions as those in NS, the solution to NS is the only attractor solution to GMT, as shown theoretically and numerically by Mani & Park (2021). With this interpretation, (1.1) is in fact a closure operator to the ensemble-averaged GMT and not the Reynolds-averaged Navier–Stokes (RANS) equation. Therefore, ${\mathsf{D}}_{ijkl}$ and ${\mathsf{D}}_{jikl}$ are not required to be equal, since $\overline {u_i'v_j'}\ne \overline {u_j'v_i'}$. The present study addresses this issue, by examining the raw eddy viscosity operator without any symmetry averaging. We confirm that while the eddy viscosity kernel of channel flow is not symmetric, it still results in symmetric Reynolds stresses when it acts on the mean velocity gradient of the same flow from which the eddy viscosity data are obtained.

As mentioned previously, Mani & Park (2021) provide a statistical technique called the MFM, which allows direct measurement of a flow's eddy viscosity ${\mathsf{D}}_{ijkl}$ with data gathered from DNS of the NS equation and GMT. More generally speaking, MFM allows precise computation of RANS closure operators via applying various macroscopic forcing to the GMT equations which can be utilised to extract the eddy viscosity operator. It is worth noting that macroscopic forcing is not limited to delta functions, which reveal Green's functions as outputs. For instance, Shirian & Mani (2022) employed harmonic forcing to efficiently unveil the eddy diffusivity operator for homogeneous isotropic turbulence. They successfully fitted this operator with an analytical expression. An alternative approach by Mani & Park (2021), which is more relevant to this study, is the inverse macroscopic forcing method (IMFM), in which forcing to constrain mean polynomial fields was shown to reveal non-local moments of the underlying eddy diffusivity operator in an economical way compared with the Green's function approach. We examine a systematic procedure for obtaining a local operator approximation of the full eddy viscosity operator by considering a Kramers–Moyal expansion (Van Kampen 1992) of the eddy viscosity operator and quantifying its leading term. This approach does not only enable estimation of the eddy viscosity in an economical fashion, but it also separates out the easy-to-comprehend local eddy viscosity by utilising this established expansion, which we believe was a missing piece in the analysis of Hamba (2005).

The rest of this paper is organised as follows. In § 2, we define the flow system and the model used which involves the fourth-order tensorial eddy viscosity kernel, and review the computational methodology. In § 3, we begin with evaluating the isotropy assumption in the Boussinesq's approximation. For simplicity, we conduct the leading-order (local-limit) approximation to the eddy viscosity kernel to solely focus on the anisotropy of the eddy viscosity. With the measured local eddy viscosity tensor, we discuss the following: the standard eddy viscosity, the quantified anisotropy, the dependency of Reynolds stress on the rate of rotation, the leading-order Reynolds stress and the positive definiteness of the leading-order eddy viscosity operator. In § 4, we extend our study to non-local effects, by computing the full eddy viscosity kernel representing the non-local effects in the wall-normal direction. In § 5, we summarise our results and discuss potential extensions to this study.

2. Problem set-up and governing equations

Figure 1 shows the schematics of the channel flow and its coordinate system where the flow is bounded by top and bottom walls spaced $2\delta$ apart. We denote the Cartesian coordinates $x_i$, where $x_1$ is the streamwise direction, $x_2$ is the wall-normal direction, and $x_3$ is the spanwise direction. The dimensionless equations expressing mass and momentum conservation are as follows:

(2.1)$$\begin{gather} \frac{\partial u_i}{\partial t} + \frac{\partial u_j u_i}{\partial x_j} ={-} \frac{\partial p}{\partial x_i}+\frac{1}{{Re}} \frac{\partial^2 u_i}{\partial x_j \partial x_j} + r_i , \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial u_j}{\partial x_j} = 0 , \end{gather}$$

where $u_i$ is the flow velocity, $p$ is the pressure normalised by the density, $t$ is time, and $r_i=(1,0,0)$ represents normalised mean pressure gradient. The dimensionless spatial coordinates are normalised by $\delta$ and ${Re}$ represents the Reynolds number defined based on $\delta$ and the friction velocity $u_\tau = \sqrt {\tau _w / \rho }$ where $\tau _w$ is the mean wall shear stress balancing the force due to the mean pressure gradient and $\rho$ is the fluid density.

Figure 1. Schematics of the channel flow.

The RANS equations can be obtained by taking the ensemble-average of (2.1) and (2.2), yielding

(2.3)$$\begin{gather} \frac{\partial U_i}{\partial t} + \frac{\partial U_j U_i}{\partial x_j} ={-} \frac{\partial \bar{p}}{\partial x_i} + \frac{1}{{Re}} \frac{\partial^2 U_i}{\partial x_j \partial x_j} - \frac{\partial \overline{u'_j u'_i}}{\partial x_j} + \bar{r}_i, \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial U_j}{\partial x_j}=0, \end{gather}$$

where $U_{i}$ is the mean velocity, $u'_{i}$ is the velocity fluctuation around the mean velocity and $\overline {({\cdot })}$ implies ensemble-averaged quantities. To close this system, the divergence of Reynolds stresses, ${\partial \overline {u'_j u'_i}}/{\partial x_j}$, needs to be modelled in terms of the primary variable $U_i$. This can be generally expressed as an operator acting on the ensemble-averaged field, $- {\partial \overline {u'_j u'_i}}/{\partial x_j} \equiv \bar {\mathcal {L}}( U_i )$. One form of such operators is expressed in (1.1).

A DNS solution to the channel flow does not provide enough information to fully quantify the non-local eddy viscosity kernel, ${\mathsf{D}}_{ijkl}(\boldsymbol {x}, \boldsymbol {y})$. A full characterisation of $\boldsymbol{\mathsf{D}}$ requires quantification of Reynolds stresses in response to all possible independent flow gradients scenarios. Following Mani & Park (2021), we next describe the procedure of obtaining $\boldsymbol{\mathsf{D}}$. In this paper, however, we limit the scope of our analysis to the one-dimensional RANS context, in which the wall-normal coordinate is the only independent variable since the flow is statistically homogeneous in all other space–time coordinates. In other words, we assume the form ${\mathsf{D}}_{ijkl} = {\mathsf{D}}_{ijkl} (x_2, y_2)$, and the other dimensions are integrated out in (1.1). However, the employed macroscopic forcing methodology is in principle generalisable to multi-dimensional cases, and with higher computational expense can capture the full behaviour of ${\mathsf{D}}_{ijkl} (\boldsymbol {x}, \boldsymbol {y})$.

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 2021). 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 (2021), 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 2021). 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 (2005) 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 (2005) 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 (2023) that such a brute force approach is theoretically equivalent to Hamba's Green's function approach (Hamba 2005).

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 1992), as investigated by Mani & Park (2021). 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 (2021), 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 (1949), 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 (2005), 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 (2021), 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 (2010) and modified by Seo, García-Mayoral & Mani (2015). The present DNS uses the fractional step method with semi-implicit time advancement (Kim & Moin 1985). 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. 1998). 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. Anisotropy analysis

In this section, we compute the leading-order eddy viscosity tensor ${\mathsf{D}}^0_{ijkl}$ and focus on the analysis on anisotropy of the eddy viscosity and specifically contrast it to the standard eddy viscosity implied by the Boussinesq model. In addition, we assess dependency of Reynolds stresses on the rate of rotation, examine reconstruction of Reynolds stresses using the leading-order eddy viscosity, and lastly discuss positive definiteness of the leading-order eddy viscosity tensor.

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 (2021) 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 (2005) 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 2005). 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. 2015), 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 (2022) 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 2000; Mani et al. 2013; Rumsey et al. 2020). 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 1992) and (2) corrections to the Boussinesq eddy viscosity based on Cayley–Hamilton theorem (Pope 1975), which are also incorporated in quadratic constitutive relation (QCR) models (Spalart & Allmaras 1994). 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 (2021), 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 1977; Schumann 1977). We often require eddy viscosity to satisfy the same condition as done in the Boussinesq approximation with $\nu _T \geqslant 0$ (Speziale, Abid & Durbin 1994) 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. 2019) 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 (2020) 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 2020). 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. Non-locality analysis

In § 3.3, we concluded that aside from $\overline {u'_2 u'_1}$, capturing other components of the Reynolds stress field requires inclusion of non-local terms in the eddy viscosity operator. To better understand the non-local effect, in this section, we investigate the full eddy viscosity kernel ${\mathsf{D}}_{ijkl}(x_2, y_2 )$ in (2.7).

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.

5. Conclusions

This study presents a quantification of non-Boussinesq effects in eddy viscosity in a subclass of turbulent wall-bounded flows. The presented analysis is systematically focused on two aspects: anisotropy and non-locality of momentum mixing. To assess these effects and quantify the deviation from Boussinesq limit, we calculated the eddy viscosity of the turbulent channel flow at ${Re}_\tau =180$ using a statistical technique that we recently developed called MFM. Using MFM, we quantified the leading-order eddy viscosity tensor for the analysis of anisotropy and expanded our study to quantify the eddy viscosity tensorial kernel for the analysis of non-locality.

Our results indicate the following: (1) eddy viscosity is highly anisotropic with some elements orders of magnitude larger than the nominal eddy viscosity; (2) the Reynolds stresses reconstructed from this eddy viscosity depends not only on the mean rate of strain but also on mean rate of rotation; (3) leading-order eddy viscosity, which is obtained by neglecting higher spatial moments of the closure kernel, generates a non-symmetric Reynolds stress tensor; and (4) aside from the shear component of the Reynolds stress, $\overline {u^\prime _2 u^\prime _1}$, which showed a limited level of non-locality, the dependence of other components of Reynolds stress on mean velocity gradient is highly non-local at the level where some components of the Reynolds stress are influenced by the velocity gradient on the other half of the channel.

The exact measurement of the eddy viscosity of the channel flow has different implication for RANS modelling of parallel flows and that of the spatially developing attached boundary layers. For parallel flows, only one Reynolds stress component and one velocity gradient are important; hence anisotropy does not influence the prediction of the mean flow as long as ${\mathsf{D}}_{2121}$ is properly modelled. At the same time, not only the anisotropy but also non-locality may be omitted for the channel flow. This outcome is in part due to relatively narrower kernel of ${\mathsf{D}}_{2121}(x_2, y_2)$, as shown in our MFM measurement of the eddy viscosity kernel, but also due to coincidental error cancellations that render reasonable estimation of the Reynolds stresses based on a single-point quadrature relegating the entire weight of the kernel on the local mean velocity gradient.

These two findings may explain why the Boussinesq approximation works well for prediction of mean parallel flows. However, our quantification suggests that this conclusion does not hold for normal components of the Reynolds stress, as well as for spatially developing wall-bounded flows where the non-parallel effects become important. For instance, even a small gradient in the streamwise direction can have a non-negligible effect since ${\mathsf{D}}_{1111}^0$ is very large compared with most of other eddy viscosity components. Our measurements reveal that the eddy viscosity is highly anisotropic and highly non-local, when it comes to components other than ${\mathsf{D}}_{2121}$, indicating a clear need to include non-Boussinesq effects in RANS models.

While we focused on full non-local analysis in the $x_2$ direction, we did not consider non-local spatial effects in other directions and non-local temporal effects. Equation (2.7) is a reduced version of (1.1) using leading-order moments in $x_1$, $x_3$ and $t$. These leading-order reductions are justified for channel flow since it is statistically homogeneous in these directions, and are expected to be qualitatively valid for systems with slow variation of turbulence in these directions. While it is possible to quantitatively assess such effects with MFM, we defer analysis of streamwise and spanwise non-locality in eddy viscosity to a future study.