Macroscopic forcing method: a tool for turbulence modeling and analysis of closures

TL;DR

The paper introduces the macroscopic forcing method (MFM), a numerical technique to identify differential operators governing turbulence closures, enabling precise analysis of eddy diffusivity and momentum transport in complex flows.

Contribution

This work develops MFM as a novel, quantitative tool for revealing turbulence closure operators, including nonlocal and anisotropic effects, across various flow configurations.

Findings

MFM accurately determines eddy diffusivity operators in canonical flows.

The method captures nonlocal mixing effects through Fourier space operator representation.

Extension of MFM to wall-bounded flows and momentum transport provides exact RANS solutions.

Abstract

This study presents a numerical procedure, which we call the macroscopic forcing method (MFM), which reveals the differential operators acting on the mean fields of quantities transported by underlying fluctuating flows. Specifically, MFM can precisely determine the eddy diffusivity operator, or more broadly said, it can reveal differential operators associated with turbulence closure for scalar and momentum transport. We present this methodology by considering canonical problems with increasing complexity. Starting from the well-known problem of dispersion of passive scalars by parallel flows we elucidate the basic steps in quantitative determination of the eddy viscosity operators using MFM. Utilizing the operator representation in Fourier space, we obtain a stand-alone compact analytical operator that can accurately capture the nonlocal mixing effects. Furthermore, a cost-effective generalization of MFM for analysis of nonhomogeneous and wall-bounded flows is developed and is comprehensively discussed through a demonstrative example. Extension of MFM for analysis of momentum transport is theoretically constructed through the introduction of a generalized momentum transport equation. We show that closure operators obtained through MFM analysis of this equation provide the exact RANS solutions obtained through averaging of the Navier-Stokes equation. We introduce MFM as an effective tool for quantitative understanding of non-Boussinesq effects and assessment of model forms in turbulence closures, particularly, the effects associated with anisotropy and nonlocality of macroscopic mixing.

Equations

$$\mathcal{L}c\left(x_1,...,x_n,t\right)=0,$$

$$\mathcal{L}=\frac{\partial }{\partial t}+u_j\frac{\partial }{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_M\frac{\partial }{\partial x_j}\right),$$

$$\overline{\mathcal{L}}\overline{c}=0$$

$$\frac{\partial \overline{c}}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{u_{j}c}\right)-\frac{\partial }{\partial x_j}\left(D_M\frac{\partial \overline{c}}{\partial x_j}\right)=0.$$

$$\overline{\mathcal{L}}=\frac{\partial }{\partial t}+\overline{u_j}\frac{\partial }{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_M\frac{\partial }{\partial x_j}\right)+\overline{{\mathcal{L}}^{\prime }}.$$

$$u_1=U\cos \left(\frac{2\pi }{L_2}{x}_2\right),u_2=0,$$

$$\frac{\partial c}{\partial t}+u_j\frac{\partial c}{\partial x_j}=D_M\frac{{\partial }^{2}c}{\partial x_1^2}+D_M\frac{{\partial }^{2}c}{\partial x_2^2}.$$

$$\overline{c}(x_1)=\frac{1}{L_2}{\int }_0^{L_2}c\left(x_1,x_2\right)\text{d}{x}_2.$$

$$\frac{\partial c}{\partial t}+\cos \left(x_2\right)\frac{\partial c}{\partial x_1}=\frac{{\partial }^{2}c}{\partial x_2^2}+{ϵ}^2\frac{{\partial }^{2}c}{\partial x_1^2},$$

$$\frac{\partial \overline{c}}{\partial t}+\frac{\partial }{\partial x_1}\left(\overline{\cos (x_2)c}\right)=0.$$

$$\frac{{\partial }^2{c}^{\prime }}{\partial x_2^2}-\cos \left(x_2\right)\frac{\partial \overline{c}}{\partial x_1}=\frac{\partial c^{\prime }}{\partial t}+\frac{\partial }{\partial x_1}{\left[\cos \left(x_2\right)c^{\prime }\right]}^{\prime },$$

$$\frac{\partial \overline{c}}{\partial t}=\frac{1}{2}\frac{{\partial }^2\overline{c}}{\partial x_1^2}.$$

$$\frac{\partial }{\partial t}\left(\overline{c}+\frac{1}{2}\frac{{\partial }^2\overline{c}}{\partial x_1^2}\right)=\frac{1}{2}\frac{{\partial }^2\overline{c}}{\partial x_1^2}+\frac{1}{32}\frac{{\partial }^4\overline{c}}{\partial x_1^4}.$$

$$[\mathcal{L}][c]=0.$$

$$\left[\overline{\mathcal{L}}\right]\left[\overline{c}\right]=0,$$

$$\overline{\mathcal{L}}\overline{c}\left(s\right)=s,$$

$$\mathcal{L}\left[c\left(x_1,...,x_n,t\right)\right]=s.$$

$$\frac{\partial c}{\partial t}+\cos \left(x_2\right)\frac{\partial c}{\partial x_1}=\frac{{\partial }^{2}c}{\partial x_2^2}+s(x_1,t).$$

$$s(x_1,t)=\exp \left(i\omega t+ik{x}_1\right).$$

$$\overline{c}\left(x_1,t\right)=\overline{\hat{c}}\exp \left(i\omega t+ik{x}_1\right).$$

$$\widehat{\overline{\mathcal{L}}}=1/\overline{\hat{c}}\left(\omega ,k\right).$$

$$c=\hat{c}(\omega ,k,x_2)\exp \left(i\omega t+ik{x}_1\right).$$

$$\left[i\omega -\frac{{\partial }^2}{\partial x_2^2}+ik\cos \left(x_2\right)\right]\hat{c}\left(\omega ,k;x_2\right)=1,$$

$$\hat{c}\left(\omega ,k;x_2\right)=\sum _{n=0}^{\infty }{a}_n\left(\omega ,k\right)\cos (n{x}_2).$$

$$\hat{\overline{{\mathcal{L}}^{\prime }}}=1/\overline{\hat{c}}\left(\omega ,k\right)-i\omega$$

$$\hat{\overline{{\mathcal{L}}^{\prime }}}\simeq \sqrt{{\left(1+i\omega \right)}^2+k^2}-\left(1+i\omega \right).$$

$$\left[\sqrt{{\left(\mathcal{I}+\frac{\partial }{\partial t}\right)}^2-\frac{{\partial }^2}{\partial x_1^2}}-\mathcal{I}\right]\overline{c}\left(x_1,t\right)=s\left(x_1,t\right),$$

$$\hat{\overline{{\mathcal{L}}^{\prime }}}=\frac{0.5k^2}{\sqrt{1+0.25k^2}}.$$

$$\overline{{\mathcal{L}}^{\prime }}=-\frac{\partial }{\partial x_1}\left[\left(\frac{0.5}{\sqrt{\mathcal{I}-\frac{1}{4}\frac{{\partial }^2}{\partial x_1^2}}}\right)\frac{\partial }{\partial x_1}\right],$$

$$\mathcal{D}=\frac{D}{\sqrt{\mathcal{I}-l^2\frac{{\partial }^2}{\partial x_1^2}}},$$

Full text

I Introduction

The kinetic theory of gases has been an invaluable tool for developing continuum models governing transport of mass, momentum, and energy. Without such models, prediction of flows and heat transfer in natural and industrial systems would have required tracking of individual fluid molecules by solving an enormous system of coupled molecular dynamics equations. To date, the largest of such calculations can only span domains of size of order $100nm$ which is by far much smaller than many engineering scales. Kinetic theory bridged this gap by analyzing interactions of molecules in a statistical sense. These analyses provide transport coefficients to partial differential equations (PDEs) describing transport phenomena in a macroscopic senseVincenti1965 . Unlike the underlying chaotic Brownian dynamics, these scaled-up Eulerian PDEs result in smooth solutions with characteristic length scales multiple orders of magnitude larger than molecular scales. These PDEs therefore offer drastic computational savings by lowering the required number of degrees of freedom in the systems to be solved.

The diffusion equation and the Navier-Stokes equation governing the evolution of mass, heat, and momentum are examples of such transport PDEs. These equations however, do not completely close the gap of length scales all the way to the engineering scales. Most notably, as many engineering applications involve flows at high Reynolds numbers, the Navier-Stokes equation itself admits chaotic solutions spanning a wide spectrum of length scalesMoin1997 , rendering the system-level computations prohibitively expensive. Again, as primary quantities of interest involve only the averaged fields, researchers have actively sought approaches for further scale up of the mathematical models over the past decades.

Central in this context is the RANS (Reynolds-Averaged Navier-Stokes) methodology Pope2001 ; Speziale1991 ; Speziale1998 ; Durbin2018 . The goal of RANS models is to predict the average velocity field in a turbulent flow without the need of tracking individual chaotic eddies. In this case, the average is the mean over many ensembles of realizations of the same system under the same boundary conditions. However, as most systems involve statistically time-stationary turbulence, temporal averages are equivalent to ensemble averagesPope2001 . Both averages result in highly smooth mean fields with a much smaller number of degrees of freedom suitable for pursuing scaled-up models.

RANS modeling has, from its inception, been inspired by approaches from kinetic theory. Most notably, the majority of the RANS modelsSpalart1992 ; Hanjalic1972 ; Wilcox2008 ; Menter1994 ; Durbin1995 are based on the Boussinesq approximationBoussinesq1877 suggesting closures in terms of diffusion operators with the coefficient referred to as eddy diffusivity. By making an analogy to the kinetic theory, this approximation implies that in the same way that mixing by Brownian motion leads to a scaled-up continuum model involving a diffusion equation with an isotropic coefficient tensor, turbulent mixing should be expressible in terms of diffusion operators when statistically averaged. However, unlike such Brownian motion, turbulent fluctuations are often non-isotropic in practice. Aside from this issue, a major criticism of this analogy is that the diffusion approximation requires separation of scales between the scaled-up fields and the underlying microscopic, sub-scale, physics. This is the case for molecular-to-continuum modeling since the continuum scale is often orders of magnitude larger than the mean free path of molecules. However, in turbulent flows the underlying chaotic eddies span a wide range of scales; large eddies, which have the highest statistical significance, are often as large as the scale of the mean fields. Therefore, the Boussinesq approximation, while being useful for qualitative analysis and scaling assessments, is not accurate for quantitative predictions.

