Macroscopic forcing method: a tool for turbulence modeling and analysis of closures
MMacroscopic forcing method: a tool for turbulence modeling andanalysis of closures
Ali Mani ∗ and Danah Park Department of Mechanical Engineering,Stanford University, Stanford, CA 94305, USA
Abstract
This study presents a numerical procedure, which we call the macroscopic forcing method(MFM), which reveals the differential operators acting upon the mean fields of quantities trans-ported by underlying fluctuating flows. Specifically, MFM can precisely determine the eddy diffu-sivity operator, or more broadly said, it can reveal differential operators associated with turbulenceclosure for scalar and momentum transport. We present this methodology by considering canonicalproblems with increasing complexity. As an example demonstrating the usefulness of the devel-oped methodology, we show that an eddy diffusivity operator, i.e. model form, obtained froman MFM analysis of homogeneous isotropic turbulence leads to significant improvement in RANSprediction of axisymmetric turbulent jets. We show a cost-effective generalization of MFM foranalysis of non-homogeneous and wall-bounded flows, where the eddy diffusivity is found to be aconvolution acting on the macroscopic gradient of transported quantities. We introduce MFM asan effective tool for quantitative understanding of non-Boussinesq effects and assessment of modelforms in turbulence closures, particularly, the effects associated with anisotropy and non-localityof macroscopic mixing. ∗ [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] M a y . INTRODUCTION The kinetic theory of gases has been an invaluable tool for developing continuum modelsgoverning transport of mass, momentum, and energy. Without such models, prediction offlows and heat transfer in natural and industrial systems would have required tracking ofindividual fluid molecules by solving an enormous system of coupled molecular dynamicsequations. To date, the largest of such calculations can only span domains of size of order100 nm which is by far much smaller than many engineering scales. Kinetic theory bridgedthis gap by analyzing interactions of molecules in a statistical sense. These analyses providetransport coefficients to partial differential equations (PDEs) describing transport phenom-ena in a macroscopic sense[1]. Unlike the underlying chaotic Brownian dynamics, thesescaled-up Eulerian PDEs result in smooth solutions with characteristic length scales mul-tiple orders of magnitude larger than molecular scales. These PDEs therefore offer drasticcomputational savings by lowering the required number of degrees of freedom in the systemsto 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 notcompletely 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-Stokesequation itself admits chaotic solutions spanning a wide spectrum of length scales[2], ren-dering the system-level computations prohibitively expensive. Again, as primary quantitiesof interest involve only the averaged fields, researchers have actively sought approaches forfurther scale up of the mathematical models over the past decades.Central in this context is the RANS (Reynolds-Averaged Navier-Stokes) methodology[3–6]. The goal of RANS models is to predict the average velocity field in a turbulent flowwithout the need of tracking individual chaotic eddies. In this case, the average is the meanover many ensembles of realizations of the same system under the same boundary conditions.However, as most systems involve statistically time-stationary turbulence, temporal averagesare equivalent to ensemble averages[3]. Both averages result in highly smooth mean fieldswith 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 the-ory. Most notably, the majority of the RANS models[7–11] are based on the Boussinesq2pproximation[12] suggesting closures in terms of diffusion operators with the coefficientreferred to as eddy diffusivity. By making an analogy to the kinetic theory, this approxi-mation implies that in the same way that mixing by Brownian motion leads to a scaled-upcontinuum model involving a diffusion equation with an isotropic coefficient tensor, turbu-lent mixing should be expressible in terms of diffusion operators when statistically averaged.However, unlike Brownian motion, turbulent fluctuations are often non-isotropic in practice.Aside from this issue, a major criticism of this analogy is that the diffusion approximationrequires separation of scales between the scaled-up fields and the underlying microscopicphysics. This is the case for molecular-to-continuum modeling since the continuum scale isoften orders of magnitude larger than the mean free path of molecules. However, in tur-bulent flows the underlying chaotic eddies span a wide range of scales; large eddies, whichhave 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 andscaling 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 approximation[7–11]. Among the models thatdo not use the Boussinesq approximation, some still retain diffusion-type operators in whichthe mean momentum flux is expressed as a function of local velocity gradients[13]. Even themodels that do not use any closure at the momentum transport level, such as the Reynoldsstress closure models[14], still retain similar approximations in closure of higher moments.Model-form inaccuracies in turbulence models are to some degree compensated by tuning ofthe 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 byany underlying flow field. MFM acts similar to the way that molecular dynamics simulationsreveal transport coefficients for a continuum model, with the exception that in MFM bothinput and output spaces are continua and that MFM does not make any simplifying assump-tion such as isotropy or separation of scales between turbulence and mean fields (scaled upfields). MFM precisely reveals the degree of non-locality and anisotropy of the differentialoperators governing the mean transport.In the midst of the challenge of determining accurate RANS models, decades ago, super-3omputers allowed direct numerical simulation (DNS) of canonical turbulent flows[15] andprovided unprecedented access to their detailed underlying transport processes. Despitedrastic improvements in computational power, DNS continues to remain inaccessible to in-dustrial applications, and is only viable for canonical and often academic settings. Sinceits advent, it has been expected that understanding based on DNS will inspire novel RANSmodels with superior accuracy compared to those developed prior to the DNS era. However,thus far DNS, at best, has been used to provide reference mean fields for the tuning of RANSmodel coefficients, while the model form has been set a priori. A missed opportunity hasbeen the use of DNS for revealing the model form itself. MFM closes this gap by combiningDNS with a novel statistical analysis.Lastly, the approach, which we shall introduce, unifies the problem of macroscopic model-ing 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 laminarparallel flow, first introduced by G. I. Taylor[16]. By recasting this problem in wavenum-ber space, we show that while MFM correctly reproduces Taylor’s approximation for thelow wavenumber limit, it provides a natural framework for extending his solution to allwavenumbers. We then report that the application of MFM to homogeneous turbulence,results in similar operators as those obtained from analysis of parallel flows. The paperthen continues to pave the way for the development of MFM for more complex flows, suchas those involving interactions with walls and boundary conditions, and ends by providinguseful insights about MFM from both physical and mathematical point of views. A. Problem Statement
The purpose of this section is to introduce the basic terminology and mathematicalnotation used throughout the paper. While we ultimately develop a general approach thatcan treat a wide range of transport phenomena, in this section we introduce the problem byconsidering the passive scalar advection-diffusion equation. Later, we introduce extensionsto momentum fields governed by the Navier-Stokes equation. The starting point is thecontinuum microscopic equation: L c ( x , ..., x n , t ) = 0 , (1)4here c represents the transported field, while x to x n and t are the independent variables,which correspond to the spatial and temporal coordinates, respectively. L is a differen-tial operator representing the transport physics. For example, when the advective flow isincompressible, one can write L = ∂∂t + u j ∂∂x j − ∂∂x j (cid:18) D M ∂∂x j (cid:19) , (2)where u j ( x , ..., x n , t ) represents the advecting velocity field, and D M represents moleculardiffusion. Equation (1) is accompanied by initial and boundary conditions. When (1) iswritten in discretized form, these conditions can be embedded in the operator L and theright hand side of (1).We refer to Equation (1) as the “microscopic equation” or the “microscopic transportequation,” where L is called the microscopic operator acting on the microscopic field c .Next, we define the “macroscopic” or scaled-up field, c , which represents the quantity ofengineering interest. Mathematically, c is defined as the average of c over a certain sub-set of spatiotemporal coordinates. For example, in dispersion of scalars in laminar pipeflows[16], averaging is performed in the pipe cross-section. In RANS modeling, the averag-ing 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 c , but outside ofthe scope of this report, utilizes filters or weighted averages.The problem is to determine the“macroscopic” operator, L , with an appropriate set ofboundary conditions such that the solution to equation L c = 0 (3)matches exactly the c obtained from averaging the microscopic solution. Since the macro-scopic space, has a lower dimension than the microscopic space, or better said, involves asignificantly lower number of degrees of freedom, access to the macroscopic operator allowstremendous savings in computational cost. As we shall see, L will depend on the statisticsof the microscopic advective field u j , D M , the microscopic boundary conditions, and alsothe definition of averaging that converts c to c .The immediate difficulty with the above problem statement is that there are numerousmacroscopic operators that satisfy the stated requirements. In other words, the answer to the5roblem is not unique. In what follows we will utilize additional physics-based constraintssuch that L will be uniquely determined.A common method for seeking L , referred to as the Reynolds decomposition[17], is toapply the averaging operator directly to Equation (1) and analyze the resulting terms.For some terms in L , averaging commutes with the operator in that term. For example,in Equation (2), ensemble averaging commutes with the time derivative and the diffusionoperator, but not with multiplication by u j in the advective term. Therefore, after averagingEquation (1), and noting that ∂u j /∂x j = 0, one can write ∂c∂t + ∂∂x j ( u j c ) − ∂∂x j (cid:18) D M ∂c∂x j (cid:19) = 0 . (4)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 ∂ ( u j c ) /∂x j , and therefore thejob of finding L would have been complete. However, commutation is rarely permissiblewhen the underlying microscopic flow field is non-uniform. With this line of thinking, theproblem of finding the macroscopic operator is then focused on obtaining closures to u j c .Consistent with this mindset, we sometimes use the notation L (cid:48) to refer to the portion of L that cannot be closed due to commutation error. In other words, for the case of advection-diffusion problem one can write: L = ∂∂t + u j ∂∂x j − ∂∂x j (cid:18) D M ∂∂x j (cid:19) + L (cid:48) . (5)Use of this notation might mistakenly imply that L (cid:48) is solely dependent on the microscopicflow velocity fluctuations, and independent of other processes such as molecular diffusion. Wemust remember that this is generally not the case; L (cid:48) depends on all microscopic processesincluding those involving commuting operators.When the microscopic advective process is limited to length scales much smaller than themacroscopic field, for example, a spatially periodic velocity profile with effective “mean freepath” much smaller than the domain size, the macroscopic operator L (cid:48) can be approximatedby a diffusion operator[18]. However, this condition is rarely met in turbulent flows.Before proceeding to the next section, we introduce some additional terminology. We useΩ to refer to the microscopic phase space, and Ω to refer to the macroscopic phase space.With this definition, c ∈ Ω and c ∈ Ω ⊂ Ω. The dimension of Ω, equal to the number ofdegrees of freedom of c , is much smaller than the dimension of Ω. Averaging is a projection6perator mapping points (fields) from Ω to Ω. The microscopic operator, L , maps pointswithin Ω, and the macroscopic operator maps points within Ω. While we refer to L as themacroscopic operator, we use “macroscopic closure operator” to refer to L (cid:48) . In general, themacroscopic closure operator can be anisotropic, may involve non-constant coefficients, mayinvolve 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. InSection I B we present an approximate solution to this problem consistent with the method-ology of G. I. Taylor[16]. Then after introducing the macroscopic forcing method (MFM)in Section II, we will revisit this problem and assess the performance of the macroscopicoperators obtained using MFM and those obtained from Taylor’s solution. B. Simple Example: Dispersion by a Parallel Flow
A macroscopic model for dispersion of passive scalars by parallel flows was first introducedby G. I. Taylor[16], where he considered the evolution of a passive scalar by a laminarpipe 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 diffusionequation to the leading order with an effective diffusivity quadratically dependent on thevelocity magnitude and inversely proportional to the molecular diffusivity. Since then, thisproblem has been extensively revisited[19–22] with some studies offering corrections to theleading-order model[21, 23] including extensions to non-parallel flows, such as those in porousmedia[24–26]. We will show that MFM results are consistent with these corrections, whileproviding a systematic and robust computational approach for the determination of themacroscopic operators in more general settings including macroscopically inhomogeneousflows.Without loss of key physical ingredients and lesson points, we consider a simplified parallelflow in a two-dimensional domain as depicted in Figure 1.a u = U cos (cid:18) πL x (cid:19) , u = 0 , (6)where U is the velocity amplitude and L is the domain length in the x -direction. Weconsider a microscopic transport equation given by ∂c∂t + u j ∂c∂x j = D M ∂ c∂x + D M ∂ c∂x . (7)7e consider the physical domain −∞ < x < + ∞ , and 0 ≤ x < L with periodic boundaryconditions in the x direction. Here, we define the averaging operator c ( x ) = 1 L ˆ L c ( x , x ) d x . (8)Our task is to determine the macroscopic differential operator acting only in the x directionthat describes the evolution of c without requiring a direct solution to (7).It is worth starting the analysis with a non-dimensionalization of coordinates. Scalingthe spanwise length by L / (2 π ), the streamwise length by U L / (4 π D M ), and time by L / (4 π D M ), results in the following dimensionless PDE: ∂c∂t + cos ( x ) ∂c∂x = ∂ c∂x + (cid:15) ∂ c∂x , (9)defined over a domain −∞ < x < + ∞ , and 0 ≤ x < π , with periodic boundary conditionsin the x direction. (cid:15) = 2 πD M / ( L U ) is the only dimensionless parameter of the problem.For the present, we consider the case of (cid:15) = 0, corresponding to the limit of large Pecletnumber, which is adequate for capturing the essence of the problem. Figure 1.c shows anexample solution to this problem. The macroscopic average of this solution, c , is shown inFigure 1.d. The evolution equation for c can be derived by applying the averaging operatorto Equation (9). In the averaging process, the first term on the right hand side of (9)vanishes due to periodic boundary conditions. Considering (cid:15) = 0, the averaged equation is ∂c∂t + ∂∂x (cid:16) cos( x ) c (cid:17) = 0 . (10)The first term in (10) is already closed as seen before, while the second term involves anunclosed product. Next, we derive an approximate closure to this equation following Taylor’sapproach. Although Taylor did not examine zero-mean harmonic parallel flows, we still referto the solutions presented in this section as Taylor’s solution.In the limit that the concentration field involves features with large streamwise length l (cid:29)
1, often associated with long development time, it is possible to derive an approximatesolution to Equation (9) via a perturbation approach. In this case, one can conclude thatthe time scale of mixing in the streamwise direction, t ∼ O ( l ), is much longer than the timescale of mixing in the spanwise ( x ) direction via diffusion, t ∼ O (1). Therefore, concentra-tion fields can be considered almost well-mixed in the spanwise direction (or more accuratelysaid, “rapidly developed”) while undergoing slower evolution by the streamwise flow. Intro-ducing the decomposition c ( x , x ) = c ( x ) + c (cid:48) ( x , x ), fast mixing in the spanwise direction8mplies that c (cid:48) (cid:28) ∆ c , where ∆ c represents variation of the averaged concentration in thestreamwise direction. More fundamentally, rapid development in the spanwise direction,allows a quasi-steady approximation for the evolution of c (cid:48) . By subtracting Equation (10)from Equation (9) one can obtain an exact evolution equation for c (cid:48) : ∂ c (cid:48) ∂x − cos ( x ) ∂c∂x = ∂c (cid:48) ∂t + ∂∂x [cos ( x ) c (cid:48) ] (cid:48) , (11)where the right-most prime ( (cid:48) ) on the square bracket implies deviation from the mean ofthe term inside the bracket. Given c (cid:48) (cid:28) ∆ c , ∂/∂x ∼ O (1 /l ) (cid:28)
1, and the quasi-steadyassumption, the leading-order balance must be between the two terms on the left hand side.This results in a leading order solution for c (cid:48) as c (cid:48) = − cos ( x ) ( ∂c/∂x ). Substituting thisexpression into (10) results in the following macroscopic equation for the evolution of c∂c∂t = 12 ∂ c∂x . (12)We refer to Equation (12) as Taylor’s solution (or model). Dimensionally, this equation im-plies that the cross-sectional averaged concentration field experiences a macroscopic diffusiv-ity equal to U L / (8 π D M ). The most notable outcome of this expression is the unintuitiveinverse dependence on the molecular diffusivity, which is well explained in the literature andconfirmed experimentally.Before introducing our work, we briefly note that for the specific problem of parallelflows, improvements to Taylor’s model have been investigated[21, 23]. In what is presentedabove, the only approximation made in the solution to (11) was ignoring the right handside terms. One may improve this approximation by considering a series expansion for c (cid:48) as c (cid:48) = (cid:80) ∞ i =0 c (cid:48) i ( x , x ) in which c (cid:48) i ∼ O (∆ c/l i +1 ), and substitution in (11). Considering ∂/∂x ∼ O (1 /l ) and ∂/∂t ∼ O (1 /l ), as inferred from (12), and equating terms of similarorder, one may obtain a recursive relation through which higher order solutions to c (cid:48) canbe successively obtained by substitution of lower order solutions into the right hand sideof (11). Substituting the improved c (cid:48) in (10), results in an improved macroscopic model.For example, the next correction results in the following macroscopic PDE: ∂∂t (cid:18) c + 12 ∂ c∂x (cid:19) = 12 ∂ c∂x + 132 ∂ c∂x . (13)Equation (13) is a perturbative correction to Taylor’s solution. The unfamiliar term, ∂ c/∂t∂x can be simplified as (1 / ∂ c/∂x + O (1 /l ) as suggested by (12). Substituting this9esult in (13) leads to a simpler macroscopic PDE, ∂c/∂t = (1 / ∂ c/∂x − (7 / ∂ c/∂x 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 de-scribed by a microscopic operator, L , and seek to determine the macroscopic operator, L ,such that the solution to Equation (3) results in the same c as that obtained from averagingof the microscopic solution. Given the linearity of the microscopic operator, it is straightfor-ward to deduce that the macroscopic operator, L , must also be linear. In discretized spacewe can therefore write the macroscopic equation as (cid:2) L (cid:3) [ c ] = 0 , (14)where brackets denote matrices and vectors. In order to determine L , we examine how themacroscopic system above responds to forcing by setting up the input-output relation L c ( s ) = s, (15)where c ( s ) denotes the response due to the forcing s . One can reveal columns of the matrixrepresenting L − by obtaining c in response to activation of different elements in s . Bycombining these columns, it is possible to construct both L − and L .Note that operations above are defined in the macroscopic space and thus s ∈ Ω, meaningthat the forcing must be macroscopic ( s = s ). However, in the absence of L , one can obtain c ( s ) by directly solving the microscopic equation, (1), with the macroscopic forcing addedto its right hand side, L [ c ( x , ..., x n , t )] = s. (16)MFM is the procedure of determining L by obtaining c in response to different macroscopicforcing scenarios. With the description given above, however, the immediate concern isthe expense of the procedure. Brute force MFM is indeed very expensive, since it oftenrequires many DNS calculations. We will present in sections IV B and IV C 2 methods tosubstantially reduce the cost of MFM. 10n the description above, we skipped discussion of initial condition and boundary con-ditions. One may naturally incorporate these by extending the discretized system (14) toinclude both initial and boundary points in [ c ], properly adjust the operator [ L ], and theright-hand-side of (14). In doing so, we implicitly assume that both initial condition andboundary conditions are macroscopic, i.e. their description in the macroscopic space is thesame as that in the microscopic space. This constraint, however, does not limit the ap-plicability of MFM to turbulence modeling since most practical problems are insensitiveto initial conditions, and involve boundary conditions that are indeed macroscopic (e.g.,boundary condition remains intact between Navier-Stokes and RANS descriptions). Keep-ing this in mind, we revisit the problem of dispersion by parallel flows by constraining theinitial condition to be macroscopic, c | t =0 = c ( x ), as shown in Figure 1b. A. Revisiting Dispersion by a Parallel Flow
As a pedagogical example, we revisit the problem of dispersion by parallel flows. Thistime, we apply MFM to reveal the exact macroscopic operator. The equation to be solvedis ∂c∂t + cos ( x ) ∂c∂x = ∂ c∂x + s ( x , t ) . (17)Note that since the macroscopic space does not involve x , s is chosen to depend only on x and t . Given that the problem is homogeneous in x and t , it makes sense to continue theanalysis in Fourier space. Therefore, we consider s ( x , t ) = exp ( iωt + ikx ) . (18)After averaging the direct solution to Equation (17), we will find c ( x , t ) = (cid:98) c exp ( iωt + ikx ) . (19) L then can be expressed in Fourier space as (cid:98) L = 1 / (cid:98) c ( ω, k ) . (20)With this procedure in mind, we proceed to solve Equation (17) assuming a forcing in theform of (18). We assume a solution of the form c = (cid:98) c ( ω, k, x ) exp ( iωt + ikx ) . (cid:20) iω − ∂ ∂x + ik cos ( x ) (cid:21) (cid:98) c ( ω, k ; x ) = 1 . (21)For a given k and ω , one may expand (cid:98) c as (cid:98) c ( ω, k ; x ) = ∞ (cid:88) n =0 a n ( ω, k ) cos( nx ) . (22)Substituting this expansion into (21) 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 Fourierspace, (cid:98) L , can be determined by substituting (cid:98) c = a from this solution into (20). From (10)we remember that the sole closed term in the full operator is the ∂/∂t term, and thus themacroscopic closure operator can be computed as (cid:98) L (cid:48) = 1 / (cid:98) c ( ω, k ) − iω (23)Figures 1.e and 1.f show the solution for (cid:98) L (cid:48) ( ω, k ) obtained from the above procedure.We first note that in the limit of small k and ω the plotted solution matches with Taylor’ssolution given in (12), which can be expressed in Fourier space as (cid:98) L (cid:48) = 0 . k (see Figure 1.e).Further examination indicates that higher-order terms in the Taylor series of this numericalsolution also match the correction to the Taylor’s solution given by (13), which is expressedin Fourier space as (cid:98) L (cid:48) = 0 . k − . iωk − (1 / k .Unlike these perturbative approaches, the present MFM solution discloses the macro-scopic closure operator over the entire spectrum of scales. As shown in Figures 1.e, forwavenumbers or frequencies of O (1) or higher, (cid:98) L (cid:48) becomes significantly different from theoperators from the perturbative models. In fact, even for k slightly larger than one, thesuccessive perturbative corrections lead to divergence from the true (cid:98) L (cid:48) (not shown here).This is because perturbative corrections tend to approximate (cid:98) L (cid:48) ( k, ω ) in terms of a poly-nomial expansion where the order of the dominant term naturally increases as k increases.In contrast, the true (cid:98) L (cid:48) transitions from a second-order power-law (cid:98) L (cid:48) ∼ k at small k to alower-order power-law, (cid:98) L (cid:48) ∼ k , in the large k limit as shown in the figure.To demonstrate this point, we compare the performance of macroscopic closure mod-els against a reference solution obtained from averaging of a DNS solution involving highwavenumbers. For the reference solution, we numerically solved the two-dimensional trans-port problem described by Equation (9). We considered c ( x , x ,
0) = exp ( − x / . -2 -4 -2 -2 -10 -10 -10 -10 -2 -10 -4 k k (e) (f ) L’^ ωω -1 0 1 u x (a) -2 0 2 x (b) -2 0 2 x (c) -2 0 2 x (d) FIG. 1: Schematic of a parallel flow described by u = cos( x ) (a), and an example initialcondition for a scalar field to be dispersed by this flow (b) according to Equation (9) with (cid:15) = 0. Solution field at t = 1 is shown in (c), and its projection to the macroscopic space isshown in (d). (e) and (f) respectively show the real and imaginary parts of themacroscopic closure operator, (cid:98) L (cid:48) , versus k for ω = 10 {− , − , − , , , , } . Solid lines show thethe exact values obtained through the procedure discussed in Appendix A, and the dottedsymbols show the fitted operator of Equation (24). In (e) the three first frequency casesare almost indistinguishable. The gray dashed line shows (cid:98) L (cid:48) = k / x = 0 . x = 0 . π , the temporal evolution of the scalar field is obtained using the forth-orderRunge-Kutta method with time step ∆ t = 0 . c ( x , t ) from 2D DNS isshown in Figure 2.a. Figure 2.b shows the solution obtained from Taylor’s model, (12), forthe same setup. It is evident that the early-time dispersion is not captured properly. Moreremarkably, the perturbative correction to Taylor’s model, (13), predicts a worse behavior13 t -4 -2 0 2 4 x t -4 -2 0 2 4 x (a) (b)(c) (d)DNS Leading-order TaylorHigher-order Taylor MFM-based FIG. 2: Space-time evolution of x -averaged concentration field with an initial condition c = exp ( − x / . (cid:98) L (cid:48) in the limits of small and large ω and k , we obtained14he following fitted expression for (cid:98) L (cid:48) (cid:98) L (cid:48) (cid:39) (cid:113) (1 + iω ) + k − (1 + iω ) . (24)Specifically, Equation (24) is exact in the limits of small and large k and ω as shown inFigure 1.e and f. Before, developing an intuitive comprehension of this macroscopic closuremodel, we assess its performance against DNS as done for previous closure operators (seenumerical details in Appendix B). Figure 2.d shows the space-time evolution of c obtainedfrom a numerical solution to the macroscopic equation with the closure operator as in (24). Itis evident that this MFM-inspired closure operator offers remarkable improvement comparedto Taylor’s model and its perturbative correction.Given that the expression in (24) is not a polynomial expression, one may wonder howthe macroscopic operator looks in physical space. Based on Equation (24), the macroscopicPDE in physical space can be obtained by taking the inverse Fourier transform of (24), andadding back the closed ∂/∂t term: (cid:115)(cid:18) I + ∂∂t (cid:19) − ∂ ∂x − I c ( x , t ) = s ( x , t ) , (25)where I represents the identity operator. Equation (25) represents an MFM-inspired macro-scopic model for dispersion of scalars by the considered parallel flow. This mathematicaloutcome may seem surprising at first, because the resulting operator looks very differentfrom those describing macroscopic dispersion, such as those in Taylor’s solution and its per-turbative corrections. This is one aspect of our study that we will continue to highlightrepeatedly since the observed difference offers opportunity to improve closure models in usetoday. However, before moving forward with this mission, we build more intuition aboutthis result by making the following observations:1- As a sanity check, we note that the resulting model in (24) offers a spatially conservativeoperator. Given the form of the closure in (10), the macoscopic operator must be divergenceof a flux. Here the divergence is implicitly embedded in the expression, and can be verifiedby checking that the closure operator vanishes in the limit of k = 0.2- The resulting model involves a square root acting on an operator. The most straight-forward interpretation of the square root in physical space is to consider a discretized systemfor which operators are represented by matrices. The square root simply implies the square15oot of the matrix. We will build more physical intuition about the meaning of this squareroot in Section IV A. Use of fractional-order operators for closure models has been suggestedbefore[27, 28]. 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. Addi-tionally, the obtained operator differs from common fractional-order closure operators, sincethe closure term is not solely ∂ α c/∂x α , with α being a non-integer constant. Neither do weobtain such behavior in any asymptotic limit. As we examine more cases and build moreintuition, our initial excitement about fractional-order operators fades away and is replacedwith a more general concept involving non-locality of operators as discussed in Section IV A.3- In the limit of small k and ω , Equation (24) reduces to (cid:98) L (cid:48) = k /
2. In other words,the macroscopic model in physical space is ∂c/∂t = (1 / ∂ c/∂x . This matches the Tay-lor’s model, presented in Equation (12), which is appropriate for long-term dispersion ofcontaminants 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 theirperformance in predicting early-time dispersion of contaminants as shown in Figure 2. Thelinear propagation of the contaminant front with time, shown in Figure 2.a, is due to thefact that the underlying mechanism of transport is advective. Taylor’s solution, however,results in a diffusion equation, (12), predicting contaminant propagation as x ∼ √ t , whichimplies an unbounded dispersion speed at early times proportional to t − / . In contrast,the MFM-based model, (25), honors the physically expected dispersion speeds at both earlyand long times. This difference in physical space, can be traced in Fourier space by notingthe character of (cid:98) L (cid:48) in the high wavenumber limit. Unlike Taylor’s model, true (cid:98) L (cid:48) , andthe proposed model in 24 exhibit a high-wavenumber scaling as (cid:98) L (cid:48) ∼ k . Although theconsidered example is highly simplified, this takeaway message is directly generalizable tomore complex flows.In the next section we build more intuition about these results by focusing on the steady-state limit of this problem. One motivation is that most practical turbulent flows are sta-tistically time-stationary, and thus the RANS operator is a steady operator. Secondly, evenin the context of semi-parallel laminar flows, sometimes the steady macroscopic operatoris of interest. Though in such cases, instead of an initial value problem, one is interestedin a boundary value problem. In some cases, Taylor’s solution is not satisfactory as one is16 -1 L’^ k FIG. 3: Macroscopic closure operator, (cid:98) L (cid:48) , versus k for the steady limit, ω = 0, of themicroscopic equation (9). The solid line shows the exact operator obtained from the MFMprocedure discussed in Appendix A; the dashed line is (cid:98) L (cid:48) = √ k −
1; thedashed-dotted line shows (cid:98) L (cid:48) = k / √ k . The thin dotted lines show the asymptoticlimits of 0 . k , and k , respectively in the small and large wavenumbers. The imaginarypart of all solutions is zero.interested in solutions beyond the small wavenumber limit[29, 30]. B. Steady Limit of Dispersion by a Parallel Flow
We examine the steady limit of the same problem discussed in Section II A. Given thehomogeneity of this problem in the x direction, we are not yet examining a boundary valueproblem. We will consider statistically inhomogeneous problems with boundary conditionsin Section IV.Figure 3 shows the macroscopic closure operator for the steady limit ( ω = 0) in wavenum-ber space. In the small wavenumber limit, where one may assume separation of scalesbetween mean fields and underlying flow scale, the k scaling is recovered. In the largewavenumber limit, as discussed earlier, a linear k scaling prevails, even in the steady limit.Converting the problem back to dimensional units, the transition wavenumber scales as4 π D M / ( U L ). The inverse of this scale, U L / (4 π D M ), can be interpreted as a “meanfree path”: it is the product of contaminant’s streamwise characteristic velocity, U , and the17ime scale L / (4 π D M ) for contaminants to diffusively cross a distance of order L henceswitching the sign of their streamwise advective motion.With this interpretation, the analogy with molecular to continuum macroscopic modelsbecomes more clear. Taylor’s diffusion model, (12), suffices for macroscopic scales muchlarger than the “mean free path.” In this case the macroscopic diffusion coefficient is pro-portional to the product of “mean free path,” U L / (4 π D M ) and the microscopic velocity U . Conversely, for macroscopic scales smaller than the “mean free path”, the boundednessof dispersion speed requires the macroscopic operator to change character to a k scaling.However, unlike the advection operator, which also scales as k , the nature of (cid:98) L (cid:48) is dissi-pative and not advective. In other words, the pre-factor to this proportionality is real andpositive, not imaginary. Therefore, a first derivative operator L (cid:48) ∼ ∂/∂x does not capturethe correct trends. This is why we used a square root to capture the expected asymptoticlimit as L (cid:48) ∼ (cid:112) ∂ /∂x . However, this should be viewed only as a quick remedy; in Sec-tion IV A we present a generalized operator form after building more intuition about theproblem.The fitted model presented in Equation (24) reduces to (cid:98) L (cid:48) = √ k − (cid:98) L (cid:48) = 0 . k √ . k . (26)With this model the macroscopic closure operator in physical space can be written as L (cid:48) = − ∂∂x . (cid:113) I − ∂ ∂x ∂∂x , (27)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, D ∂/∂x , which is also a familiar form.We emphasize that the macroscopic diffusivity is a non-local operator, D = 0 . (cid:113) I − ∂ ∂x . (28)In the low wavenumber limit, the second derivative in the denominator is negligible, andthus the macroscopic diffusivity is a number D = 0 .
5. Conversely, in the high-wavenumberlimit, macroscopic diffusivity vanishes inversely proportional to k . We remind that the eddydiffusivity operator introduced in (28) is non-singular. The denominator of this operator isdiagonally dominant and invertible. 18 . Dispersion of a Solenoidal Passive Vector Field The MFM methodology described in Section II can be applied to equations for transportof vector fields. For example, the transport equation for a solenoidal passive vector field, v i ,can be forced by the macroscopic vector field, s i , in order to reveal the averaged response ∂v i ∂t + ∂u j v i ∂x j = − ∂q∂x i + ν ∂ v i ∂x j ∂x j + s i , (29) ∂v j ∂x j = 0 . (30)Here, u j is a known transporter velocity field that transports v i , ν is the molecular (micro-scopic) diffusivity, and q is the scalar potential to satisfy the divergence-free condition for v i .Passive vector fields have been studied before in the context of the vorticity equation and thescalar gradient transport equation, and also to generate insights about turbulence[31, 32].Here we examine a different evolution equation (29, 30) for passive vector fields that is moreappropriate as a stepping stone for extending MFM to the Navier-Stokes equation, whichwe will discuss in Section IIIIn this section, we apply MFM to parallel flows similar to what we considered earlier forscalar transport. In the simplest form, we consider u = cos( x ), u = 0, molecular diffusivityonly in the x direction, and ω = 0 to focus on the steady response. Given that averagesare with respect to the x direction, we seek a macroscopic operator that describes thebehavior of v i ( x ) in response to a macroscopic forcing s i . Due to the continuity constraint, ∂v /∂x = 0, the solution for v reduces to a trivial form. Thus it is sufficient to study theresponse of v ( x ) to a macroscopic forcing s ( x ). The details of the solution methodologyare discussed in Appendix C.Figure 4.d shows the resulting macroscopic operator in wavenumber space, and its asymp-totic behavior in the limit of small and large k . Interestingly, here the leading-order macro-scopic closure operator (limit of k (cid:39)
0) is fourth-order, indicating a macroscopic equationof the form 0 . ∂ v /∂x = s ( x ). This is different from the second-order diffusion model ofTaylor obtained for dispersion of scalars by the same flow. However, as we shall see, this isan exceptional case due to the specific transporter field, u j . In other cases, the leading-orderoperator is often second-order.Another important observation is the behavior of the macroscopic operator in the largewavenumber limit. Similar to the scalar problem we observe that (cid:98) L (cid:48) ∼ k in this limit. This is19 -1 -4 -2 L’^ k x x x -1-0.500.51 x x x (a) (b)(c) (d) FIG. 4: An example macroscopic force field s = cos(2 x ), s = 0, corresponding to k = 2plotted in physical space (a), and response of v to this force (b) as a steady solution toEquation (29). The solid lines in (b) show the streamlines of the v field, while thetransporter velocity u is a parallel flow given by u = cos ( x ), u = 0. Projection of the v response to the macroscopic space is obtained after averaging in the x direction and isshown in (c). Panel (d) shows the macroscopic closure operator, (cid:98) L (cid:48) = s /v , versus k forthe steady limit, ω = 0. The solid line shows the exact operator obtained from the MFMprocedure discussed in Appendix C; the dashed line is (cid:98) L (cid:48) = 0 . k ; the dashed-dotted lineshows (cid:98) L (cid:48) = k .again a macroscopically non-local behavior. For the scalar problem we justified this behaviorintuitively based on the boundedness of the microscopic dispersion speed. For transportedsolenoidal vector fields, however, the pressure potential allows instant propagation of infor-mation. Nevertheless, the k scaling of the macroscopic operator persists.In the next sections we extend MFM by departing from the simplified limit of dispersionby parallel flows which we used here for pedagogy. While we consider more complex flowsinvolving 3D and unsteady velocity fields, we arrive at more general but qualitatively similar20onclusions. In this context we find that the large-eddy size in turbulence plays a roleanalogous to the mean free path in parallel flows. For macroscopic scales on the order ofor smaller than the large eddy, macroscopic operators depict considerable non-local effects,which are missing in the mainstream eddy diffusivity closure models. An additional issueis anisotropy of the macroscopic operators in multi-dimensions which persists even whenthe macroscopic scale is large. We show how MFM can measure both non-locality andanisotropy of macroscopic operators in a more general setting. Before these extensions, wefirst formally introduce MFM for the Navier-Stokes equation. III. MACROSCOPIC FORCING METHOD FOR THE NAVIER-STOKES EQUA-TION
Extension of MFM to the Navier-Stokes equation (NS) requires discussion of two remain-ing issues. First is the fact that the microscopic solution is often time dependent, inho-mogeneous, and turbulent. This issue can be treated by using ensemble averages. Giventhat in many practical scenarios, statistical fields are time-stationary, it is sufficient to em-ploy MFM for obtaining the steady macroscopic operator by considering time-independentforcing terms, while averaging is performed over time. The spatial inhomogeneity of suchproblems results in macroscopic operators that have space-dependent coefficients. Statisticalanisotropy 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 generalizedequation as follows.Consider a velocity field, u i , which is the solution to the incompressible Navier-Stokesequation ∂u i ∂t + ∂u j u i ∂x j = − ∂p∂x i + ν ∂ u i ∂x j ∂x j + r i , (31) ∂u j ∂x j = 0 , (32)where p is the pressure field normalized by fluid density, and ν 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 = r i . For example, in turbulent channel flows r i represents the imposed mean pressuregradient, which is a unity vector field in the streamwise direction when reported in units21f shear velocity and channel half height. r i may also represent the known inhomogeneousterms associated with the boundary conditions for discretized problems.Next, following the previous section, we define the Generalized Momentum Transport(GMT) equation for a given u i as ∂v i ∂t + ∂u j v i ∂x j = − ∂q∂x i + ν ∂ v i ∂x j ∂x j + s i , (33) ∂v j ∂x j = 0 , (34)where q is the generalized pressure to ensure incompressibility. The boundary conditionsare the same type as those for u i , but the forcing s i is allowed to be different from r i .Equation (33), subject to constraint (34), describes a linear system for v i , and thus itsmacroscopic form can be obtained using MFM. In this case, s must be selected from Ω andthus s = s . Let us assume that this reduced form can be written as L ( v ) = −∇ q + s subject to the macroscopic incompressibility constraint, ∇ . v = 0. The question is whetherthe resulting operator can produce the correct RANS solution. In other words, we shouldshow that L ( u ) = −∇ p + r , (35)where r is the same body force as in (31), and subject to the same boundary conditions asthose for the Navier-Stokes system. To show this condition, it is sufficient to conclude thatwhen Equation (33) is supplemented with the same body force and boundary conditions asthose in the Navier-Stokes, the resulting microscopic solution, v will have the same mean asthe mean of Navier-Stokes solution, i.e., v = u .It turns out, an even stronger condition holds: when the macroscopic body force of GMTis matched with that of NS, not only will the mean solution of GMT match that of NS, butthe 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 whoselong term behavior is highly sensitive to the initial condition. How could the solution tosuch systems merge to identical fields regardless of the initial condition? The answer lies inthe linearity of GMT. A quick qualitative proof can be made by introducing the vector field w = v − u . The evolution equation for w can be obtained by subtracting NS from GMT, ∂w i ∂t + ∂u j w i ∂x j = − ∂φ∂x i + ν ∂ w i ∂x j ∂x j (36)22
50 100 150 200 t -4 -2 | w | m a x (a) (b) (c) FIG. 5: Instantaneous iso-contour of the Q-criterion at t = 0 (a) and t = 5 δ/u τ (b) for w = v − u representing the difference between the NS solution to a turbulent channel flowat Re τ = 180 and a GMT solution with matching boundary condition and forcing. Theinitial condition for GMT is v = . Maximum magnitude of the difference field, | w | max versus time is shown in (c). ∂w j ∂x j = 0 , (37)where the body force is identically zero, given s = r . Additionally, since the boundaryconditions of NS and GMT match, this equation is subject to homogeneous boundary con-ditions. Without the presence of any body force or energy injection from the boundaries,the energy in the field w shall decline over time. One may show this formally by de-riving an evolution equation for the total kinetic energy of this difference field defined as K = (1 / ´ ´ ´ w i w i dx dx dx . Contracting (36) with w and integrating in space resultsin ∂K∂t = − ν ˆ ˆ ˆ (cid:20) ∂w i ∂x j ∂w i ∂x j (cid:21) dx dx dx , (38)where integration by parts in conjunction with the homogeneous boundary conditions for w are used to arrive at the above equation. The only long term solution to this systemis w = 0. Therefore, the macroscopic operator for the generalized momentum transportequation, obtained through MFM, is a correct RANS operator in the sense that it admitsthe true RANS solution.The decay rate of field w , indicates an important time scale that needs to be consideredin MFM. When performing averaging of the microscopic solutions, one should considertime windows much larger than this decay time as a condition for convergence of statisticsand independence to the initial conditions. In practice, w decays exponentially after anearly transient. As an example, we performed simulations of a turbulent channel flow at23e τ = 180, and examined the solution to the GMT, Equation (33), when it is supplementedwith the same macroscopic forcing and boundary conditions as those of the NS solution butusing an arbitrary initial condition. Figure 5 shows the computed norm of w versus timeconfirming its exponential decay. A fitted line to this plot suggests a decay time scale ofabout τ mix = 16 . δ/u τ , where δ is the channel half width, and u τ is the shear scale velocity.Practically, τ mix is the time scale that denotes slowest mode of mixing and can be used asa mixing index; systems with smaller τ mix are better mixers. Generally, τ mix for mixing ofscalars can be different than those for mixing of momentum as suggested by the resultsshown in [33]. A. A Philosophical Discussion on GMT
The most novel aspect of this work is probably recasting the Navier-Stokes equation interms of a generalized momentum transport equation, which we called GMT (33). Thisgeneralization formalizes two interpretations of velocity as a vector quantity. The firstinterpretation is a kinematic interpretation; it represents the time rate of material volumecrossing an interface per unit area. The second interpretation of velocity is a dynamicone; it represents momentum per unit mass which is the ability of the material to resistagainst forces. These two interpretations are not new. They are in fact emphasized often ingraduate teaching of momentum analysis in terms of the Reynolds Transport Theorem. InRANS modeling, the goal is often prediction of forces, and thus one thinks of the dynamicinterpretation of velocity. The mission of RANS models is therefore determining how themean momentum is being transported by the underlying turbulent flow. When dealing withthe nonlinear advective term, ∂/∂x j ( u j u i ), the first velocity, u j is the kinematic velocity; itdoes not represent momentum, but it is a geometric “transporter” of momentum. The secondvelocity, u i , is the quantity of interest; it is the “transportee,” representing momentumper unit mass. The RANS problem can be recast into the following question: how doesturbulence mix momentum? MFM plus GMT answers this question precisely by interpretingthe turbulence as a transporter of momentum.In systems governed by the Navier-Stokes equation, the transporter and transportee fieldsare further constrained to be equal. However, this is only a special case of GMT. By obtainingan operator representing GMT, we also cover the special case of NS.24hile MFM uses linear techniques, the final macroscopic RANS equation is indeed anonlinear equation; it involves an operator acting linearly on the velocity but the operatoritself is dependent on the velocity. The caveat is that MFM obtains the macroscopic operatorin terms of the independent coordinates for a given flow. An analytical MFM solutionwould allow expression of the macroscopic operator directly in terms of the statistics of theunderlying flow and would resolve the RANS closure problem. Numerical MFM solutions,however, fall short in accomplishing this final goal. Nevertheless, numerical MFM allowssubstantial advancement towards achieving this goal by determining the RANS operatorsin terms of independent coordinates and thus quantifying model-form errors in existingRANS models. In other words, numerical MFM acts as an advanced rheometer determininghow a given flow mixes momentum. Whether this mixing is describable by local meangradients or non-local operators, and whether the tensorial coefficients are isotropic or thereare substantial non-Boussinesq effects, MFM reveals these details quantitatively. One hope isthat much of these missing pieces in current RANS models be universal and thus applicationof MFM on canonical problems could inform models that will be predictive of practicalscenarios.We show a preliminary example of this practice next. We demonstrate that the appli-cation of MFM to homogeneous and isotropic turbulence guides a RANS model correction,leading to significant improvement in RANS prediction of turbulent round jets without theneed of applying MFM to the turbulent jet flow itself. B. MFM Applied to Homogeneous Isotropic Turbulence
Here we briefly include results of a companion study by Shirian and Mani[33] in which theMFM methodology is applied to determine the macroscopic closure operators for transportof passive scalars and momentum in stationary homogeneous isotropic turbulence (HIT).Remarkably, this study also reveals a macroscopic closure operator with similar asymptoticbehavior to what was observed for the case of scalar dispersion by a parallel flow. In thecase of turbulence, the transition between (cid:98) L (cid:48) ∼ k and (cid:98) L (cid:48) ∼ k occurs when k is on theorder of inverse large eddy size, which plays a role analogous to the role of mean free path inprevious cases. Using their MFM data, Shirian and Mani [33] obtained fitted macroscopic25losure operators in physical space of the form L (cid:48) = −∇ . (cid:34) D ( I − l ∇ ) / ∇ (cid:35) , (39)where for the highest Reynolds number considered, they reported D = 0 . u /ε , l =1 . u /ε for the case of scalar transport and D = 0 . u /ε , l = 0 . u /ε for the caseof momentum transport. Here, u rms and ε are respectively the root mean square of a singlecomponent velocity fluctuations and the turbulence dissipation rate. Equation (39) impliesthat the macroscopic diffusivity, which from now on we refer to as “eddy diffusivity”, is nota constant number. Instead, for HIT it is a non-local but isotropic operator, here expressedas D = D ( I − l ∇ ) / . (40) C. Impact on Prediction of Practical Flows
We next show the predictive impact of the obtained RANS operator by applying it toa practical flow and comparing its prediction against available experimental data. For thistest case, we consider the self-similar turbulent round jet at a high Reynolds number. Wethen select an available eddy-diffusivity-based turbulence model, and replace its closureoperator with Equation (39). For this purpose, we adopt the Prandtl Mixing Length Model(PMLM)[34] since it provides D and l as explicit functions of spatial coordinates. We obtainsolutions to both the original and modified RANS models. Figure 6.b shows a remarkableimprovement in RANS predictions, almost coinciding with the experimental data[35]. Moredetails on this analysis are given in Appendix D. IV. MFM FOR INHOMOGENEOUS FLOWS
So far we have demonstrated how MFM can be applied to flows with homogeneous di-rections. In these cases, given the smoothness of L (cid:48) when transformed to Fourier space, wewere able to save on the number of required MFM simulations. For non-homogeneous prob-lems, however, formulation in Fourier space is inappropriate. Nevertheless, the methodologypresented in sections II and III is extendable to non-homogeneous problems if one appliesMFM in physical space. In what follows, we will first develop intuition about inhomogeneous26 r/(x-x ) -1 u / u c __ r x u _ (b)(a) FIG. 6: Schematic of the self similar turbulent round jet (a), and its mean velocity profile(b). Symbols show the LDA measured data of Hussein et al.[35], the dashed line showsRANS prediction using Prandtl mixing length model (PMLM), and the solid line showsRANS prediction using a closure operator of the form shown in Equation (39) withcoefficients scaled consistently with the PMLM prescription.macroscopic operators while ignoring the computational expense of MFM. In the followingsections we will remedy the high expense of this brute force MFM, by introducing a morecomputationally feasible MFM, which focuses on computation of moments of eddy diffusivitykernel.
A. Linear Algebra Interpretation
In order to generate some mathematical insights, it is worth developing an understandingof MFM from a linear algebra point of view. The starting point is the microscopic equation,[ L ] [ c ] = [ s ] , (41)where [ L ] represents the discretized linear microscopic operator, [ c ] is a vector representingthe microscopic field, and [ s ] is a vector representing the forcing. Square brackets indicate amatrix or vector involving a discrete set of numbers. We refer to matrix and vector elementsby subscripts after the brackets. One may represent the generalized momentum transportequation in a similar fashion by considering the field vector to represent momentum andpressure fields and the operator matrix to represent the discretized form of (33) and (34)where the boundary conditions can be embedded in [ L ].27ext, we define the averaging operator as[ c ] = [ P ] [ c ] , (42)where [ P ] is the non-square project matrix with number of rows and columns respectivelyequal to the dimension of Ω and Ω. Given that the forcing is macroscopic, we have s = s incontinuum space. In discrete space, however, [ s ] and [ s ] have different dimensions. Therefore,before applying MFM, one needs to extend (or interpolate) [ s ] to the microscopic space[ s ] = [ E ] [ s ] , (43)where [ E ] is the non-square extension matrix with number of rows and columns respectivelyequal to the dimension of Ω and Ω. [ E ] and [ P ] satisfy the relation [ P ] [ E ] = I .With these definitions at hand, we can obtain the macroscopic operator in terms of thedefined matrices above. Combining equations (41), (42) and (43), one can write[ c ] = [ P ] [ L ] − [ E ] [ s ] . (44)From the definition of macroscopic operator in (15), we have that [ c ] = (cid:2) L (cid:3) − [ s ]. Com-paring this with (44) we conclude that (cid:2) L (cid:3) = (cid:8) [ P ] [ L ] − [ E ] (cid:9) − . (45)Equation (45) is a foundational equation in the sense that it represents the macroscopicoperator defined directly in terms of the microscopic operator and the definition of averagewhich determines the projection and extension operators. Unlike the brute force approach,introduced in Section II, this linear-algebra-based approach does not explicitly involve theforce field s . However, computation of L from this approach can be more expensive due to thecost of direct computation of [ L ] − . The computational approach introduced in Section II,utilizes DNS of the microscopic equations in response to explicit forcing in order to avoiddirect computation of [ L ] − . Nevertheless, both approaches remain 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-based MFM. We consider a 2D domain representing a channel with the left and right wallsat x = ± π with a Dirichlet condition c = 0, and the top and bottom walls at x = 0 , π x x FIG. 7: Flow streamlines in a 2D channel. Velocity field is given by Equation (47).with a no-flux condition ∂c/∂x = 0. The concentration field c ( x, y ) is governed by thesteady advection diffusion equation u ∂c∂x + u ∂c∂x = 0 . ∂ c∂x + ∂ c∂x (46)One may interpret the unequal diffusivities in the two directions to be an outcome of di-rectional nondimensionalization as seen in Section I B, i.e., (cid:15) = 0 .
05. For this example,we consider an incompressible flow satisfying the no-penetration conditions on all wallsdescribed by (see Figure 7) u = [1 + cos ( x )] cos ( x ) , u = sin ( x ) sin ( x ) . (47)Discretization of the governing PDE is performed by utilizing the second-order centraldifference scheme on a uniform staggered mesh using N and N interior grid points in the x and x directions, respectively. The advective fluxes on the faces are computed using asecond-order interpolation and then multiplication by the analytical normal velocity at theface centers. Fluxes on the left and right boundaries are computed using the interior schemeafter extrapolating the concentration field to ghost points half a grid size outside of thedomain using a second-order extrapolation scheme that uses knowledge of the concentrationfield 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 thoseassociated with the concentration field itself. Therefore, the forcing term is also defined atthe cell centers corresponding to the same locations as the interior concentration fields.For this problem we define c as the average of the scalar field in the x direction. Con-sidering N = 40 and N = 10, Figure 8.a shows the matrix associated with (cid:2) L (cid:3) obtained29 -2-1.5-1-0.500.511.52
10 20 30 4010203040 (a) (b)
FIG. 8: (a) (cid:2) L (cid:3) for the example problem described in Section IV A using N = 40 and N = 10. The color bar is intentionally truncated at − D ], for the same problem.from the linear-algebra-based MFM described by Equation (45). One can see qualitativesimilarity between this matrix and the standard second-order diffusion operator. Quantita-tively however, the obtained macroscopic operator involves a non-local eddy diffusivity. Toshow this, we confirmed computationally that (cid:2) L (cid:3) can be written as the product of threematrices: (cid:2) L (cid:3) = − [ ∂/∂x ] [ D + D M I ] [ ∂/∂x ] , (48)where D M = 0 .
05 is the molecular diffusivity in the x direction. The right-most [ ∂/∂x ]is a gradient operator utilizing the aforementioned ghost point scheme, and the left-most[ ∂/∂x ] is a divergence operator acting on the fluxes defined on the cell faces. The middleoperator is the total diffusivity, and [ D ] is the eddy-diffusivity matrix whose entries areshown in Figure 8.b. The non-locality of this operator is quantified by the off-diagonalentries in Figure 8.b. It physically implies that the macroscopic flux at a location is notjust dependent on the macroscopic gradients at the same location. Instead, one needs tocombine the macroscopic gradients from a neighborhood.Translating this result back to continuum space, we conclude that eddy diffusivity isdescribed by a convolution kernel acting on the mean gradients, − u (cid:48) c (cid:48) ( x ) = ˆ y D ( x , y ) ∂c∂x | y dy . (49)Extending this to multi-dimensions, eddy diffusivity is described by a tensorial kernel,30 y D ( x = , y ) -3 -2 -1 0 1 2 3 y -x D ( x , y ) -2 -1 y -x D ( x , y ) -10 -5 0 5 10 y -x -6 -4 -2 D ( x , y ) (a) (b)(c) (d) FIG. 9: (a) Numerically computed kernel at x = 0 associated with D obtained fromlinear-algebra-based MFM for the described inhomogeneous problem. Different plotsrepresent mesh refinement: squares ( N = 40 , N = 10), circles ( N = 80 , N = 20), solidline ( N = 160 , N = 40). (b) eddy diffusivity associated with the homogeneous operator D = 1 / (cid:112) − ∂ /∂x , also plotted in semi-log scales (c) and (d). D ji ( x , y ), that expresses the macroscopic fluxes in terms of the macroscopic gradients as − u (cid:48) j c (cid:48) ( x ) = ˆ y D ji ( x , y ) ∂c∂x i | y dy . (50)We predict that a similar result holds for the case of transported vector fields with thedifference that the eddy-diffusivity operator will be described by a kernel involving a fourth-order tensor as − u (cid:48) j v (cid:48) i ( x ) = ˆ y D jilk ( x , y ) ∂v k ∂x l | y dy . (51)For the numerical example considered above, we show that the implied kernel from ourdiscrete solution, D ( x ,i , y ,j ) = [ D ] ij / ∆ x , converges by mesh refinement (here, the problemis 1D, and subscripts after the bracket imply discrete mesh points). As an example, we showin Figure 9.a the plot of D ( x = 0 , y ) versus y obtained after successive mesh refinements.31elow we report the major observations, and some speculations regarding their validityfor more general settings:1- For cases with zero normal velocity near a boundary and with boundary conditions thatremain invariant under averaging, the eddy-diffusivity kernel vanishes near the boundary(see Figure 8.b). Therefore, the macroscopic operator reduces to the molecular diffusivityoperator near a boundary. This implies there is no need for additional boundary conditionswhen scaling up to the macroscopic fields.2- In the present example, the macroscopic space was 1D and the eddy-diffusivity involveda positive kernel. For multidimensional systems the kernel is tensorial and entries can bepositive or negative. It will be interesting to examine however, whether the diagonal entriesof the tensor are positive-valued, and whether the tensor is positive definite.3- In the present example, the computed [ D ] was symmetric to the machine precision. Thisimplies that if a mean gradient at x A result in a mean flux at x B , the same mean gradientapplied to x B will result in same mean flux at x A . In other words, there is reciprocity ineddy diffusivity operator. It will be interesting to investigate under what conditions thisreciprocity holds, and whether it is generalizable to cases with mean advection as well ascases with tensorial eddy diffusivity.4- When the macroscopic dimensions are statistically homogeneous the eddy-diffusivitykernel can be simplified as D ( x , y ) = D ( x − y ), and thus convolution in the physical spacewill result in multiplication in the Fourier space. The specific fitted form for the eddy-diffusivity operator identified as D = D/ (cid:112) I − l ∂ /∂x (see Equations 28 and 40) resultsin a qualitatively similar kernel as that obtained for the inhomogeneous problem (see Fig-ure 9.b). Specifically, the kernel is positive for all values of y and the kernel width is on theorder of l , which itself scales with the characteristic eddy size.5- In the limit that the eddy size (i.e. kernel width) is much smaller than the macroscopiclength, the eddy-diffusivity kernel can be safely approximated by a Dirac delta function. Inother words, the convolution can be approximated as multiplication by the area under thekernel. This is what Taylor’s model, presented in Section I B, achieves. Locality of thekernel width is also one of the key conditions necessary for the validity of the Boussinesqapproximation. However, even in this limit, the eddy diffusivity is not necessarily an isotropictensor further challenging the validity of the Boussinesq approximation.6- In earlier sections we saw universally in all examples that in the limit of large wavenum-32er, the macroscopic closure operator is proportional to k , implying that (cid:98) D ∼ k − . Inversetransforming this relation to physical space implies that the eddy-diffusivity kernel, D ( x , y ),must be proportional to log (1 / | x − y | ) for small separation distances in macroscopic 1Dsettings as shown in Figure 9.c. Consistently D ( x , y ) ∼ / | x − y | in macroscopic 2D and D ( x , y ) ∼ / | x − y | in macroscopically 3D settings. These are all singular but integrablekernels.7- In earlier sections we also observed that the eddy-diffusivity operator is a smooth curvein Fourier space. We here assume an infinitely differentiable curve in k -space. This impliesthat the kernel associated with the eddy diffusivity, D ( x , y ), must vanish exponentially inthe limit of large | x − y | as shown in Figure 9.d. B. Inverse MFM
The insights gained from examination of the full macroscopic operator in Section IV Aallows development of a more economical method for approximate determination of themacroscopic 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 c in response to a specific macroscopic forcing. In contrast, IMFMobtains the macroscopic forcing required for sustaining a pre-specified c . For problems inwhich the microscopic PDE involves a time stepping process, IMFM is straightforward.Given the PDE ∂c∂t = f ( c, ∇ c, t, ... ) + s, (52)one simply can take the average of this equation to arrive at ∂c∂t = f + s. (53)By expanding the time discrete form of the left hand side term, and substituting c ( t + ∆ t )with the pre-specified c one obtains an explicit expression for s = s which can be sub-stituted in the equation. For example, for the forward Euler time-stepping method s =[ c ( t + ∆ t ) − c ( t )] / ∆ t − f ( t ). Substituting this expression in the microscopic PDE con-strains c to be the pre-specified value, but the non-macroscopic modes in c are allowed toevolve. 33or problems with a steady microscopic PDE, this treatment can be used in conjunctionwith artificial time stepping, where the true s is evaluated in the long term limit. Forturbulent flow problems that are statistically time stationary, the microscopic PDE is stilltime dependent. In such cases, s is time-independent but is needed instantaneously. Onemay 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 runsimulations 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 andconverge quickly with number of ensembles in our MFM analysis of turbulent channel flows(not presented here). While running simultaneous ensembles may increase computationalcost per time step, it does not affect the total cost since convergence can be achieved overa proportionally lower number of time-steps.Lastly we note that when performing IMFM, the initial condition and boundary condi-tions must be adjusted to be consistent with the input c . For example, while the left-rightboundary conditions in the example discussed in Section IV A remain of Dirichlet type, thespecified values of c at the boundary must be consistent with those of c evaluated at theboundary.The introduced IMFM methodology can be used as another brute force method to revealcolumns of matrices associated with closure operators, or even in Fourier space for homoge-neous problems. However, this method can also be used to reveal the low-order moments ofthe eddy-diffusivity with few DNS simulations as discussed below. We present this methodspecifically because it can offer a cost-effective MFM.With the IMFM methodology described above one may use a small number of simulationsto determine the low-order moments of the eddy-diffusivity kernel, D ( x , y ). For example,for the macroscopically 1D problem discussed in Section IV A, by selecting c = x andperforming one IMFM DNS, and post-processing the u (cid:48) c (cid:48) data, one finds the zeroth momentof the kernel as (see Equation 49), D ( x ) = ˆ y D ( x , y ) dy = − u (cid:48) c (cid:48) | c = x . (54)As we shall see, D is the coefficient of the leading-order macroscopic closure operator,which is the usual diffusion operator. In the limit that the “mean free path” or eddy sizeis much smaller than the macroscopic scale this will be the dominant operator. For multi-34imensional problems, D is replaced by a tensor due to its anisotropy, i.e. the Boussinesqlaw is not followed.To quantify the non-local effects, higher order moments must be computed. The firstmoment of the eddy-diffusivity kernel can be computed by selecting c = x / D ( x ) = ˆ y ( y − x ) D ( x , y ) dy = − (cid:16) u (cid:48) c (cid:48) | c = x / − x u (cid:48) c (cid:48) | c = x (cid:17) , (55)where the last term on the right hand side can be substituted using D . Likewise, the secondmoment of D can be computed as D ( x ) = ˆ y
12 ( y − x ) D ( x , y ) dy = − (cid:16) u (cid:48) c (cid:48) | c = x / − x u (cid:48) c (cid:48) | c = x / + (cid:0) x / (cid:1) u (cid:48) c (cid:48) | c = x (cid:17) . (56)One may generalize this recursive methodology to compute higher moments of D . Deter-mination of the leading-order moments requires few IMFM DNSs, and thus it is much lesscostly than the brute force MFM.The additional advantage of IMFM, either in the context of moments or even when usedas a brute force method, is that it allows direct probing of the desired unclosed terms in themacroscopic operator. For example, in this case, one may directly determine the operatorassociated with the unclosed flux u (cid:48) c (cid:48) , without requiring first obtaining the full macroscopicoperator L (as done in FMFM). Aside from cost saving, closure of individual fluxes, as op-posed to the full unclosed term, has the advantage of guaranteeing conservativeness of themacroscopic model. For more complicated microscopic systems that may involve multipleunclosed terms, IMFM allows probing of different unclosed terms independently, and thusprovides more insightful guidelines for modeling and understanding of the underlying dy-namics. Note that in IMFM, the forcing field s is not explicitly used when examining themacroscopic operator.Next, we discuss how the macroscopic operators can be constructed with reasonableaccuracy from the knowledge of these moments. C. Construction of the Macroscopic Operators from Kernel Moments
Having moments of D at hand, the task is to approximate the convolution integralsin (50). We first consider a non-convergent, but conceptually insightful, approximation.35 . Approximation of kernel integral using Taylor series The integral in Equation (50) can be rewritten by expanding ∂c/∂x i in terms of itsspatial Taylor series around y = x . Substitution and utilization of the defined momentsleads to representations of u (cid:48) j c (cid:48) in terms of a series involving spatial derivatives of c with thecomputed moments as pre-factors to each term. For example, for the test case considered inSection IV A, the unclosed flux, and thus the macroscopic closure operator, can be modeledas − u (cid:48) c (cid:48) ( x ) = (cid:20) D ( x ) + D ( x ) ∂∂x + D ( x ) ∂ ∂x + ... (cid:21) ∂c∂x . (57)The leading term in this expansion, D , is the local approximation to the eddy-diffusivityoperator; it approximates the kernel with a Dirac delta function with matching integral.Higher order terms in (57) are approximations of the non-local effects based on local Taylorseries of c (!). These terms involve the higher moments of D which can become significantwhen the macroscopic scale becomes comparable to the kernel’s characteristic width, i.e.,eddy size.One might expect that a truncated expansion based on a few moments beyond D couldresult in considerable model improvement. Such expansion would still lead to a well-poseddifferential equation, since coefficients (to higher derivatives) vanish near the boundary as thekernel itself vanishes, and thus a model based on this expansion does not impose requirementof additional boundary conditions. However, the expansion presented in Equation (57) isnot useful in practice since addition of higher moments does not lead to convergence of thesolution. For example, considering the test case discussed in Section IV A, we use DNS ofthe following microscopic system to assess performance of macroscopic models, u ∂c∂x + u ∂c∂x = 0 . ∂ c∂x + ∂ c∂x + 1 . (58)Figure 10.a shows the 2D DNS results and Figure 10.b shows the macroscopic fields. Amacroscopic model that retains only the leading operator, D , performs a reasonable job inpredicting the qualitative features of the mean field, while still lacking quantitatively. Whenexpansion (57) is used, addition of the new terms did not lead to prediction improvementdespite relatively smooth mean fields.This observation is consistent with the issue seen in Section II A when we examinedhomogeneous parallel flow systems. In fact Taylor’s model represented the zeroth moment36 x x -3 -2 -1 0 1 2 3 x c _ FIG. 10: (a) contour plots representing concentration field computed from 2D DNS ofEquation (58). (b) concentration averaged in the x direction: solid, DNS; dashed, eddydiffusivity represented by Equation (57) but including only D ; dash-dotted, eddydiffusivity represented by Equation (57) including D and D terms; circle, eddy diffusivityobtained from the kernel reconstruction method.of the macroscopic diffusivity kernel for the considered flow while its perturbative correctionwas indeed an analytical derivation of the higher moments in the expansion of convolutionintegral (with the exception that in the analytical derivation in Section II A we consideredthe unsteady problem and thus captured both spatial and temporal moments). As previouslyshown in Figure 2.c, addition of terms in the moment expansion does not result in improvedmacroscopic prediction.It turns out, the expansion in (57) is analogous to the Kramers-Moyal expansion[36]previously formulated for the statistical description of stochastic processes associated withBrownian dynamics. Lack of convergence of this expansion has been previously investigatedas described by the Pawula theorem. Likewise, in these applications, while the leadingmoment promisingly captures a significant portion of the macroscopic behavior, additionof higher moments do not improve quantitative prediction. We are aware of one previouswork[26] reporting analogy to the Kramers-Moyal expansion in the continuum context intheir analytical study of laminar transport of reacting scalars in porous media.37 . Approximate kernel reconstruction Although Equation (57) is not a useful expansion in practice, the computed momentsfrom IMFM are still valuable in the sense that one can utilize them to construct convergingmacroscopic operators. We next show an example of such construction. Our approach hasbeen guided by the fact that all higher-order operators beyond D are inherently non-local.We therefore construct approximate convolution kernels that match the moments obtainedfrom MFM. The details of the procedure for the example problem in Figure 10 are as follows.We first used IMFM to compute the first three moments of D , i.e., D , D , D . We thenconstructed positively valued kernels that matched the computed moments at each x . Oneapproach would be to assume known kernel profiles, and tune their parameters to matchthe low-order moments (see Appendix E). We then used the constructed kernel to setup themacroscopic equation as, ∂∂x (cid:20) ˆ y D ( x , y ) (cid:18) ∂c∂x (cid:19) | y dy + 0 . ∂c∂x (cid:21) + 1 = 0 . This macroscopic equation can be solved numerically after discretization. For practicalproblems involving 3D convolution integrals, one may need to design iterative techniques forfast approximate calculation of the convolution integral. Figure 10.b shows the macroscopicsolution obtained through the reconstructed kernel method. Compared to prediction of themodel that retains only the leading operator, D , the solution from the reconstructed kernel,not only shows convergence, but also significant quantitative improvement.Due to the limited scope of this report, we skip providing examples of the applicationof IMFM for analysis of inhomogeneous and wall-bounded turbulent flows. However, thisis our main intention in developing the presented methodology. Extension to turbulentflows and its advantages are straightforward and clear. Such extensions will help turbulencemodels in two ways. First, by considering only linear mean momentum profiles in GMT,and performing IMFM, one can measure the full tensorial form of the leading-order eddydiffusivity D jilk ( x ) = − u (cid:48) j v (cid:48) i | v m = x n δ nl δ mk . (59)In other words, with only nine DNS simulations, the local (leading) term in the eddy-diffusivity operator can be measured precisely. One can then use this knowledge to quan-titatively assess the anisotropy in momentum mixing by the underlying turbulent flow. To38his end, IMFM acts as a non-intrusive rheometer determining the macroscopic “materialproperties” of the underlying turbulence.Second, by computing the higher moments of D one can quantitatively assess the levelof non-locality involved in turbulent mixing. Quantitative and independent assessments ofanisotropy and non-locality provide invaluable insights in determining the missing pieces intoday’s RANS models. V. CONCLUSIONS
We presented a systematic numerical procedure, which we call MFM, for determination ofscaled up operators describing the macroscopic influence of microscopic transport processes.MFM is applicable to both laminar and turbulent flows, homogeneous and inhomogeneousflows, and can be used to examine transport of scalar fields as well as vector fields such asthose describing momentum as a transported quantity.For a given setting, MFM can obtain the exact macroscopic differential operator gov-erning the evolution of mean fields. However, MFM does not solve the closure problem inturbulence modeling. This is because MFM obtains the operator coefficients explicitly interms of the independent coordinates for each specific setup. Tackling the turbulence closurequestion requires expressing these coefficients in terms of the mean field quantities them-selves resulting in nonlinear operators. The advantage of MFM is therefore in quantificationof the model-form error and identification of improved model forms. Whether the RANSmodel should involve an anisotropic but local eddy diffusivity or a non-local operator, thepresented methodology can provide this information quantitatively.Additionally, MFM provides a robust framework for RANS model verification. In today’spractice, RANS model forms are postulated based on physical intuition, and their coefficientsare 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 whatdegree 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 allowsa more rigorous verification of RANS models.We point out that MFM is an expensive procedure (either forward or inverse), involvingat least multiple DNS solutions. Therefore, we expect this method to be mostly used to39nvestigate canonical problems. The hope is that an understanding of closure operators incanonical flows will reveal some universality, and that this understanding can be extendableto flows with more complex geometries. As an example, we showed how application of MFMin homogeneous and isotropic turbulence can inform new RANS operators which then led tosignificant improvements in prediction of axisymmetric jet flows without requiring an MFManalysis of jets.Examination of the macroscopic operators for homogeneous flows in Fourier space re-vealed a smooth curve with k scaling in the low-wavenumber limit and k scaling in thehigh-wavenumber limit. The k scaling was universally observed in a wide range of condi-tions for both laminar and turbulent problems and for both scalar and momentum trans-port. While the k scaling is consistent with the intuitive expectation of bounded dispersionspeed, the pre-factor to this scaling is real and positive, and thus cannot be captured by anadvective operator. Instead, we showed that this scaling is indicative of a non-local, but dis-sipative mechanism of transport, as suggested by the fitted operator, ∇ (cid:2) D/ √ − l ∇ (cid:3) ∇ .Due to presence of two asymptotic scalings in the small and large k limits, the non-localitydoes not conform to a simple fractional-order differential operator expressed as a Laplacianto a non-integer power. Extending our observations to inhomogeneous cases, we concludethat 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 intro-duced kernel D ( x , y ) that expresses mean fluxes in terms of a weighted superposition of thesurrounding mean field gradients.We point out that MFM can also enrich LES modeling in the long term. One fundamentalapproach is to recast the LES modeling problem in the MFM context by modifying thedefinition of averages that are used here. A more immediately accessible approach is toutilize knowledge derived from MFM in the RANS context. For example, the reported k scaling of macroscopic operator for the high-wavenumber limit is also expected in the caseof LES, implying the LES eddy diffusivity should vanish in the high-wavenumber limit.Additionally, the required anisotropy treatments in wall-modeled LES can be informed fromanisotropy of D to be understood in the RANS context.While in this report we considered scalar and momentum transport in the context ofincompressible flows, it is possible to envision extension of MFM for broader applicationsincluding reacting scalars, compressible systems, shock-turbulence interactions, and buoy-40ncy driven flows. These are examples where anisotropy and/or coupling with other physicsis known to impose challenge in modeling. MFM is a valuable tool that can shed light intoeffects of multiphysics coupling on closure operators in a quantitative manner.MFM can be viewed as a numerical rheometer, similar to the way that molecular dynamicssimulations predict continuum-level transport coefficients. Experimental rheometers mea-sure momentum diffusivity due to the underlying chaotic Browning dynamics. When usingthem, the key fundamental assumption is that placing of the material under the rheome-ter 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. Inthis sense, a novel aspect of our work is recasting the Navier-Stokes equation in a moregeneral form by separating the roles of velocity as a transporter and transportee. Rigor-ous rheometry must be non-intrusive to the transporter mechanism but it can be intrusiveto the transportee field. Given emerging supercomputing power, MFM provides a freshopportunity to study turbulence by measuring its “material properties”. VI. ACKNOWLEDGEMENTS
This work was supported by the Boeing Company under grant number SPO 134136.and the National Science Foundation under award number 1553275. Danah Park was alsosupported by a Kwanjeong Educational Foundation Fellowship. The authors are gratefulfor critical comments from Drs. Chad Winkler, Philippe Spalart, and Mori Mani throughthe course of the research resulting in this report. The authors are in debt to Jessie Liu andDr. Karim Shariff for their technical discussions and elaborate comments on a draft of thisreport. The phrase “Generalized Momentum Transport” (GMT) as a name for Equation (33)was suggested in discussion with Dr. Karim Shariff.
Appendix A: Solution to Equation (21)
Substitution of Equation (22) into (21) results in a series involving cos ( nx ). The productcos ( nx ) cos( x ) can be expanded as (1 /
2) cos [( n + 1) x ] + (1 /
2) cos [( n − x ]. Using thisexpansion and balancing the left-hand side and right-hand side of the equation for each term41esults in iωa + ik a = 1 , (A1) ika + ( iω + 1) a + ik a = 0 , (A2) ik a n − + ( iω + n ) a n + ik a n +1 = 0 , (A3)where the last equation applies for n ≥
2. The recursive relation in (A3) needs two initialconditions, i.e., a and a . Depending on the ratio of a and a , it admits two modesof solutions, one increasing in magnitude exponentially, and one decreasing in magnitudeexponentially. To obtain a physical solution, we impose the condition that the exponentiallyincreasing mode must be zero. This condition sets a unique value for the ratio a /a . Weinvestigated this condition by rewriting (A3) in terms of r n = a n /a n − and re-arranging as r n +1 = 2 ( − n − iω ) ik − r n . The question is to determine the value of r such that this recursive relation leads to | r n →∞ | smaller than unity. A randomly selected r leads to unbounded | r n →∞ | . This can be un-derstood better when thinking in terms of Equation (A3). A random combination of a and a will have both exponentially increasing and decreasing modes present. Even with asmall pre-factor, the exponentially increasing mode eventually prevails for sufficiently large n . Only when the ratio of a and a are set to precisely remove the increasing mode, one willobtain the physically correct solution. Our numerical trick is to use the recursive relationbackward in n . By choosing a random finite number for r N where N is a sufficiently largenumber, we solve the above recursive relation backward in n to find a value for r . Movingbackward 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 to be converged within many significant digits. Once the ratio a /a is at hand, one cansolve Equations (A1) and (A2) to obtain a . Appendix B: Numerical Solution of Equation (24)
Even though in other sections we have solved macroscopic equations directly in physicalspace using techniques of linear algebra, for the specific problem discussed in Section II A, itwas easier to solve the problem in the Fourier space for both spatial and temporal dimensions.42o do this, we represent the initial condition as a source term that is activated at time t = 0via a Dirac delta function, s ( x , t ) = exp (cid:0) − x / . (cid:1) δ ( t ) . The above expression is Fourier transformed in space and time using a physical mesh with∆ x = 0 .
05 and ∆ t = 0 .
025 and a truncated domain − ≤ x ≤ (cid:98) s ( k, ω ) is then dividedby (cid:98) L ( k, ω ) = (cid:113) (1 + iω ) + k − c ( x , t ) in thephysical space.Next, we discuss the treatment of mode k = ω = 0. Given that the temporal domainis not infinite, this procedure results in a time-periodic solution corresponding to repeatedspontaneous injections of scalar in the domain. To eliminate the artifacts associated withthis periodicity, we applied the following remedies. First, the temporal domain is taken tobe substantially large (0 ≤ t ≤ c to avoid a linear growth due to the periodic cycles. To this end, the mode k = 0 is solvedanalytically in the physical space using the trivial differential equation ∂/∂t = 0.To verify our method, we repeat this procedure using the exact (cid:98) L , obtained throughMFM discussed in Appendix A and confirm that the obtained c ( x , t ) matches that of thetwo-dimensional DNS solved in the physical space under resolution refinement. Appendix C: Macroscopic Operator for a Solenoidal Passive Vector Field Subjectto a Parallel Flow
As discussed in Section II C, we consider the simple limit of zero streamwise microscopicdiffusivity ( (cid:15) = 0), and steady limit of the microscopic equations. By considering s =exp ( ikx ), v i = (cid:98) v i ( x ) exp ( ikx ) and q = (cid:98) q ( x ) exp ( ikx ), we seek the solution in Fourierspace. Substituting these expressions in (29) and (30) and simplification results in ik cos ( x ) (cid:98) v = d (cid:98) v dx − ik (cid:98) q, (C1) ik cos ( x ) (cid:98) v = d (cid:98) v dx − (cid:98) q (cid:48) + 1 , (C2)43 k (cid:98) v + d (cid:98) v dx = 0 . (C3)One may eliminate variables by combining these equations and obtain a single equation forthe unknown (cid:98) v . ik cos ( x ) (cid:98) v + 1 k d (cid:98) v dx + 1 ik ddx (cid:20) cos ( x ) d (cid:98) v dx (cid:21) − d (cid:98) v dx = 1 . (C4)Next, we consider the following series solution for (cid:98) v (cid:98) v = ∞ (cid:88) n =0 a n cos ( nx ) . (C5)Substitution in (C4) and rearranging leads to ika = 2 (C6)2 ik a + 2 (cid:0) k + 1 (cid:1) a + ik (cid:0) k + 2 (cid:1) a = 0 (C7) ik (cid:0) k + n ( n − (cid:1) a n − + 2 (cid:0) n k + n (cid:1) a n + ik (cid:0) k + n ( n + 1) (cid:1) a n +1 = 0 , n ≥ . (C8)The solution methodology is similar to what we discussed in Appendix A since the systemof equations are similar to those obtained for the scalar field. Once this solution at hand,one can obtain the macroscopic closure operator as (cid:98) L (cid:48) ( k ) = (cid:98) s / (cid:98) v = 1 /a ( k ). Appendix D: Details of RANS Solutions for the Axisymmetric Jet
For this problem we follow the solution of Tollmien[37] by seeking a similarity solutionof the form u ( x, r ) = 1 x f ( η ) = 1 x f (cid:16) rx (cid:17) , (D1)where u is the axial velocity and x and r are the axial and radial coordinates respectively.Following approximations for semi-parallel flows, the RANS equation can be written as: u ∂u∂x + v ∂u∂r = 1 r ∂∂r (cid:0) − ru (cid:48) v (cid:48) (cid:1) , (D2)where v = (1 /x ) g ( η ) is the mean radial velocity and can be computed from f via integralof the continuity equation as η g ( η ) = ˆ η η ddη [ ηf ( η )] dη. (D3)44n Prandtl Mixing Length Model (PMLM), the following closure model is used − u (cid:48) v (cid:48) = D ∂u∂r , (D4)where the eddy diffusivity, D , is approximated as the product of a mixing length and avelocity scale, D = l p u p . The Prandtl mixing length, l p , is assumed to scale with the jet halfwidth l p = α r = α xη , where η is defined such that f ( η ) = f (0) /
2, and α is an order unity constant to bedetermined. The Prandtl velocity scale, u p , is assumed to scale with the product of themixing length and mean velocity gradient as u p = α l p | ∂u/∂r | . With this description, the closure model becomes − u (cid:48) v (cid:48) = α α r (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂r (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂r . Substitution of this model in (D2) results in the following similarity ODE − f ddη ( ηf ) + g dfdη = 1 η ddη (cid:20) α α η η (cid:12)(cid:12)(cid:12)(cid:12) ∂f∂η (cid:12)(cid:12)(cid:12)(cid:12) ∂f∂η (cid:21) . (D5)In PMLM, it is only necessary to prescribe the combined constant α α ; the model remainsambiguous about the independent values of α and α .In the MFM-inspired closure model, we replace the eddy diffusivity in (D4) with anoperator similar to that in Equation (40). Following simplifications for semi-parallel flows,this leads to the following closure model − u (cid:48) v (cid:48) = (cid:26) I − l p (cid:20) r ∂∂r (cid:18) r ∂∂r (cid:19) − r (cid:21)(cid:27) − D ∂u∂r . (D6)Given both l p in PMLM and l in MFM-based operators scale with large eddy size (alsoquantified as u rms /ε ), and given α is still a free parameter, we can simplify the analysisby defining l p = l . Substituting this model into (D2) and rearranging terms results in thefollowing similarity ODE − f ddη ( ηf ) + g dfdη = 1 η ddη (cid:32) η (cid:26) I − α η (cid:20) η ∂∂η (cid:18) η ∂∂η (cid:19) − r (cid:21)(cid:27) − (cid:20) α α η (cid:12)(cid:12)(cid:12)(cid:12) ∂f∂η (cid:12)(cid:12)(cid:12)(cid:12) ∂f∂η (cid:21)(cid:33) . (D7)45e solved both equations (D5) and (D7) each coupled with (D3), to respectively repre-sents solutions of the original and MFM-inspired RANS models. The numerical procedureused a second-order central difference scheme. The square root and inverse operators in (D7)were computed using MATLAB’s commands after casting the operators in matrix format.The computational domain was 0 ≤ η ≤ .
9. At η = 0 a normalization condition f = 1 wasused, while η = 0 . f = 0 would be justifiable. An iterative procedure was designed toconverge to the steady solution.The remaining task is determination of coefficients α and α . In PMLM the combinedcoefficient α α is the only coefficient that must be determined, but this coefficient is tunedrather than being derived. In our analysis of this model we selected α α = 0 .
05 so that thecomputed η = 0 . α (instead oftuning).Given the MFM analysis of HIT, presented in Section III B, and reported in more detailin [33], l p must be l p = l = Cu /ε, (D8)where C is a constant and we seek its value in the limit of high Re. Shirian and Manireported C to be 0.69 for Re λ = 26 and 0.67 for Re λ = 40 (see their Table II). We anticipate C to be a bit smaller than 0.67 in the much higher Re limit, but given all uncertainties, wecontinue our analysis by assuming C = 0 . u rms and ε in (D8) using reported experimental measurements ofthe jet centerline. A recent high-quality data is provided by Darisse et al. [38]. While thisstudy reports a mean centerline velocity profile consistent to that of Hussein [35], which weused in Figure 6, they provide more accuracy in reporting of kinetic energy and dissipationrates. From this study we estimate u rms = 0 .
22 (see Figure 5 in [38]), and ε = 0 .
017 (seeFigure 13 in [38]), both in units of centerline velocity and r . Substituting these resultsin (D8) yields, l p = 0 . r . Hence, α = 0 .
38 or equivalently α η = 0 . α α such that the computed η = 0 . α α = 0 . D near the axis of symmetry. Appendix E: Kernel Reconstruction from Moments
For the problem considered in Section IV C 2, we considered an approximate kernel de-scribed by D ( x , y ) = h ( x ) exp (cid:16) x − y l ( x ) (cid:17) , x < y ; h ( x ) exp (cid:16) y − x l ( x ) (cid:17) , x ≥ y . (E1)For each x -position, the moments of the above kernel form can be computed analyticallyin terms of h , l and l . By matching these moments to those obtained from IMFM, weobtained profiles of h , l and l versus x . Aside from kernel shape approximation anddiscretization errors, we commit two additional errors when matching the moments. First,the analytical computation of the moments in terms of h , l and l does not consider domaintruncation by walls. In other words, we assume that l , l , and h become sufficiently smallnear the wall to allow neglecting of this error. Second, in rare points where the matchedvalue of l or l becomes negative, we replace its value with zero. By utilizing such kernelswe show in Figure 10 significant improvements in prediction of mean concentration fieldin comparison to that obtained from the local eddy diffusivity. This kernel reconstructionmethod also resolves the convergence issue observed when corrections to the leading operatorwas incorporated using expansions based on Taylor series and moments.For the specific problem considered here, we avoided using the canonical homogeneouskernel, D = D/ (cid:112) I − l ∂ /∂x because this kernel is symmetric in space, while the actualkernel involves spatial asymmetry near the boundaries. However, we here briefly note that47n a homogeneous system the moments of this kernel are D = D , D = 0, and D = Dl / [1] Walter Guido Vincenti and Charles H Kruger, Introduction to physical gas dynamics , Vol. 246(Wiley New York, 1965).[2] Parviz Moin and John Kim, “Tackling turbulence with supercomputers,” Scientific American , 62–68 (1997).[3] Stephen B Pope, “Turbulent flows,” (2001).[4] Charles G Speziale, “Analytical methods for the development of reynolds-stress closures inturbulence,” Annual review of fluid mechanics , 107–157 (1991).[5] C Speziale, “Turbulence modeling for time-dependent rans and vles: a review,” AIAA journal , 173–184 (1998).[6] Paul A Durbin, “Some recent developments in turbulence closure modeling,” Annual Reviewof Fluid Mechanics , 77–103 (2018).[7] PRaA Spalart and S1 Allmaras, “A one-equation turbulence model for aerodynamic flows,”in (1992) p. 439.[8] K Hanjali´c and BE Launder, “A reynolds stress model of turbulence and its application tothin shear flows,” Journal of fluid Mechanics , 609–638 (1972).[9] David C Wilcox, “Formulation of the k-omega turbulence model revisited,” AIAA journal ,2823–2838 (2008).[10] Florian R Menter, “Two-equation eddy-viscosity turbulence models for engineering applica-tions,” AIAA journal , 1598–1605 (1994).[11] Paul A Durbin, “Separated flow computations with the k-epsilon-v-squared model,” AIAAjournal , 659–664 (1995).[12] J Boussinesq, “Essai sur la th´eorie des eaux courantes,” in M´emoires pr´esent´es par diverssavants `a l’Acad´emie des Sciences (Impr. nationale, 1877).[13] SB Pope, “A more general effective-viscosity hypothesis,” Journal of Fluid Mechanics ,331–340 (1975).[14] Brian Edward Launder, G Jr Reece, and W Rodi, “Progress in the development of a reynolds-stress turbulence closure,” Journal of fluid mechanics , 537–566 (1975).[15] John Kim, Parviz Moin, and Robert Moser, “Turbulence statistics in fully developed channel ow at low reynolds number,” Journal of fluid mechanics , 133–166 (1987).[16] Geoffrey Ingram Taylor, “Dispersion of soluble matter in solvent flowing slowly through atube,” Proc. R. Soc. Lond. A , 186–203 (1953).[17] Osborne Reynolds, “On the dynamical theory of incompressible viscous fluids and the deter-mination of the criterion,” Philosophical Transactions of the Royal Society of London. A ,123–164 (1895).[18] Andrew J Majda and Peter R Kramer, “Simplified models for turbulent diffusion: theory,numerical modelling, and physical phenomena,” Physics reports , 237–574 (1999).[19] R Aris, “On the dispersion of a solute in pulsating flow through a tube,” Proc. R. Soc. Lond.A , 370–376 (1960).[20] Armand Ajdari, Nathalie Bontoux, and Howard A Stone, “Hydrodynamic dispersion in shal-low microchannels: the effect of cross-sectional shape,” Analytical Chemistry , 387–392(2006).[21] NG Barton, “On the method of moments for solute dispersion,” Journal of Fluid Mechanics , 205–218 (1983).[22] Søren Vedel and Henrik Bruus, “Transient taylor–aris dispersion for time-dependent flows instraight channels,” Journal of Fluid Mechanics , 95–122 (2012).[23] WN Gill and R Sankarasubramanian, “Exact analysis of unsteady convective diffusion,” Proc.R. Soc. Lond. A , 341–350 (1970).[24] Donald L Koch and John F Brady, “A non-local description of advection-diffusion with ap-plication to dispersion in porous media,” Journal of Fluid Mechanics , 387–403 (1987).[25] I Frankel and H Brenner, “On the foundations of generalized taylor dispersion theory,” Journalof Fluid Mechanics , 97–119 (1989).[26] Roberto Mauri, “Dispersion, convection, and reaction in porous media,” Physics of Fluids A:Fluid Dynamics , 743–756 (1991).[27] Wen Chen, “A speculative study of 2/ 3-order fractional laplacian modeling of turbulence:Some thoughts and conjectures,” Chaos: An Interdisciplinary Journal of Nonlinear Science , 023126 (2006).[28] Brenden P Epps and Benoit Cushman-Roisin, “Turbulence modeling via the fractional lapla-cian,” arXiv preprint arXiv:1803.05286 (2018).[29] Andriy Yaroshchuk, Emiliy Zholkovskiy, Sergey Pogodin, and Vladimir Baulin, “Coupled oncentration polarization and electroosmotic circulation near micro/nanointerfaces: Taylor–aris model of hydrodynamic dispersion and limits of its applicability,” Langmuir , 11710–11721 (2011).[30] I Rubinstein and B Zaltzman, “Convective diffusive mixing in concentration polarization: fromtaylor dispersion to surface convection,” Journal of Fluid Mechanics , 239–278 (2013).[31] P Go Saffman, “On the fine-scale structure of vector fields convected by a turbulent fluid,”Journal of Fluid Mechanics , 545–572 (1963).[32] HK Moffatt, “Some developments in the theory of turbulence,” Journal of Fluid Mechanics , 27–47 (1981).[33] Mani A. Shirian Y., “Eddy diffusivity in homogeneous and isotropic turbulence,” In review(2019).[34] L Prandtl, “Report on investigation of developed turbulence,” Mechanik (1925).[35] Hussein J Hussein, Steven P Capp, and William K George, “Velocity measurements in ahigh-reynolds-number, momentum-conserving, axisymmetric, turbulent jet,” Journal of FluidMechanics , 31–75 (1994).[36] NG Van Kampen, “Stochastic processes in physics and chemistry,” (1981).[37] Walter Tollmien, “Berechnung turbulenter ausbreitungsvorg¨ange,” ZAMM-Journal of AppliedMathematics and Mechanics/Zeitschrift f¨ur Angewandte Mathematik und Mechanik , 468–478 (1926).[38] A Darisse, J Lemay, and A Benaissa, “Budgets of turbulent kinetic energy, reynolds stresses,variance of temperature fluctuations and turbulent heat fluxes in a round jet,” Journal ofFluid Mechanics , 95–142 (2015)., 95–142 (2015).