Despite these shortcomings, the majority of turbulence models in use employ eddy-diffusivity closures based on the Boussinesq approximationSpalart1992 ; Hanjalic1972 ; Wilcox2008 ; Menter1994 ; Durbin1995 . Among the models that do not use the Boussinesq approximation, some still retain diffusion-type operators in which the mean momentum flux is expressed as a function of local velocity gradientsPope1975 . Even the models that do not use any closure at the momentum transport level, such as the Reynolds stress closure modelsLaunder1975 , still retain similar approximations in closure of higher moments. Model-form inaccuracies in turbulence models are to some degree compensated by tuning of the model coefficients and tailoring them for specific regimes and applications. As a result, turbulence models are still far from being universal and truly predictive.

In this report we present a method, which we call the macroscopic forcing method (MFM) that allows determination of the “closure operators” that govern the mean-field-mixing by any underlying flow field. MFM acts similar to the way that molecular dynamics simulations reveal transport coefficients for a continuum model, with the exception that in MFM both input and output spaces are continua and that MFM does not make any simplifying assumption such as isotropy or separation of scales between turbulence and mean fields (scaled up fields). MFM precisely reveals the degree of non-locality and anisotropy of the differential operators governing the mean transport.

Lastly, the approach, which we shall introduce, unifies the problem of macroscopic modeling between low-Reynolds-number laminar flows and high-Reynolds-number turbulent flows. As a starting point, we consider the problem of dispersion of passive scalars by a laminar parallel flow, first introduced by G. I. TaylorTaylor1953 . By recasting this problem in wavenumber space, we show that while MFM correctly reproduces Taylor’s approximation for the low wavenumber limit, it provides a natural framework for extending his solution to all wavenumbers. The paper then continues to pave the way for the development of MFM for more complex flows, involving statistical inhomogeneity and wall boundary conditions. Lastly, an extension of MFM to the analysis of momentum transport is theoretically formalized by introducing a microscopic equation representing momentum transport in the generalized sense, and providing a theoretical proof regarding convergence of the statistical solution of this equation to the solution obtained from Reynolds-averaging of the Navier-Stokes equation. The paper ends by providing useful insights about MFM from both physical and mathematical point of views.

I.1 Problem Statement

The purpose of this section is to introduce the basic terminology and mathematical notation used throughout the paper. While we ultimately develop a general approach that can treat a wide range of transport phenomena, in this section we introduce the problem by considering the passive scalar advection-diffusion equation. Later, we introduce extensions to momentum fields governed by the Navier-Stokes equation. The starting point is the continuum microscopic equation:

[TABLE]

where $c$ represents the transported field, while $x_{1}$ to $x_{n}$ and $t$ are the independent variables, which correspond to the spatial and temporal coordinates, respectively. $\mathcal{L}$ is a differential operator representing the transport physics. For example, when the advective flow is incompressible, one can write

[TABLE]

where $u_{j}\left(x_{1},...,x_{n},t\right)$ represents the advecting velocity field, and $D_{M}$ represents molecular diffusion. Equation (1) is accompanied by initial and boundary conditions. When (1) is written in discretized form, these conditions can be embedded in the operator $\mathcal{L}$ and the right hand side of (1).

We refer to Equation (1) as the “microscopic equation” or the “microscopic transport equation,” where $\mathcal{L}$ is called the microscopic operator acting on the microscopic field $c$. Next, we define the “macroscopic” or scaled-up field, $\overline{c}$, which represents the quantity of engineering interest. Mathematically, $\overline{c}$ is defined as the average of $c$ over a certain subset of spatiotemporal coordinates. For example, in dispersion of scalars in laminar pipe flowsTaylor1953 , averaging is performed in the pipe cross-section. In RANS modeling, the averaging dimensions are those over which the system is statistically homogeneous. In most cases, the temporal dimension is statistically stationary. In the absence of statistical stationarity, averaging is defined over many ensembles. A more advanced definition of $\overline{c}$, but outside of the scope of this report, utilizes filters or weighted averages.

The problem is to determine the“macroscopic” operator, $\overline{\mathcal{L}}$, with an appropriate set of boundary conditions such that the solution to equation

[TABLE]

matches exactly the $\overline{c}$ obtained from averaging the microscopic solution. Since the macroscopic space, has a lower dimension than the microscopic space, or better said, involves a significantly lower number of degrees of freedom, access to the macroscopic operator allows tremendous savings in computational cost. As we shall see, $\overline{\mathcal{L}}$ will depend on the statistics of the microscopic advective field $u_{j}$, $D_{M}$, the microscopic boundary conditions, and also the definition of averaging that converts $c$ to $\overline{c}$.

The immediate difficulty with the above problem statement is that there are numerous macroscopic operators that satisfy the stated requirements. In other words, the answer to the problem is not unique. In what follows we will utilize additional physics-based constraints such that $\overline{\mathcal{L}}$ will be uniquely determined.

A common method for seeking $\overline{\mathcal{L}}$, referred to as the Reynolds decompositionReynolds1895 , is to apply the averaging operator directly to Equation (1) and analyze the resulting terms. For some terms in ${\mathcal{L}}$, averaging commutes with the operator in that term. For example, in Equation (2), ensemble averaging commutes with the time derivative and the diffusion operator, but not with multiplication by $u_{j}$ in the advective term. Therefore, after averaging Equation (1), and noting that $\partial u_{j}/\partial x_{j}=0$, one can write

[TABLE]

If averaging had commuted with multiplication by $u_{j}$, for example in the case of a constant $u_{j}$, then one could have written the advective term in (4) as ${\partial\left(\overline{u_{j}}\ \overline{c}\right)}/{\partial x_{j}}$, and therefore the job of finding $\overline{\mathcal{L}}$ would have been complete. However, commutation is rarely permissible when the underlying microscopic flow field is non-uniform. With this line of thinking, the problem of finding the macroscopic operator is then focused on obtaining closures to $\overline{u_{j}c}$. Consistent with this mindset, we sometimes use the notation $\overline{\mathcal{L^{\prime}}}$ to refer to the portion of $\overline{\mathcal{L}}$ that cannot be closed due to commutation error. In other words, for the case of advection-diffusion problem one can write:

[TABLE]

Use of this notation might mistakenly imply that $\overline{\mathcal{L^{\prime}}}$ is solely dependent on the microscopic flow velocity fluctuations, and independent of other processes such as molecular diffusion. We must remember that this is generally not the case; $\overline{\mathcal{L^{\prime}}}$ depends on all microscopic processes including those involving commuting operators.

When the microscopic advective process is limited to length scales much smaller than the macroscopic field, the macroscopic operator $\overline{\mathcal{L^{\prime}}}$ can be approximated by a diffusion operatorMajda1999 . However, this condition is rarely met in turbulent flows.

Before proceeding to the next section, we introduce some additional terminology. We use $\Omega$ to refer to the microscopic phase space, and $\overline{\Omega}$ to refer to the macroscopic phase space. With this definition, $c\in\Omega$ and $\overline{c}\in\overline{\Omega}\subset\Omega$. In other words, the number of the degrees of freedom in $\overline{c}$ is (much) smaller than the those in $c$ owing to elimination of degrees of freedom by the averaging operator. Averaging is a projection operator mapping points (fields) from $\Omega$ to $\overline{\Omega}$. The microscopic operator, $\mathcal{L}$, maps points within $\Omega$, and the macroscopic operator maps points within $\overline{\Omega}$. While we refer to $\overline{\mathcal{L}}$ as the macroscopic operator, we use “macroscopic closure operator” to refer to $\overline{\mathcal{L^{\prime}}}$ . In general, the macroscopic closure operator can be anisotropic, may involve non-constant coefficients, may involve spatial derivatives higher than second order, and may even be non-local.

Next, we review, as a pedagogical example, the problem of dispersion by parallel flows. In Section I.2 we present an approximate solution to this problem consistent with the methodology of G. I. TaylorTaylor1953 . Then after introducing the macroscopic forcing method (MFM) in Section II, we will revisit this problem and assess the performance of the macroscopic operators obtained using MFM and those obtained from Taylor’s solution.

I.2 Simple Example: Dispersion by a Parallel Flow

A macroscopic model for dispersion of passive scalars by parallel flows was first introduced by G. I. TaylorTaylor1953 , where he considered the evolution of a passive scalar by a laminar pipe flow. Consistent with experimental observations, his analysis predicted that the long-term evolution of the cross-sectional averaged scalar field can be predicted by a diffusion equation to the leading order with an effective diffusivity quadratically dependent on the velocity magnitude and inversely proportional to the molecular diffusivity. Since then, this problem has been extensively revisitedAris1960 ; Ajdari2006 ; Barton1983 ; Vedel2012 with some studies offering corrections to the leading-order modelGill1970 ; Barton1983 including extensions to unsteady flowsChatwin77 , and non-parallel flows, such as those in porous mediaKoch1987 ; Frankel1989 ; Mauri1991 . We will show that MFM results are consistent with these corrections, while providing a systematic and robust computational approach for the determination of the macroscopic operators in more general settings including macroscopically inhomogeneous flows.

Without loss of key physical ingredients and lesson points, we consider a simplified parallel flow in a two-dimensional domain as depicted in Figure 1.a

[TABLE]

where $U$ is the velocity amplitude and $L_{2}$ is the domain length in the $x_{2}$-direction. We consider a microscopic transport equation given by

[TABLE]

We consider the physical domain $-\infty<x_{1}<+\infty$, and $0\leq x_{2}<L_{2}$ with periodic boundary conditions in the $x_{2}$ direction. Here, we define the averaging operator

[TABLE]

Our task is to determine the macroscopic differential operator acting only in the $x_{1}$ direction that describes the evolution of $\overline{c}$ without requiring a direct solution to (7).

It is worth starting the analysis with a non-dimensionalization of coordinates. Scaling the spanwise length by $L_{2}/(2\pi)$, the streamwise length by $UL_{2}^{2}/(4\pi^{2}D_{M})$, and time by $L_{2}^{2}/(4\pi^{2}D_{M})$, results in the following dimensionless PDE:

[TABLE]

defined over a domain $-\infty<x_{1}<+\infty$, and $0\leq x_{2}<2\pi$, with periodic boundary conditions in the $x_{2}$ direction. $\epsilon=2\pi D_{M}/(L_{2}U)$ is the only dimensionless parameter of the problem, representing an inverse Peclet number, Pe. Here, we introduce our mathematical tools and analyses considering only the regime of large Pe, and for a simplified version of the microscopic equation, in which we neglect the last term on the right hand side of (9) with the intention of keeping this pedagogical toy problem brief. Later, we will present a more complex example in which finite $\epsilon$ effects are considered.

Figure 1.c shows an example solution to this problem. The macroscopic average of this solution, $\overline{c}$, is shown in Figure 1.d. The evolution equation for $\overline{c}$ can be derived by applying the averaging operator to Equation (9). In the averaging process, the first term on the right hand side of (9) vanishes due to periodic boundary conditions. Considering the aforementioned simplification, the averaged equation can be written as

[TABLE]

The first term in (10) is already closed as seen before, while the second term involves an unclosed product. Next, we derive an approximate closure to this equation following Taylor’s approach. Although Taylor did not examine zero-mean harmonic parallel flows, we still refer to the solutions presented in this section as Taylor’s solution.

In the limit that the concentration field involves features with large streamwise length $l\gg 1$, often associated with long development time, it is possible to derive an approximate solution to Equation (9) via a perturbation approach. In this case, one can conclude that the time scale of mixing in the streamwise direction, $t_{1}~{}\sim O(l)$, is much longer than the time scale of mixing in the spanwise ($x_{2}$) direction via diffusion, $t_{2}\sim O(1)$. Therefore, concentration fields can be considered almost well-mixed in the spanwise direction (or more accurately said, “rapidly developed”) while undergoing slower evolution by the streamwise flow. Introducing the decomposition $c(x_{1},x_{2})=\overline{c}(x_{1})+c^{\prime}(x_{1},x_{2})$, fast mixing in the spanwise direction implies that $c^{\prime}\ll\Delta\overline{c}$, where $\Delta\overline{c}$ represents variation of the averaged concentration in the streamwise direction. More fundamentally, rapid development in the spanwise direction, allows a quasi-steady approximation for the evolution of $c^{\prime}$. By subtracting Equation (10) from Equation (9) one can obtain an exact evolution equation for $c^{\prime}$:

[TABLE]

where the right-most prime (′) on the square bracket implies deviation from the mean of the term inside the bracket. Given $c^{\prime}\ll\Delta\overline{c}$, $\partial/\partial x_{1}\sim O(1/l)\ll 1$, and the quasi-steady assumption, the leading-order balance must be between the two terms on the left hand side. This results in a leading order solution for $c^{\prime}$ as $c^{\prime 0}=-\cos\left(x_{2}\right)\left({\partial\overline{c}}/{\partial x_{1}}\right)$. Substituting this expression into (10) results in the following macroscopic equation for the evolution of $\overline{c}$

[TABLE]

We refer to Equation (12) as Taylor’s solution.According to this approximate solution, the macroscopic operator is $\overline{\mathcal{L}}\simeq\frac{\partial}{\partial t}-\frac{1}{2}\frac{\partial^{2}}{\partial x_{1}^{2}}$, while the second term alone represents the macroscopic closure operator, $\overline{\mathcal{L}^{\prime}}=-\frac{1}{2}\frac{\partial^{2}}{\partial x_{1}^{2}}$.

Dimensionally, this equation implies that the cross-sectional averaged concentration field experiences a macroscopic diffusivity equal to $U^{2}L_{2}^{2}/(8\pi^{2}D_{M})$. The most notable outcome of this expression is the unintuitive inverse dependence on the molecular diffusivity, which is well explained in the literature and confirmed experimentally.

Before introducing our work, we briefly note that for the specific problem of parallel flows, improvements to Taylor’s solution have been investigatedGill1970 ; Barton1983 . In what is presented above, the only approximation made in the solution to (11) was ignoring the right hand side terms. One may improve this approximation by considering a series expansion for $c^{\prime}$ as $c^{\prime}=\sum_{i=0}^{\infty}c^{\prime i}(x_{1},x_{2})$ in which $c^{\prime i}\sim O(\Delta\overline{c}/l^{i+1})$, and substitution in (11). Considering $\partial/\partial x_{1}\sim O(1/l)$ and $\partial/\partial t\sim O(1/l^{2})$, as inferred from (12), and equating terms of similar order, one may obtain a recursive relation through which higher order solutions to $c^{\prime}$ can be successively obtained by substitution of lower order solutions into the right hand side of (11). Substituting the improved $c^{\prime}$ in (10), results in an improved macroscopic model. For example, the next correction results in the following macroscopic PDE:

[TABLE]

Equation (13) is a perturbative correction to Taylor’s solution, with the macroscopic closure operator represented as $\overline{\mathcal{L}^{\prime}}=-\frac{1}{2}\frac{\partial^{2}}{\partial x_{1}^{2}}+\frac{1}{2}\frac{\partial^{3}}{\partial t\partial x_{1}^{2}}-\frac{1}{32}\frac{\partial^{4}}{\partial x^{4}}$. The unfamiliar term, $\partial^{3}\overline{c}/\partial t\partial x_{1}^{2}$ can be simplified as $(1/2)\partial^{4}\overline{c}/\partial x_{1}^{4}+O(1/l)^{6}$ as suggested by (12). Substituting this result in (13) leads to a simpler macroscopic PDE, $\partial\overline{c}/\partial t=(1/2)\partial^{2}\overline{c}/\partial x_{1}^{2}-(7/32)\partial^{4}\overline{c}/\partial x_{1}^{4}$ with coefficient signs consistent with dissipation mechanism.

By following the necessary algebraic steps one may realize that even in this simple setting, the correction procedure is cumbersome and quickly becomes analytically intractable.

II The Macroscopic Forcing Method

We now resume the general problem by envisioning an arbitrary transport process described by a microscopic operator, $\mathcal{L}$. In discretized space we can write the microscopic equation as

[TABLE]

This equation is identical to Equation (1), however we use brackets to denote matrices and vectors, reminding that the problem at hand is essentially a linear algebra problem. Our goal is to determine the macroscopic operator $[\overline{\mathcal{L}}]$ that acts directly on the macroscopic field $[\overline{c}]$. The answer is straightforward once one sees the problem from a linear algebra perspective: the macroscopic solution is simply projection of the microscopic solution to the mean space, which itself can be represented as a linear operation. Proper combination of the projection operator and the microscopic operator should therefore provide the macroscopic operator. We will discuss an approach based on this mindset later in Section III.1, however, we point that this direct approach would require determination of the inverse of the microscopic operator, which is prohibitively expensive for practical problems. We therefore, present the macroscopic forcing method as an alternative method for obtaining the same $[\overline{\mathcal{L}}]$ without requiring inversion of the microscopic operator.

Given the linearity of the microscopic operator, it is straightforward to deduce that the macroscopic operator, $\overline{\mathcal{L}}$, must also be linear. We can therefore write the macroscopic equation as

[TABLE]

In order to determine $\overline{\mathcal{L}}$, we examine how the macroscopic system above responds to forcing by setting up the input-output relation

[TABLE]

where $\overline{c}\left(s\right)$ denotes the output response due to the input forcing $s$. One can reveal columns of the matrix representing $\overline{\mathcal{L}}^{-1}$ by obtaining $\overline{c}$ in response to activation of different elements in $s$. By combining these columns, it is possible to construct both $\overline{\mathcal{L}}^{-1}$ and $\overline{\mathcal{L}}$.

Note that operations above are defined in the macroscopic space and thus $s\in\overline{\Omega}$, meaning that the forcing must be macroscopic ($s=\overline{s}$). However, in the absence of $\overline{\mathcal{L}}$, one can obtain $\overline{c}\left(s\right)$ by directly solving the microscopic equation, (1), with the macroscopic forcing added to its right hand side,

[TABLE]

MFM is the procedure of determining $\overline{\mathcal{L}}$ by obtaining $\overline{c}$ in response to different macroscopic forcing scenarios. With the description given above, we obtain $\overline{\mathcal{L}}$ while bypassing the expensive step of inverting $\overline{\mathcal{L}}$, however, the expense of the procedure still remains a critical concern. Brute force MFM, as discussed above, is still expensive, since it requires many DNS calculations. We will present in sections III.3 and III.4.2 methods to substantially reduce the cost of MFM owing to certain smoothness properties of $\overline{\mathcal{L}}$.

In the description above, we skipped discussion of initial condition and boundary conditions. One may naturally incorporate these by extending the discretized system (15) to include both initial and boundary points in $[\overline{c}]$, and properly adjusting the operator $[\overline{\mathcal{L}}]$, and the right-hand-side of (15). In doing so, we implicitly assume that both initial condition and boundary conditions are macroscopic, i.e. their description in the macroscopic space is the same as that in the microscopic space. This constraint, however, does not limit the applicability of MFM to turbulence modeling since most practical problems are insensitive to initial conditions, and involve boundary conditions that are indeed macroscopic (e.g., boundary condition remains intact between Navier-Stokes and RANS descriptions). Keeping this in mind, we revisit the problem of dispersion by parallel flows by constraining the initial condition to be macroscopic, $c|_{t=0}=c\left(x_{1}\right)$, as shown in Figure 1b.

II.1 Revisiting Dispersion by a Parallel Flow

As a pedagogical example, we revisit the problem of dispersion by parallel flows. This time, we apply MFM to reveal the exact macroscopic operator. The equation to be solved is

[TABLE]

Note that since the macroscopic space does not involve $x_{2}$, $s$ is chosen to depend only on $x_{1}$ and $t$. Selection of $s$ as a delta function in space and time will result in a $\overline{c}$ field representing the Green’s function associated with the macroscopic operator. However, here, a more insightful practice would be to continue the analysis in Fourier space. As we shall see, an additional advantage of forcing in Fourier space is the smooth dependence of the macroscopic operator on the spatial wavenumber, $k$, and temporal frequency, $\omega$, allowing its accurate representation with a coarse sampling of the forcing space. Therefore, we consider

[TABLE]

After obtaining the direct solution to Equation (18) and averaging the result, we will find

[TABLE]

Given (16), $\overline{\mathcal{L}}$ then can be expressed in Fourier space as

[TABLE]

With this procedure in mind, we proceed to solve Equation (18) assuming a forcing in the form of (19). We assume a solution of the form

[TABLE]

Substituting this expression into Equation (18) results in

[TABLE]

which is the Mathieu equation. However, to better establish intuition about the results and specifically ease interpretation between Fourier and physical spaces, we seek solutions that are asymptotically expressible in terms of polynomials or simple operations acting on polynomials. Therefore, instead of expressing the solution in terms of Mathieu function we take the following semi-analytical approach. For a given $k$ and $\omega$, one may expand $\widehat{c}$ as

[TABLE]

Substituting this expansion into (22) results in a system of linear equations for the $a_{n}$s. We present a solution to this system in Appendix A. The macroscopic operator in Fourier space, $\widehat{\overline{\mathcal{L}}}$, can be determined by substituting $\overline{\widehat{c}}=a_{0}(k,\omega)$ from this solution into (21). From (10) we remember that the sole closed term in the full operator is the $\partial/\partial t$ term, and thus the macroscopic closure operator can be computed as

[TABLE]

Figures 1.e and 1.f show the solution for $\widehat{\overline{\mathcal{L}^{\prime}}}\left(\omega,k\right)$ obtained from the above procedure. We first note that in the limit of small $k$ and $\omega$ the plotted solution matches with Taylor’s solution given in (12), which can be expressed in Fourier space as $\widehat{\overline{\mathcal{L}^{\prime}}}=0.5k^{2}$ (see Figure 1.e). Further examination indicates that higher-order terms in the Taylor series of this numerical solution also match the correction to the Taylor’s solution given by (13), which is expressed in Fourier space as $\widehat{\overline{\mathcal{L}^{\prime}}}=0.5k^{2}-0.5i\omega k^{2}-(1/32)k^{4}$.

Unlike these perturbative approaches, the present MFM solution discloses the macroscopic closure operator over the entire spectrum of scales. As shown in Figures 1.e, for wavenumbers or frequencies of $~{}O(1)$ or higher, $\widehat{\overline{\mathcal{L}^{\prime}}}$ becomes significantly different from the operators from the perturbative models. In fact, even for $k$ slightly larger than one, the successive perturbative corrections lead to divergence from the true $\widehat{\overline{\mathcal{L}^{\prime}}}$ (not shown here). This is because perturbative corrections tend to approximate $\widehat{\overline{\mathcal{L}^{\prime}}}(k,\omega)$ in terms of a polynomial expansion where the order of the dominant term naturally increases as $k$ increases. In contrast, the true $\widehat{\overline{\mathcal{L}^{\prime}}}$ transitions from a second-order power-law $\widehat{\overline{\mathcal{L}^{\prime}}}\sim k^{2}$ to a lower-order power-law, $\widehat{\overline{\mathcal{L}^{\prime}}}\sim k^{1}$, as the magnitude of $k$ becomes larger as shown in the figure.

To demonstrate this point, we compare the performance of macroscopic closure models against a reference solution obtained from averaging of a DNS solution involving high wavenumbers. For the reference solution, we numerically solved the two-dimensional transport problem described by Equation (9). We considered $c(x_{1},x_{2},0)=\exp\left(-x_{1}^{2}/0.025\right)$ as an initial condition representing a thin zone of contaminant in a 2D space (see Figure 1.b). Using a second-order central difference code, and mesh resolutions $\Delta x_{1}=0.05$ and $\Delta x_{2}=0.04\pi$, the temporal evolution of the scalar field is obtained using the forth-order Runge-Kutta method with time step $\Delta t=0.002$. The resulting $\overline{c}\left(x_{1},t\right)$ from 2D DNS is shown in Figure 2.a. Figure 2.b shows the Taylor’s solution, given by Equation (12) for the same setup. It is evident that the early-time dispersion is not captured properly. More remarkably, the perturbative correction to Taylor’s solution, given by (13), predicts a worse behavior as shown in Figure 2.c. This lack of convergence has also been noted by othersFrankel1989 , but it is insightful to revisit it here by examining the operators in Fourier space.

It is needless to say that the full macroscopic closure operator, obtained from MFM, predicts a solution matching the mean of the DNS solution. While we have verified this claim numerically, this exercise is of little practical utility unless we can express the macroscopic closure operator in terms of a closed mathematical equation. In the absence of analytical solutions, we attempted to obtain analytical curve fits to the plots shown in Figure 1.e and 1.f. In doing so, we learned that capturing the true asymptotic behavior of the curves is more crucial than capturing the higher order terms in their Taylor series expansion. By examining the asymptotic behavior of $\widehat{\overline{\mathcal{L}^{\prime}}}$ in the limits of small and large $\omega$ and $k$, we obtained the following fitted expression for $\widehat{\overline{\mathcal{L}^{\prime}}}$

[TABLE]

Specifically, Equation (25) is exact in the limits of small or large $k$ or $\omega$ as shown in Figure 1.e and f. Before, developing an intuitive comprehension of this macroscopic closure model, we assess its performance against DNS as done for previous closure operators (see numerical details in Appendix B). Figure 2.d shows the space-time evolution of $\overline{c}$ obtained from a numerical solution to the macroscopic equation with the closure operator as in (25). It is evident that this MFM-inspired closure operator offers remarkable improvement compared to Taylor’s solution and its perturbative correction.

Given that the expression in (25) is not a polynomial expression, one may wonder how the macroscopic operator looks in physical space. Based on Equation (25), the macroscopic PDE in physical space can be obtained by replacing $k^{2}$ and $i\omega$ with their appropriate operator representation in physical space. These are respectively, $-\partial^{2}/\partial x_{1}^{2}$, and $\partial/\partial t$. After adding back the closed $\partial/\partial t$ term, one obtains the full macroscopic equation as:

[TABLE]

where $\mathcal{I}$ represents the identity operator. Equation (26) represents an MFM-inspired macroscopic model for dispersion of scalars by the considered parallel flow. This mathematical outcome may seem surprising at first, because the resulting operator looks very different from those describing macroscopic dispersion, such as those in Taylor’s solution and its perturbative corrections. This is one aspect of our study that we will continue to highlight repeatedly since the observed difference offers opportunity to improve closure models in use today. However, before moving forward with this mission, we build more intuition about this result by making the following observations:

1- As a sanity check, we note that the resulting model in (25) offers a spatially conservative operator. Given the form of the closure in (10), the macoscopic operator must be divergence of a flux. Here the divergence is implicitly embedded in the expression, and can be verified by checking that the closure operator vanishes in the limit of $k=0$, noting that the intended root for the expression under the square root expands as $1+i\omega+O(k^{2})$.

2- The resulting model involves a square root acting on an operator. The mathematical interpretation of this square root in Fourier space is straightforward: it implies algebraic square root acting on a number that depends on $k$ and $\omega$. In physical space, the mathematical interpretation can be eased if one thinks about how to implement this square root numerically: in the discrete space, operators are represented by matrices. The square root simply implies the square root of the matrix. We will build more physical intuition about the meaning of this square root in Section III.1. Use of fractional-order operators for closure models has been suggested beforeChen2006 ; Epps2018 . However, we here directly obtain the operator from the governing equation, as opposed to inferring it from fitting solutions to experimental data of mean fields. Additionally, the obtained operator differs from common fractional-order closure operators, since the closure term is not solely $\partial^{\alpha}\overline{c}/\partial x_{1}^{\alpha}$, with $\alpha$ being a non-integer constant. Neither do we obtain such behavior in any asymptotic limit. As we examine more cases and build more intuition, our initial excitement about fractional-order operators fades away and is replaced with a more general concept involving non-locality of operators as discussed in Section III.1.

3- In the limit of small $k$ and $\omega$, Equation (25) reduces to $\widehat{\overline{\mathcal{L}^{\prime}}}=k^{2}/2$. In other words, the macroscopic model in physical space is $\partial\overline{c}/\partial t=(1/2)\partial^{2}\overline{c}/\partial x_{1}^{2}$. This matches the Taylor’s solution, presented in Equation (12), which is appropriate for long-term dispersion of contaminants by parallel flows.

4- The key difference between the MFM-inspired closure operator and those in (12) and (13), is in their character in the high-wavenumber limit, which has close connection to their performance in predicting early-time dispersion of contaminants as shown in Figure 2. The linear propagation of the contaminant front with time, shown in Figure 2.a, is due to the fact that the underlying mechanism of transport is advective. Taylor’s solution, however, results in a diffusion equation, (12), predicting contaminant propagation as $x_{1}\sim\sqrt{t}$, which implies an unbounded dispersion speed at early times proportional to $t^{-1/2}$. In contrast, the MFM-based model, (26), honors the physically expected dispersion speeds at both early and long times. This difference in physical space, can be traced in Fourier space by noting the character of $\widehat{\overline{\mathcal{L}^{\prime}}}$ in the high wavenumber limit. Unlike Taylor’s solution, true $\widehat{\overline{\mathcal{L}^{\prime}}}$, and the proposed model in (25) exhibit a high-wavenumber scaling as $\widehat{\overline{\mathcal{L}^{\prime}}}\sim k^{1}$. While there exist other models in the literature, such as the telegraph equation Goldstein51 , that satisfy the property of finite dispersion speed, here the proposed MFM technique allows a systematic determination of the closure operator over all spatio-temporal scales. The obvious advantage is achieving quantitative accuracy, in addition to satisfying the expected qualitative trends. As a demonstration, in Appendix C we discuss a comparison of the solution obtained by the telegraph equation to our solution presented in Figure 2. Furthermore, in the steady limit ($\omega=0$), and in contrast to the telegraph equation, the true closure operator does not reduce to the standard diffusion operator, as we shall see next.

One motivation for examining the steady limit of the closure operator is that most practical turbulent flows are statistically time-stationary, and thus the RANS operator is a steady operator. Secondly, even in the context of semi-parallel laminar flows, sometimes the steady macroscopic operator is of interest. Though in such cases, instead of an initial value problem, one is interested in a boundary value problem. In some cases, Taylor’s solution is not satisfactory as one is interested in solutions beyond the small wavenumber limitYaroshchuk2011 ; Rubinstein2013 .

II.2 Steady Limit of Dispersion by a Parallel Flow

We examine the steady limit of the same problem discussed in Section II.1. Given the homogeneity of this problem in the $x_{1}$ direction, we are not yet examining a boundary value problem. We will consider statistically inhomogeneous problems with boundary conditions in Section III.

Figure 3 shows the macroscopic closure operator for the steady limit ($\omega=0$) in wavenumber space. In the small wavenumber limit, where one may assume separation of scales between mean fields and underlying flow scale, the $k^{2}$ scaling is recovered. In the large wavenumber limit, as discussed earlier, a linear $k^{1}$ scaling prevails, even in the steady limit. The transition wavenumber between these two scalings is $k=2$, which is obtained from the intersection of (dotted) asymptotic lines in Figure 3. Converting back to dimensional units, this transition wavenumber is $8\pi^{2}D_{M}/(UL_{2}^{2})$. The inverse of this scale, $UL_{2}^{2}/(8\pi^{2}D_{M})$, plays a role analogous to the “mean free path” in kinetic theory of gases: it is a length that is scaled by the product of contaminant’s streamwise characteristic velocity, $U$, and the time scale $L_{2}^{2}/(8\pi^{2}D_{M})$ for contaminants to diffusively cross a distance of order $L_{2}/2$ and switching the sign of their streamwise advective motion, hence completing a “Brownian walk”. Following PrandtlPrandtl1925 , we will refer to this length as the “mixing length” instead of a mean free path.

With this observation, the analogy with molecular to continuum macroscopic models becomes more clear. Taylor’s diffusion model, (12), suffices for macroscopic scales much larger than the mixing length, and thus analogue Knudsen number should be defined based on the mixing length. Furthermore, the macroscopic diffusion coefficient is equal to the product of the mixing length (mean free path), $UL_{2}^{2}/(8\pi^{2}D_{M})$ and the microscopic velocity $U$.

For macroscopic scales smaller than the mixing length, the boundedness of dispersion speed requires the macroscopic operator to change character to a $k^{1}$ scaling. However, unlike the advection operator, which also scales as $k^{1}$, the nature of $\widehat{\overline{\mathcal{L}^{\prime}}}$ is dissipative and not advective. In other words, the pre-factor to this proportionality is real and positive, not imaginary. Therefore, a first derivative operator $\overline{\mathcal{L}^{\prime}}\sim\partial/\partial x_{1}$ does not capture the correct trends. This is why we used a square root to capture the expected asymptotic limit as $\overline{\mathcal{L}^{\prime}}\sim\sqrt{\partial^{2}/\partial x_{1}^{2}}$. However, this should be viewed only as a quick remedy; in Section III.1 we present a generalized operator form after building more intuition about the problem.

The fitted model presented in Equation (25) reduces to $\widehat{\overline{\mathcal{L}^{\prime}}}=\sqrt{1+k^{2}}-1$ in the steady limit. In this limit, however, we were able to obtain a better fit (see Figure 3 dash-dotted),

[TABLE]

With this model the macroscopic closure operator in physical space can be written as

[TABLE]

which is a divergence of a flux as expected. The flux itself is written as a macroscopic “diffusivity” times the gradient of the mean field, $\mathcal{D}\partial/\partial x_{1}$, which is also a familiar form. We emphasize that the macroscopic diffusivity is a non-local operator,

[TABLE]

where $D=0.5$ is the diffusion coefficient from Taylor’s analysis, and $l=0.5$ is the mixing length. In the low wavenumber limit, the second derivative in the denominator is negligible, and thus the macroscopic diffusivity is a number $\mathcal{D}=D=0.5$. Conversely, in the high-wavenumber limit, macroscopic diffusivity vanishes inversely proportional to $k$. We remind that the eddy diffusivity operator introduced in  (29) is non-singular. The denominator of this operator is diagonally dominant and invertible.

Although the parallel flow problem considered here is a simplified example, the learning lessons presented have broad implications both regarding the method applied and regarding the obtained closure operator. For example, in Shirian2019 we show application of MFM to homogeneous isotropic turbulence, yet revealing a closure operator form identical to (29).

Analysis in Fourier space limits the applicability of the presented method so far to problems in which the macroscopic space is spatially homogeneous. In such cases it is appropriate to utilize harmonic functions to construct a basis representing the solution space. In the next sections we extend MFM by departing from this simplified limit by considering inhomogeneous flows and domains involving wall boundary conditions. We demonstrate that an alternative solution-space basis in terms of polynomial functions provide a cost-manageable approach to MFM analysis.

III MFM for Inhomogeneous Flows

So far we have demonstrated how MFM can be applied to flows with homogeneous directions. In these cases, given the smoothness of $\overline{\mathcal{L}^{\prime}}$ when transformed to Fourier space, we were able to save on the number of required MFM simulations. For non-homogeneous problems, however, formulation in Fourier space is inappropriate. Nevertheless, the methodology presented in Section II is extendable to non-homogeneous problems if one applies MFM in physical space. In what follows, we will first adopt this brute force approach to develop intuition about inhomogeneous macroscopic operators while ignoring the computational expense of MFM. In the following sections we will remedy the high expense of this brute force MFM, by introducing a more computationally feasible MFM, which focuses on computation of moments of eddy diffusivity kernels.

III.1 Linear Algebra Interpretation

In order to generate some mathematical insights, it is worth developing an understanding of MFM from a linear algebra point of view. The starting point is the microscopic equation,

[TABLE]

where $\left[\mathcal{L}\right]$ represents the discretized linear microscopic operator, $\left[c\right]$ is a vector representing the microscopic field, and $\left[s\right]$ is a vector representing the forcing. Square brackets indicate a matrix or vector involving a discrete set of numbers.One may represent the generalized momentum transport equation in a similar fashion by considering the field vector to represent momentum and pressure fields and the operator matrix to represent the discretized form of (49) and (50) where the boundary conditions can be embedded in $\left[\mathcal{L}\right]$.

Next, we define the averaging operator as

[TABLE]

where $\left[P\right]$ is the non-square project matrix with number of rows and columns respectively equal to the dimension of $\overline{\Omega}$ and $\Omega$. Given that the forcing is macroscopic, we have $\overline{s}=s$ in continuum space. In discrete space, however, $\left[\overline{s}\right]$ and $\left[s\right]$ have different dimensions. Therefore, before applying MFM, one needs to extend (or interpolate) $\left[\overline{s}\right]$ to the microscopic space

[TABLE]

where $\left[E\right]$ is the non-square extension matrix with number of rows and columns respectively equal to the dimension of $\Omega$ and $\overline{\Omega}$. $\left[E\right]$ and $\left[P\right]$ satisfy the relation $\left[P\right]\left[E\right]=\mathcal{I}$.

With these definitions at hand, we can obtain the macroscopic operator in terms of the defined matrices above. Combining equations (30), (31) and (32), one can write

[TABLE]

From the definition of macroscopic operator in (16), we have that $\left[\overline{c}\right]=\left[\overline{\mathcal{L}}\right]^{-1}\left[\overline{s}\right]$. Comparing this with (33) we conclude that

[TABLE]

Equation (34) is a foundational equation in the sense that it represents the macroscopic operator defined directly in terms of the microscopic operator and the definition of average which determines the projection and extension operators. Unlike MFM, this linear-algebra-based approach does not explicitly involve the force field $s$. However, computation of $\overline{\mathcal{L}}$ from this approach can be more expensive due to the cost of direct computation of $\left[\mathcal{L}\right]^{-1}$. MFM is a remedy technique that avoids costly inversion of the microscopic operator by utilizing explicit forcing in order to reveal the macroscopic operator, $\overline{\mathcal{L}}$. Nevertheless, brute force MFM is still highly expensive.

Before remedying the cost issue, we develop more insights about macroscopic operators, by examining inhomogeneous mixing of a passive scalar using the described linear-algebra method based on Equation (34). We consider a 2D domain representing a channel with the left and right walls at $x_{1}=\pm\pi$ with a Dirichlet condition $c=0$, and the top and bottom walls at $x_{2}=0,2\pi$ with a no-flux condition $\partial c/\partial x_{2}=0$. The concentration field $c\left(x,y\right)$ is governed by the steady advection diffusion equation

[TABLE]

One may interpret the unequal diffusivities in the two directions to be an outcome of directional nondimensionalization as seen in Section I.2, i.e., $\epsilon^{2}=0.05$. For this example, we consider an incompressible flow satisfying the no-penetration conditions on all walls described by (see Figure 4)

[TABLE]

Discretization of the governing PDE is performed by utilizing the second-order central difference scheme on a uniform staggered mesh using $N_{1}$ and $N_{2}$ interior grid points in the $x_{1}$ and $x_{2}$ directions, respectively. The advective fluxes on the faces are computed using a second-order interpolation and then multiplication by the divergence-free normal velocity at the face centers. Fluxes on the left and right boundaries are computed using the interior scheme after extrapolating the concentration field to ghost points half a grid size outside of the domain using a second-order extrapolation scheme that uses knowledge of the concentration field at the boundary and the two adjacent interior cells. To ensure a well-posed system, the number of degrees of freedom associated with the forcing term must be equal to those associated with the concentration field itself. Therefore, the forcing term is also defined at the cell centers corresponding to the same locations as the interior concentration fields.

For this problem we define $\overline{c}$ as the average of the scalar field in the $x_{2}$ direction. Considering $N_{1}=40$ and $N_{2}=10$, Figure 5.a shows the matrix associated with $\left[\overline{\mathcal{L}}\right]$ obtained from implementation of Equation (34) for this problem. One can see qualitative similarity between this matrix and the standard second-order diffusion operator. We next further investigate this numerical macroscopic operator to infer a continuum operator representing the eddy diffusivity as an operator, in a manner similar to Equation (29).

III.2 Eddy diffusivity as a non-local integration kernel

We confirmed computationally that $\left[\overline{\mathcal{L}}\right]$ can be written as the product of three matrices:

[TABLE]

where $D_{M}=0.05$ is the dimensionless molecular diffusivity in the $x_{1}$ direction. The right-most $\left[\partial/\partial x_{1}\right]$ is a gradient operator utilizing the aforementioned ghost point scheme, and the left-most $\left[\partial/\partial x_{1}\right]$ is a divergence operator acting on the fluxes defined on the cell faces. The middle operator is the total diffusivity, and $\left[\mathcal{D}\right]$ is the eddy-diffusivity matrix whose entries are shown in Figure 5.b. The non-locality of this operator is quantified by the off-diagonal entries in Figure 5.b. It physically implies that the macroscopic flux at a location is not just dependent on the macroscopic gradients at the same location. Instead, one needs to combine the macroscopic gradients from a neighborhood. The spatial extent of this neighborhood is scaled by the mixing length.

Translating this result back to continuum space, we conclude that eddy diffusivity is a non-local operator that outputs the (negative of the) closure flux in terms of a weighted integral acting on the mean gradient of the transported quantity. In this case $\mathcal{D}$ is the integration kernel, representing how mean gradient in one location can cause mean fluxes in another location in the domain:

[TABLE]

Extending this to multi-dimensions, eddy diffusivity is described by a tensorial kernel, $\mathcal{D}_{ji}\left({\bf x},{\bf y}\right)$, that expresses the macroscopic fluxes in terms of the macroscopic gradients as

[TABLE]

For the numerical example considered above, we show that the implied kernel from our discrete solution, $\mathcal{D}(x_{1,i},y_{1,j})=\left[\mathcal{D}\right]_{ij}/{\Delta x_{1}}$, converges by mesh refinement (here, the problem is 1D, and subscripts after the bracket imply discrete mesh points). As an example, we show in Figure 6.a the plot of $\mathcal{D}(x_{1}=0,y_{1})$ versus $y_{1}$ obtained after successive mesh refinements.

Below we report the major observations, and some speculations regarding their validity for more general settings:

1- For cases with zero normal velocity near a boundary and with boundary conditions that remain invariant under averaging, the eddy-diffusivity kernel vanishes near the boundary (see Figure 5.b). Therefore, the macroscopic operator reduces to the molecular diffusivity operator near a boundary. This implies there is no need for additional boundary conditions when scaling up to the macroscopic fields.

2- In the present example, the macroscopic space was 1D and the eddy-diffusivity involved a positive kernel. For multidimensional systems the kernel is tensorial and entries can be positive or negative. It will be interesting to examine however, whether the diagonal entries of the tensor are positive-valued, and whether the tensor is positive definite.

3- In the present example, the computed $[\mathcal{D}]$ was symmetric to the machine precision. This implies that if a mean gradient at $x_{1A}$ result in a mean flux at $x_{1B}$, the same mean gradient applied to $x_{1B}$ will result in same mean flux at $x_{1A}$. In other words, there is reciprocity in eddy diffusivity operator. It will be interesting to investigate under what conditions this reciprocity holds, and whether it is generalizable to cases with mean advection as well as cases with tensorial eddy diffusivity.

4- When the macroscopic dimensions are statistically homogeneous the eddy-diffusivity kernel can be simplified as $\mathcal{D}({\bf x},{\bf y})=\mathcal{D}({\bf x}-{\bf y})$, and thus convolution in physical space will result in multiplication in the Fourier space. The specific fitted form for the eddy-diffusivity operator identified earlier as $\mathcal{D}={D}/\sqrt{\mathcal{I}-l^{2}\partial^{2}/\partial x_{1}^{2}}$ (see Equations 29) results in a qualitatively similar kernel as that obtained for the inhomogeneous problem (see Figure 6.b). Specifically, the kernel is positive for all values of $y_{1}$ and the kernel width scales with the mixing length, $l$, as alluded to earlier.

5- In the limit that the mixing length (i.e. kernel width) is much smaller than the macroscopic length, the eddy-diffusivity kernel can be safely approximated by a Dirac delta function. Mathematically, this can be justified since the macroscopic gradient can be assumed constant over the kernel width and taken out of the integral. In other words, the integration can be approximated as multiplication by the area under the kernel. This is what Taylor’s solution, presented in Section I.2, achieves. Smallness of the mixing length is also one of the key conditions necessary for the validity of the Boussinesq approximation, but it is rarely met in realistic turbulent flows. However, even in this limit, the eddy diffusivity is not necessarily an isotropic tensor further challenging the validity of the Boussinesq approximation.

6- In earlier sections we saw universally in all examples that in the limit of large wavenumber, the macroscopic closure operator is proportional to $k^{1}$, implying that $\widehat{\mathcal{D}}\sim k^{-1}$. Inverse transforming this relation to physical space implies that the eddy-diffusivity kernel, $\mathcal{D}({\bf x},{\bf y})$, must be proportional to $\log\left(1/|{\bf x}-{\bf y}|\right)$ for small separation distances in macroscopic 1D settings as shown in Figure 6.c. Consistently $\mathcal{D}({\bf x},{\bf y})\sim 1/|{\bf x}-{\bf y}|$ in macroscopic 2D and $\mathcal{D}({\bf x},{\bf y})\sim 1/|{\bf x}-{\bf y}|^{2}$ in macroscopically 3D settings. These are all singular but integrable kernels.

7- In earlier sections we also observed that the eddy-diffusivity operator is a smooth curve in Fourier space. We here assume an infinitely differentiable curve in $k$-space. This implies that the kernel associated with the eddy diffusivity, $\mathcal{D}({\bf x},{\bf y})$, must vanish exponentially in the limit of large $|{\bf x}-{\bf y}|$ as shown in Figure 6.d.

III.3 Inverse MFM for determination of kernel moments

The insights gained from examination of the full macroscopic operator in Section III.1 allows development of a more economical method for approximate determination of the macroscopic operator. To do this, we first introduce the inverse macroscopic forcing method (IMFM). The prior macroscopic forcing method, which should be called forward MFM (FMFM), determines $\overline{c}$ in response to a specific macroscopic forcing. In contrast, IMFM obtains the macroscopic forcing required for sustaining a pre-specified $\overline{c}$. For problems in which the microscopic PDE involves a time stepping process, IMFM is straightforward. Given the PDE

[TABLE]

one simply can take the average of this equation to arrive at

[TABLE]

By expanding the time discrete form of the left hand side term, and substituting $\overline{c}\left(t+\Delta t\right)$ with the pre-specified $\overline{c}$ one obtains an explicit expression for $s=\overline{s}$ which can be substituted in the equation. For example, for the forward Euler time-stepping method $s=\left[\overline{c}\left(t+\Delta t\right)-\overline{c}\left(t\right)\right]/\Delta t-\overline{f}\left(t\right)$. Substituting this expression in the microscopic PDE constrains $\overline{c}$ to be the pre-specified value, but the non-macroscopic modes in $c$ are allowed to evolve.

For problems with a steady microscopic PDE, this treatment can be used in conjunction with artificial time stepping, where the true $s$ is evaluated in the long term limit. For turbulent flow problems that are statistically time stationary, the microscopic PDE is still time dependent. In such cases, $s$ is time-independent but is needed instantaneously. One may run multiple ensembles of the simulation and use the ensemble average at each instant. In the limit of infinite ensembles, one will obtain a time stationary $s$. However, one may run simulations using a finite number of ensembles and use the time average of the computed $s$ as an improved estimate of the true $s$. We observed this approach to be satisfactory and converge quickly with number of ensembles in our MFM analysis of turbulent channel flows (not presented here). While running simultaneous ensembles may increase computational cost per time step, it does not affect the total cost since convergence can be achieved over a proportionally lower number of time-steps.

Lastly we note that when performing IMFM, the initial condition and boundary conditions must be adjusted to be consistent with the input $\overline{c}$. For example, while the left-right boundary conditions in the example discussed in Section III.1 remain of Dirichlet type, the specified values of $c$ at the boundary must be consistent with those of $\overline{c}$ evaluated at the boundary.

The introduced IMFM methodology can be used as another brute force method to reveal columns of matrices associated with closure operators, or even in Fourier space for homogeneous problems. The immediate advantage of IMFM is that one can utilize it to directly probe a desired closure term without first encountering the full macroscopic operator. Our intention here however, is to utilize IMFM to develop a cost-effective way of determining the closure operator.

With the IMFM methodology described above one may use a small number of simulations to determine the low-order spatial moments of the eddy-diffusivity kernel, $\mathcal{D}({\bf x},{\bf y})$. For example, for the macroscopically 1D problem discussed in Section III.1, by selecting $\overline{c}=x_{1}$ and performing one IMFM DNS, and post-processing the $\overline{u_{1}^{\prime}c^{\prime}}$ data, one finds the zeroth moment of the kernel as (substitute $\partial\overline{c}/\partial x_{1}=1$ into Equation 38),

[TABLE]

In other words, with only one MFM calculation we obtain the area under the the eddy diffusivity kernel curve versus $y_{1}$ for all $x_{1}$ locations (see the plot in Figure 6.a as an example). By specifying successively higher order polynomial expressions for $\overline{c}$ versus $x_{1}$ it is possible to determine the higher spatial moments of the kernel. The first moment of the eddy-diffusivity kernel, defined below, can be computed by selecting $\overline{c}=x_{1}^{2}/2$, performing IMFM, and post processing the resulting $\overline{u_{1}^{\prime}c^{\prime}}$. Combining with the results of the zeroth moment, and analytical manipulation of Equation( 38) one can write

[TABLE]

where the last term on the right hand side can be substituted using $D^{0}$. Likewise, the second moment of $\mathcal{D}$ can be computed as

[TABLE]

One may generalize this recursive methodology to compute higher moments of $\mathcal{D}$. Determination of the leading-order moments requires few IMFM DNSs, and thus it is much less costly than the brute force MFM. Given that the kernel is smooth away from the origin ($y_{1}=x_{1}$) it is possible to accurately reconstruct it from information of the leading order moments.

Next, we discuss how the macroscopic operators can be constructed with reasonable accuracy from the knowledge of these moments.

III.4 Construction of the Macroscopic Operators from Kernel Moments

Having moments of $\mathcal{D}$ at hand, the task is to approximate the integration kernel in (39). We first consider a non-convergent, but conceptually insightful, approximation.

III.4.1 Approximation of kernel integral using Taylor series

The integral in Equation (39) can be rewritten by expanding $\partial\overline{c}/\partial x_{i}$ in terms of its spatial Taylor series around ${\bf y}={\bf x}$. Substitution and utilization of the defined moments leads to representations of $\overline{u_{j}^{\prime}c^{\prime}}$ in terms of a series involving spatial derivatives of $\overline{c}$ with the computed moments as pre-factors to each term. For example, for the test case considered in Section III.1, the unclosed flux, and thus the macroscopic closure operator, can be modeled as

[TABLE]

The leading term in this expansion, $D^{0}$, is the local approximation to the eddy-diffusivity operator; it approximates the kernel with a Dirac delta function with matching integral. Higher order terms in (45) will be present only for a non-local eddy diffusivity where the spatial moments of the kernel are finite. These terms can become significant when the macroscopic scale becomes small enough to be comparable to the kernel’s characteristic width, i.e., the mixing length.

One might expect that a truncated expansion based on a few moments beyond $D^{0}$ could result in considerable model improvement. Such expansion would still lead to a well-posed differential equation, since coefficients (to higher derivatives) vanish near the boundary as the kernel itself vanishes, and thus a model based on this expansion does not impose requirement of additional boundary conditions. However, the expansion presented in Equation (45) is not useful in practice since addition of higher moments does not lead to convergence of the solution. For example, considering the test case discussed in Section III.1, we use DNS of the following microscopic system to assess performance of macroscopic models,

[TABLE]

Figure 7.a shows the 2D DNS results and Figure 7.b shows the macroscopic fields. A macroscopic model that retains only the leading operator, $D^{0}$, performs a reasonable job in predicting the qualitative features of the mean field, while still lacking quantitatively. When expansion (45) is used, addition of the new terms did not lead to prediction improvement despite relatively smooth mean fields.

This observation is consistent with the issue seen in Section II.1 when we examined homogeneous parallel flow systems. In fact the macroscopic closure operator in the Taylor’s solution represented the zeroth moment of the macroscopic diffusivity kernel for the considered flow while its perturbative correction was indeed an analytical derivation of the higher moments in the expansion of a convolution integral (with the exception that in the analytical derivation in Section II.1 we considered the unsteady problem and thus captured both spatial and temporal moments). As previously shown in Figure 2.c, addition of terms in the moment expansion does not result in improved macroscopic prediction.

It turns out, the expansion in (45) is analogous to the Kramers-Moyal expansionVan1981 previously formulated for the statistical description of stochastic processes associated with Brownian dynamics. Lack of convergence of this expansion has been previously investigated as described by the Pawula theorem. Likewise, in these applications, while the leading moment promisingly captures a significant portion of the macroscopic behavior, addition of higher moments do not improve quantitative prediction. We are aware of one previous workMauri1991 reporting analogy to the Kramers-Moyal expansion in the continuum context in their analytical study of laminar transport of reacting scalars in porous media.

III.4.2 Approximate kernel reconstruction

Although Equation (45) is not a useful expansion in practice, the computed moments from IMFM are still valuable in the sense that one can utilize them to construct converging macroscopic operators. We next show an example of such construction. Our approach has been guided by the fact that all higher-order operators beyond $D^{0}$ are inherently non-local. We therefore construct approximate integration kernels that match the moments obtained from MFM. The details of the procedure for the example problem in Figure 7 are as follows.

We first used IMFM to compute the first three moments of $\mathcal{D}$, i.e., $D^{0},D^{1},D^{2}$. We then constructed positively valued kernels that matched the computed moments at each $x_{1}$. One approach would be to assume known kernel profiles, and tune their parameters to match the low-order moments (see Appendix D). We then used the constructed kernel to setup the macroscopic equation as,

[TABLE]

This macroscopic equation can be solved numerically after discretization. For practical problems involving 3D integrals, one may need to design iterative techniques for fast approximate calculation of the integral. Figure 7.b shows the macroscopic solution obtained through the reconstructed kernel method. Compared to prediction of the model that retains only the leading operator, $D^{0}$, the solution from the reconstructed kernel, not only shows convergence, but also significant quantitative improvement.

In the next section we discuss how the concepts developed in sections II and III can be extended to determine closure operators for momentum transport.

IV Macroscopic Forcing Method for the Navier-Stokes Equation

Extension of MFM to the Navier-Stokes equation (NS) requires discussion of two remaining issues. First is the fact that the microscopic solution is often time dependent, and turbulent. This issue can be treated by using ensemble averages. As an alternative, given that in many practical scenarios the statistical fields are time-stationary, it is sufficient to employ MFM for obtaining the steady macroscopic operator by considering time-independent forcing terms, while averaging is performed over time. The spatial inhomogeneity of such problems results in macroscopic operators that have space-dependent coefficients. Statistical anisotropy would lead to macroscopic operators of tensorial form.

The second issue, which is the more challenging one, is the nonlinearity of the Navier-Stokes equation as a microscopic model. We resolve this issue by introducing a generalized equation as follows.

Consider a velocity field, $u_{i}$, which is the solution to the incompressible Navier-Stokes equation

[TABLE]

[TABLE]

where $p$ is the pressure field normalized by fluid density, and $\nu$ is the kinematic viscosity, and $r_{i}$ is the known body force, which in most practical cases is either zero, or macroscopic, i.e., $r_{i}=\overline{r_{i}}$. For example, in turbulent channel flows $r_{i}$ represents the imposed mean pressure gradient, which is a unity vector field in the streamwise direction when reported in units of shear velocity and channel half height. $r_{i}$ may also represent the known inhomogeneous terms associated with the boundary conditions for discretized problems.

Next, we define the Generalized Momentum Transport (GMT) equation for a given $u_{i}$ as

[TABLE]

[TABLE]

where $q$ is the generalized pressure to ensure incompressibility. The boundary conditions are the same type as those for $u_{i}$, but the forcing $s_{i}$ is allowed to be different from $r_{i}$. Equation (49), subject to constraint (50), describes a linear system for $v_{i}$, and thus its macroscopic form can be obtained using MFM. In this case, ${\bf s}$ must be selected from $\overline{\Omega}$ and thus ${\bf s}=\overline{\bf s}$. Let us assume that this reduced form can be written as $\overline{\mathcal{L}}\left(\overline{\bf v}\right)=-\nabla\overline{q}+{\bf s}$ subject to the macroscopic incompressibility constraint, $\nabla.\overline{{\bf v}}=0$. The question is whether the resulting operator can produce the correct RANS solution. In other words, we should show that

[TABLE]

where ${\bf r}$ is the same body force as in (47), and subject to the same boundary conditions as those for the Navier-Stokes system. To show this condition, it is sufficient to conclude that when Equation (49) is supplemented with the same body force and boundary conditions as those in the Navier-Stokes, the resulting microscopic solution, ${\bf v}$ will have the same mean as the mean of Navier-Stokes solution, i.e., $\overline{\bf v}=\overline{\bf u}$.

It turns out, an even stronger condition holds: when the macroscopic body force of GMT is matched with that of NS, not only will the mean solution of GMT match that of NS, but the instantaneous microscopic solution of GMT also merges (exponentially) to that of NS. This condition may seem unintuitive at first, since NS often admits chaotic solutions whose long term behavior is highly sensitive to the initial condition. How could the solution to such systems merge to identical fields regardless of the initial condition? The answer lies in the linearity of GMT. A quick qualitative proof can be made by introducing the vector field ${\bf w}={\bf v}-{\bf u}$. The evolution equation for ${\bf w}$ can be obtained by subtracting NS from GMT,

[TABLE]

[TABLE]

where $\phi=q-p$ and the body force is identically zero, given ${\bf s}={\bf r}$. Additionally, since the boundary conditions of NS and GMT match, this equation is subject to homogeneous boundary conditions. Without the presence of any body force or energy injection from the boundaries, the energy in the field ${\bf w}$ shall decline over time. One may show this formally by deriving an evolution equation for the total kinetic energy of this difference field defined as $K=(1/2)\int\int\int w_{i}w_{i}dx_{1}dx_{2}dx_{3}$. Contracting (52) with ${\bf w}$ and integrating in space results in

[TABLE]

where integration by parts in conjunction with the homogeneous boundary conditions for ${\bf w}$ are used to arrive at the above equation. The only long term solution to this system is ${\bf w}=0$. Therefore, the macroscopic operator for the generalized momentum transport equation, obtained through MFM, is a correct RANS operator in the sense that it admits the true RANS solution.

The decay rate of field ${\bf w}$, indicates an important time scale that needs to be considered in MFM. When performing averaging of the microscopic solutions, one should consider time windows much larger than this decay time as a condition for convergence of statistics and independence to the initial conditions. In practice, ${\bf w}$ decays exponentially after an early transient. As an example, we performed simulations of a turbulent channel flow at $\text{Re}_{\tau}=180$, and examined the solution to the GMT, Equation (49), when it is supplemented with the same macroscopic forcing and boundary conditions as those of the NS solution but using an arbitrary initial condition. Figure 8 shows the computed norm of ${\bf w}$ versus time confirming its exponential decay. A fitted line to this plot suggests a decay time scale of about $\tau_{\text{mix}}=16.6\delta/u_{\tau}$, where $\delta$ is the channel half width, and $u_{\tau}$ is the shear scale velocity. Practically, $\tau_{\text{mix}}$ is the time scale that denotes slowest mode of mixing and can be used as a mixing index; systems with smaller $\tau_{\text{mix}}$ are better mixers. Generally, $\tau_{\text{mix}}$ for mixing of scalars can be different than those for mixing of momentum as suggested by the results shown in Shirian2019 .

One novel aspect of this work is probably recasting the Navier-Stokes equation in terms of a generalized momentum transport equation, which we called GMT (49). This generalization formalizes two interpretations of velocity as a vector quantity. The first interpretation is a kinematic interpretation; it represents the time rate of material volume crossing an interface per unit area. The second interpretation of velocity is a dynamic one; it represents momentum per unit mass which is the ability of the material to resist against forces. These two interpretations are not new. They are in fact emphasized often in graduate teaching of momentum analysis in terms of the Reynolds Transport Theorem. In RANS modeling, the goal is often prediction of forces, and thus one thinks of the dynamic interpretation of velocity. The mission of RANS models is therefore determining how the mean momentum is being transported by the underlying turbulent flow. When dealing with the nonlinear advective term, $\partial/\partial x_{j}\left(u_{j}u_{i}\right)$, the first velocity, $u_{j}$ is the kinematic velocity; it does not represent momentum, but it is a geometric “transporter” of momentum. The second velocity, $u_{i}$, is the quantity of interest; it is the “transportee,” representing momentum per unit mass. The RANS problem can be recast into the following question: how does turbulence mix momentum?

Answer to this question is dependent on the underlying turbulence. We therefore need to first specify the turbulent flow field of interest. This is done here by solving the original NS as the “donor” equation, whose solution “donates” a quantitative representation of the turbulence field of interest. The GMT equation is the “receiver” system that focuses specifically on the donated turbulence field and determines quantitatively how that field transports momentum in general. MFM analysis of GMT then provides an assessment of mean momentum transport when projected to the mean space.

In systems governed by the Navier-Stokes equation, the transporter and transportee fields are further constrained to be equal. However, this is only a special case of GMT. By obtaining an operator representing GMT, we also cover the special case of NS.

While MFM uses linear techniques, the final macroscopic RANS equation is indeed a nonlinear equation; it involves an operator acting linearly on the velocity but the operator itself is dependent on the donated velocity. The caveat is that MFM obtains the macroscopic operator in terms of the independent coordinates for a given flow. An analytical MFM solution would allow expression of the macroscopic operator directly in terms of the statistics of the underlying flow and would resolve the RANS closure problem. Numerical MFM solutions, however, fall short in accomplishing this final goal. Nevertheless, numerical MFM allows substantial advancement towards achieving this goal by determining the RANS operators in terms of independent coordinates and thus quantifying model-form errors in existing RANS models. One hope is that much of these missing pieces in current RANS models be universal and thus application of MFM on canonical problems could inform models that will be predictive of practical scenarios.

Due to the limited scope of this report, we skip providing examples of the application of MFM for analysis of turbulent flows. However, this is our main intention in developing the presented methodology.

V Discussions

We presented a systematic numerical procedure, which we call MFM, for determination of scaled up operators describing the macroscopic influence of microscopic transport processes. MFM is applicable to both laminar and turbulent flows, homogeneous and inhomogeneous flows, and can be used to examine transport of scalar fields as well as vector fields such as those describing momentum as a transported quantity.

For a given setting, MFM can obtain the exact macroscopic differential operator governing the evolution of mean fields. Additionally, MFM provides a robust framework for RANS model verification. In today’s practice, RANS model forms are postulated based on physical intuition, and their coefficients are tuned using reference data from configurations similar to those in real applications. Therefore, when a model generates accurate mean velocity profiles, it is unclear to what degree the good performance is due to cancellation of model-form error with coefficient error. Performing MFM on specific geometries, and direct comparison with RANS operators allows a more rigorous verification of RANS models.

Examination of the macroscopic operators for homogeneous flows in Fourier space revealed a smooth curve with $k^{2}$ scaling in the low-wavenumber limit and $k^{1}$ scaling in the high-wavenumber limit. The $k^{1}$ scaling was universally observed in a wide range of conditions for both laminar and turbulent problems and for both scalar and momentum transport. While the $k^{1}$ scaling is consistent with the intuitive expectation of bounded dispersion speed, the pre-factor to this scaling is real and positive, and thus cannot be captured by an advective operator. Instead, we showed that this scaling is indicative of a non-local, but dissipative mechanism of transport, as suggested by the fitted operator, $\nabla\left[D/\sqrt{1-l^{2}\nabla^{2}}\right]\nabla$. Due to presence of two asymptotic scalings in the small and large $k$ limits, the non-locality does not conform to a simple fractional-order differential operator expressed as a Laplacian to a non-integer power. Extending our observations to inhomogeneous cases, we conclude that fractional order PDEs are likely limited in scope for describing macroscopic transport. Instead one needs to consider non-locality in a broader sense, as described by the introduced kernel $\mathcal{D}\left({\bf x},{\bf y}\right)$ that expresses mean fluxes in terms of a weighted superposition of the surrounding mean field gradients.

A useful complementary study to this work would be the development of fast and low-memory computation of non-local integral operations. Ideally, to save time in both MFM and RANS calculations one should aim to utilize fast convolutions whose moments can be systematically informed by IMFM.

We acknowledge that the introduced eddy diffusivity operators in generic form as kernel integrals in Equation (38) has been recognized in previous studies in broad sense. In this context, specifically the work of HambaHamba2004 ; Hamba2005 is closest to the work presented here. However, our work involves major differences in methods of obtaining the closure operators as well as in the results developed. Specifically, the present work offers a systematic computational framework for determination of closure operators with an additional flexibility for substantial cost reduction based on method of measuring kernel moments using IMFM. Additionally, our results provided eddy diffusivity not only in terms of non-local integrals, but also we utilized the Fourier representation to arrive at expressing these operators in compact stand-alone analytical PDE form. The diffusivity operator $\mathcal{D}=D/\sqrt{1-l^{2}\nabla^{2}}$, presented for homogeneous systems, is an example of compact PDE representation of a non-local macroscopic closure which we derived and quantitatively validated. Furthermore, the present work provides a harmonious connection between closure modeling of scalar transport and that for momentum transport by extending MFM in the context of the Generalized Moment Transport equation, with provable convergence to the Navier-Stokes solution.

It might also be possible to apply MFM to the problem of LES modeling in the long term. One fundamental approach is to recast the LES modeling problem in the MFM context by modifying the definition of averages that are used here.

Lastly, MFM can be viewed as a numerical rheometer, similar to the way that molecular dynamics simulations predict continuum-level transport coefficients. Experimental rheometers measure momentum diffusivity due to the underlying chaotic Browning dynamics. When using them, the key fundamental assumption is that placing of the material under the rheometer does not affect the chaotic Browning dynamics responsible for transport of momentum. Experimental rheometry of turbulence is not possible due to violation of this condition. In this sense, a novel aspect of our work is recasting the Navier-Stokes equation in a more general form by separating the roles of velocity as a transporter and transportee. Rigorous rheometry must be non-intrusive to the transporter mechanism but it can be intrusive to the transportee field. Given emerging supercomputing power, MFM provides a fresh opportunity to study turbulence by measuring its “material properties”.

VI Acknowledgements

This work was supported by the Boeing Company under grant number SPO 134136 and by the National Science Foundation under award number 1553275. D.P. was also supported by a Kwanjeong Educational Foundation Fellowship. The authors are grateful for critical comments from Drs. Chad Winkler, Philippe Spalart, and Mori Mani through the course of the research resulting in this report. The authors are in debt to Jessie Liu and Dr. Karim Shariff for their technical discussions and elaborate comments on a draft of this report. The phrase “Generalized Momentum Transport” (GMT) as a name for Equation (49) was suggested in discussion with Dr. Karim Shariff. The terminologies of “donor” and “receiver”, in explaining MFM for momentum transport were suggested by Dr. Brandon Morgan. We are grateful to Prof. Tamer Zaki for pointing the connection between this work and the work of Dr. Fujihiro Hamba. A.M. is grateful to Dr. Parvin Mani for hosting him in San Diego during his sabbatical leave where a foundational part of this work was developed.

Appendix A Solution to Equation (22)

Substitution of Equation (23) into (22) results in a series involving $\cos\left(nx_{2}\right)$. The product $\cos\left(nx_{2}\right)\cos(x_{2})$ can be expanded as $(1/2)\cos\left[\left(n+1\right)x_{2}\right]+(1/2)\cos\left[\left(n-1\right)x_{2}\right]$. Using this expansion and balancing the left-hand side and right-hand side of the equation for each term results in

[TABLE]

where the last equation applies for $n\geq 2$. The recursive relation in (57) needs two initial conditions, i.e., $a_{1}$ and $a_{2}$. Depending on the ratio of $a_{1}$ and $a_{2}$, it admits two modes of solutions, one increasing in magnitude exponentially, and one decreasing in magnitude exponentially. To obtain a physical solution, we impose the condition that the exponentially increasing mode must be zero. This condition sets a unique value for the ratio $a_{1}/a_{2}$. We investigated this condition by rewriting (57) in terms of $r_{n}=a_{n}/a_{n-1}$ and re-arranging as

[TABLE]

The question is to determine the value of $r_{2}$ such that this recursive relation leads to $|r_{n\rightarrow\infty}|$ smaller than unity. A randomly selected $r_{2}$ leads to unbounded $|r_{n\rightarrow\infty}|$. This can be understood better when thinking in terms of Equation (57). A random combination of $a_{1}$ and $a_{2}$ will have both exponentially increasing and decreasing modes present. Even with a small pre-factor, the exponentially increasing mode eventually prevails for sufficiently large $n$. Only when the ratio of $a_{1}$ and $a_{2}$ are set to precisely remove the increasing mode, one will obtain the physically correct solution. Our numerical trick is to use the recursive relation backward in $n$. By choosing a random finite number for $r_{N}$ where $N$ is a sufficiently large number, we solve the above recursive relation backward in $n$ to find a value for $r_{2}$. Moving backward in $n$ suppresses the relative contribution of the exponentially increasing mode. Starting from sufficiently large $N$, on the order of $~{}O(100)$, we find the physical value of $r_{2}$ to be converged within many significant digits. Once the ratio $a_{2}/a_{1}$ is at hand, one can solve Equations (55) and (56) to obtain $a_{0}$.

Appendix B Numerical Solution of Equation (25)

Even though in other sections we have solved macroscopic equations directly in physical space using techniques of linear algebra, for the specific problem discussed in Section II.1, it was easier to solve the problem in the Fourier space for both spatial and temporal dimensions. To do this, we represent the initial condition as a source term that is activated at time $t=0$ via a Dirac delta function,

[TABLE]

The above expression is Fourier transformed in space and time using a physical mesh with $\Delta x_{1}=0.05$ and $\Delta t=0.025$ and a truncated domain $-5\leq x_{1}\leq 5$. $\widehat{s}\left(k,\omega\right)$ is then divided by $\widehat{\overline{\mathcal{L}}}\left(k,\omega\right)=\sqrt{\left(1+i\omega\right)^{2}+k^{2}}-1$ and then inverse transformed to obtain $\overline{c}\left(x_{1},t\right)$ in the physical space.

Next, we discuss the treatment of mode $k=\omega=0$. Given that the temporal domain is not infinite, this procedure results in a time-periodic solution corresponding to repeated spontaneous injections of scalar in the domain. To eliminate the artifacts associated with this periodicity, we applied the following remedies. First, the temporal domain is taken to be substantially large ($0\leq t\leq 40$). This time is longer than the diffusion time over the entire spatial domain, and ensures that for each period the effects of previous cycles have vanished except for their spatial mean. The second remedy is to correct the spatial mean of $\overline{c}$ to avoid a linear growth due to the periodic cycles. To this end, the mode $k=0$ is solved analytically in the physical space using the trivial differential equation $\partial/\partial t=0$.

To verify our method, we repeat this procedure using the exact $\widehat{\overline{\mathcal{L}}}$, obtained through MFM discussed in Appendix A and confirm that the obtained $\overline{c}(x_{1},t)$ matches that of the two-dimensional DNS solved in the physical space under resolution refinement.

Appendix C Comparison against solution to the telegraph equation

Following a request by anonymous referee 1, we present here the prediction of the telegraph equation for the solution to the same problem introduced in figures (1) and (2). Figure (9) shows this solution which maybe directly compared to panels (a) through (d) in Figure (2). Although this result shows finite dispersion speed at early times, indicating a performance better than the Taylor’s solution, it lacks in terms of providing key features of the solution, namely around the “edge” of the concentration profile.

The appropriate form of the telegraph equation in this case is

[TABLE]

Coincidentally, in the absence of the source term, $s$, one can show that our model presented by Equation (26) can be simplified to the above equation. However, as shown in the computational results, our model provides a significantly better solution in comparison in to the telegraph equation. The reason is that the source term cannot be ignored when manipulating models. Specifically in this case, the initial condition essentially appears as a source term with Dirac delta distribution in time and with a non-local influence on the spatio-temporal evolution of the solution.

Appendix D Kernel Reconstruction from Moments

For the problem considered in Section III.4.2, we considered an approximate kernel described by

[TABLE]

For each $x_{1}$-position, the moments of the above kernel form can be computed analytically in terms of $h$, $l_{1}$ and $l_{2}$. By matching these moments to those obtained from IMFM, we obtained profiles of $h$, $l_{1}$ and $l_{2}$ versus $x_{1}$. Aside from kernel shape approximation and discretization errors, we commit two additional errors when matching the moments. First, the analytical computation of the moments in terms of $h$, $l_{1}$ and $l_{2}$ does not consider domain truncation by walls. In other words, we assume that $l_{1}$, $l_{2}$, and $h$ become sufficiently small near the wall to allow neglecting of this error. Second, in rare points where the matched value of $l_{1}$ or $l_{2}$ becomes negative, we replace its value with zero. By utilizing such kernels we show in Figure 7 significant improvements in prediction of mean concentration field in comparison to that obtained from the local eddy diffusivity. This kernel reconstruction method also resolves the convergence issue observed when corrections to the leading operator was incorporated using expansions based on Taylor series and moments.

For the specific problem considered here, we avoided using the canonical homogeneous kernel, $\mathcal{D}=D/\sqrt{\mathcal{I}-l^{2}\partial^{2}/\partial x_{1}^{2}}$ because this kernel is symmetric in space, while the actual kernel involves spatial asymmetry near the boundaries. However, we here briefly note that in a homogeneous system the moments of this kernel are $D^{0}=D$, $D^{1}=0$, and $D^{2}=Dl^{2}/2$.