Reduced Numerical Modeling of Turbulent Flow with Fully Resolved Time Advancement. Part 1. Theory and Physical Interpretation

Simple Summary

It is too costly to compute the time advancement of many turbulent flows because the range of relevant motions is very broad and encompassing all of them is unaffordable. The typical compromise is to omit the smallest motions, replacing them by estimates of their effects, but this is not always adequate. Instead, it is proposed to fully resolve the flow advancement along selected line segments within the flow volume and to couple these segments so that the resulting formulation adequately represents all scales of motion. The modeling needed to implement this is described and the overall approach, termed autonomous microstructure evolution, is assessed conceptually and with reference to demonstrated capabilities of related methods.

Abstract

A multiscale modeling concept for numerical simulation of multiphysics turbulent flow utilizing map-based advection is described. The approach is outlined with emphasis on its theoretical foundations and physical interpretations in order to establish the context for subsequent presentation of the associated numerical algorithms and the results of validation studies. The model formulation is a synthesis of existing methods, modified and extended in order to obtain a qualitatively new capability. The salient feature of the approach is that time advancement of the flow is fully resolved both spatially and temporally, albeit with modeled advancement processes restricted to one spatial dimension. This one-dimensional advancement is the basis of a bottom-up modeling approach in which three-dimensional space is discretized into under-resolved mesh cells, each of which contains an instantiation of the modeled one-dimensional advancement. Filtering is performed only to provide inputs to a pressure correction that enforces continuity and to obtain mesh-scale-filtered outputs if desired. The one-dimensional advancement, the pressure correction, and coupling of one-dimensional instantiations using a Lagrangian implementation of mesh-resolved volume fluxes is sufficient to advance the three-dimensional flow without time advancing coarse-grained equations, a feature that motivates the designation of the approach as autonomous microscale evolution (AME). In this sense, the one-dimensional treatment is not a closure because there are no unclosed terms to evaluate. However, the approach is additionally suitable for use as a subgrid-scale closure of existing large-eddy-simulation methods. The potential capabilities and limitations of both of these implementations of the approach are assessed conceptually and with reference to demonstrated capabilities of related methods.

1. Introduction

To extend the range of applicability of numerical simulation of turbulent flows beyond the limits of affordable direct numerical simulation (DNS), under-resolved simulations with physically or algorithmically based subgrid-scale (SGS) closures are widely used, typically within the large-eddy-simulation (LES) framework. For constant-property flow, the main role of the SGS closure of LES is to dissipate kinetic energy at the physically correct rate based on the flow state. A secondary process that is important in some circumstances is the transfer of kinetic energy from unresolved scales to the LES-resolved flow solution, and accordingly, modeling of this backscatter has been incorporated into SGS closures. The resulting two-way coupling is especially important when LES is extended from constant-property flow to multiphysics flow involving combinations of buoyancy, multiphase couplings, and velocity divergence due to, for example, thermochemical processes.

Here, an approach designed to capture these and other multiphysics effects within SGS closures for LES is described. Unlike typical parameterizations, it involves idealized but physically based treatments of SGS advective and diffusive transport and other SGS processes, such as those mentioned. The approach involves DNS-level resolution of spatial structure and flow unsteadiness, which increases its cost relative to parameterizations, but it is less costly than DNS because this fine resolution is restricted to representative lines of sight within individual LES control volumes (CVs). Cost is further mitigated to the extent that the fidelity of the SGS closure allows coarsening of the LES mesh while maintaining the needed level of accuracy.

The method used for SGS closure involves one-dimensional turbulence (ODT), a form of map-based stochastic advection that has been used previously for SGS closure of LES as well as in various standalone applications, as cited below. Here, an ODT-based SGS closure is formulated that draws upon previously reported formulations but also introduces various novel elements, whose details and physical interpretations are the main focus of the present contribution.

To begin, it is important to clarify what it means to close an LES using ODT. Suppose that exact underlying DNS-level data are available at all times in all LES CVs. Then the unclosed terms in the LES equations can be evaluated exactly, yielding a fully accurate closure. However, if the LES time advances an SGS kinetic-energy equation, then the DNS data might be used to close an unclosed term in that equation. This might not be an optimal use of the DNS data, but if the available data were inexact rather than exact, then they might turn out to be more useful for closing a model equation carried by the LES rather than for closing the unclosed terms obtained by filtering the exact equation. For reacting flows, a related example is use of the data either to close the filtered chemical state or to provide the LES with reaction rates that are used to time advance filtered, hence modeled, chemical-kinetic equations.

ODT, in this context, provides surrogate DNS data that can be used to close either the filtered, but not otherwise modeled, LES equations or unclosed terms in submodels appended to these equations. Development and demonstration of either type of closure involves many choices with regard to mesh geometry, filtering technique, numerical algorithms, application cases, etc., ultimately resulting in the evaluation of only one or a few of the many options.

For the present purpose of introducing the general concepts and physical modeling approaches on which any particular modeling instantiation would be based, the following approach is adopted here. It is shown how ODT can be formulated so that, with suitable large-scale input, its time advancement can provide a local flow solution that, when filtered, is a physically plausible surrogate of the updated local filtered velocity that would otherwise be evaluated by advancing the filtered momentum equation. Thus, as an alternative to advancement of that equation, the updated filtered flow field can be obtained by time advancing and then filtering the fully resolved flow state on each ODT domain. There is one ODT domain per CV of the coarse-grained mesh, so this procedure updates the entire filtered flow field.

With the introduction of a pressure projection that adjusts the filtered flow field so as to enforce continuity, this approach becomes a self-contained ODT-based low-Mach-number flow simulation. It is a bottom-up approach that treats the evolving ODT state as the physical flow solution, albeit modeled rather than exact, that can be filtered to obtain an LES-like coarse-grained velocity field. In this context, the coarse-grained pressure is an auxiliary field that communicates large-scale information to each ODT domain, where this information reflects the collective interactions among the ODT domains and the inflow, outflow, and boundary conditions. Owing to the absence of any time advancement of the filtered velocity field, this formulation is not LES per se, so it is more appropriately termed autonomous microscale evolution (AME).

It is shown here that AME is a potentially viable alternative to LES. Its capabilities and limitations in this regard are assessed, subject to more definitive future evaluation by means of numerical implementation and model validation studies. As noted, the assessment also addresses points relevant to ODT SGS closure of LES. In what follows, operations on filtered fields are generically termed LES processes in order to discuss them in a familiar context until sufficient information specific to AME is provided so that discussion of AME becomes meaningful.

2. Overview of the Approach

2.1. ODT

In this section, the proposed SGS formulation of ODT is outlined to a sufficient extent that its coupling to other parts of the time-advancement cycle can be explained. For this purpose, the formulation is specialized to low-Mach-number constant-property flow. The filtered mesh is assumed to be Cartesian and all vector components are expressed in this fixed Cartesian frame.

Constant-property ODT time advances profiles of Cartesian velocity components $u_i$, $i=1,2,3$, on a 1D domain. This can be viewed as simulated flow evolution along a representative line of sight in a turbulent flow. Physically, such a line of sight is an open system that is continuously refreshed by fluid transport across the line of sight, but here as in most ODT formulations, the 1D evolution per se is a closed system. Other operations during the envisioned time-advancement cycle (see Section 2.3) capture large-scale and multidimensional influences that are not intrinsically represented by the evolution of the flow state on the ODT domain.

The velocity profiles $u_i$ are distinct from the velocity field ${\overline{v}}_i$ that governs coarse-grained advection in the LES sense. Here, the latter is the set of direction-i LES-prescribed CV face-normal velocities, as in a typical staggered numerical scheme. The overbar notation indicates box filtering. The box-filtering interpretation of ${\overline{v}}_i$ is explained in Section 3.

The ODT domains are Lagrangian objects, whose orientations need to be defined only for particular applications such as buoyant stratified flow and particle-laden flow. For those cases, the domain within a given CV is assumed to be aligned with the velocity vector $\overline{\mathbf{u}}$ obtained by box filtering the $u_i$ profiles, hence the domain orientation varies with time. (The unit vector in the $\overline{\mathbf{u}}$ direction will be denoted $\widehat{\overline{\mathbf{u}}}$). As noted, the components i of all vectors are referenced to the fixed Cartesian coordinates notwithstanding the varying domain orientation. A 2D visual rendering of a CV and its associated ODT domain is shown in Figure 1.

The ODT state is subject to two advancement processes. One is viscous advancement governed by

$$u_{i,t}=\nu u_{i,xx}-(1/\rho )p_{,i}.$$

(1)

In the last term, $\rho$ is density and $p_{,i}$ is the index-i directional derivative of the LES-prescribed pressure, deemed to be spatially uniform within the ODT domain. In the other terms, derivatives after a comma are denoted t for time or x for spatial location on the ODT domain. This spatial treatent is not equivalent to the 1D projection of the Laplacian of the velocity vector onto the domain direction, which cannot be evaluated in ODT owing to its reduced dimensionality. Rather, viscous transport of the vector velocity is modeled as the viscous transport of the velocity components treated as three scalar fields, which qualitatively represents the vector process.

The other advancement process is a stochastic sequence of instantaneous eddy events that punctuate the viscous advancement. Each eddy event modifies the system state within some interval $[x_0,x_0+l]$ of the ODT domain coordinate x, so the subsequent resumption of viscous advancement starts from the new system state. For present purposes, no-flux conditions are applied at the domain endpoints.

An eddy event consists of two operations. One is a rearrangement of all property profiles within the eddy interval, termed the triplet map. It is defined on the 1D continuum coordinate x (see Section 6.1), but it is intuitive to interpret it as the continuum limit of a permutation of the cells of a uniform discretization of the x coordinate, where all cell properties are carried with the cell upon its displacement. (This is one, but not the only or preferred, numerical implementation of the map.) Within ODT, this map is a Lagrangian model analog of the $\mathbf{u}·\nabla$ advection operation of the Eulerian momentum equation.

In addition to advective and viscous transport, momentum advancement includes momentum sources and sinks, at a minimum the pressure-gradient field that enforces continuity in the flow regime of present interest. In the present formulation, these are introduced in two ways. For the case of laminar flow, hence no eddies, the model reduces to viscous transport and the LES-prescribed pressure gradient as shown in Equation (1). From the ODT perspective, that term is an external source (or sink) of momentum and kinetic energy. The eddies that introduce turbulence effects imply concomitant ODT-level pressure fluctuations. The ODT eddy representation of the effects of these fluctuations is self-contained the simplest situations, so their modeled effects should not change the domain-integrated momentum or kinetic energy. A representation of pressure fluctuations that obeys this constraint is introduced as an additional eddy operation after implementation of a triplet map. It is formulated concisely in terms of chosen functions $K_i^{*}\left(x\right)$ that are added to the respective velocity components $u_i\left(x\right)$.

Momentum conservation must be obeyed by each individual velocity component $u_i\left(x\right)$ but kinetic-energy conservation applies only to the sum of component kinetic energies. For the constant-density regime under consideration, these constraints require $\int K_i^{*}\left(x\right)dx=0$ and ${\sum }_{i=1}^3\int {(K_i^{*}\left(x\right)+u_i\left(x\right))}^{2}dx={\sum }_{i=1}^3\int u_i{\left(x\right)}^{2}dx$. Other physically desirable properties are that $K_i^{*}\left(x\right)$ is continuous and is nonzero only within the eddy range $[x_0,x_0+l]$, where the latter reflects the association of the pressure fluctuation with the eddy motion. Additionally, requiring $K_i^{*}\left(x\right)$ to be piecewise linear has proven to be mathematically convenient.

On this basis, the formulation $K_i^{*}\left(x\right)=c_{i}K\left(x\right)$ is adopted, where $K\left(x^{\prime }\right)=x^{\prime }-x\left(x^{\prime }\right)$ is the map-induced displacement of the location x that is mapped to $x^{\prime }$, here expressed as the inverse-map function $x\left(x^{\prime }\right)$. For constant density, this conserves momentum by construction while kinetic-energy conservation imposes one constraint on the coefficients $c_i$.

Anticipating extension of this formulation to regimes such as variable-density flow, an additional function $J=\left|K\right|$ is introduced. Then the complete eddy event is the transformation

$$u_i\left(x\right)\to u_i^{\prime }\left(x\right)+c_{i}K\left(x\right)+b_{i}J\left(x\right).$$

(2)

Here, it is convenient to represent the mapped component profiles as new functions of the fixed coordinate x rather than as the original functions with transformed arguments, so they are denoted $u_i^{\prime }\left(x\right)$. The coefficients $b_i$ and $c_i$ are chosen so as to change component momenta and kinetic energies subject to applicable conservation laws and additonal physical modeling sufficient to specify all the coefficient values uniquely. The evaluation of $b_i$ and $c_i$ on this basis is explained in Section 6.1 and Section 6.2.

The main physical content of ODT is the stochastic process by which eddy intervals and times of occurrence are determined. The profiles $u_i$ of the velocity components within the various possible eddy intervals are the inputs to that determination, and the kernels play a key role in the specification of the eddy rate distribution $\lambda (x_0,l;t)$ based on that input. $\lambda$ has units time${}^{-1}$ length${}^{-2}$ such that $\lambda (x_0,l;t)dxdl$ is the rate of occurrence of eddies, whose left boundaries are within the range $[x_0,x_0+dx]$ and whose lengths are within the range $[l,l+dl]$, where the argument t reflects the dependence of $\lambda$ on the system state within $[x_0,x_0+l]$ at time t. The specification of $\lambda (x_0,l;t)$ and the way that it is used to sample eddy events are described in Section 6.3.

Although the profiles $u_i$ influence ODT advection by eddy events, they are not velocities operationally in the sense of appearing in a $\mathbf{u}·\nabla$ operator. Nevertheless, the statistical properties of these profiles are the principal model outputs, and it is important for these properties to faithfully represent the statistics of the physical flow field both for prediction purposes and in order for the SGS closure described here to be internally consistent. These requirements guide the formulation of ODT and its SGS implementation.

For given initial conditions, boundary conditions, and body forcings, if any, an ensemble of simulated realizations can be generated by running multiple realizations using a different random number seed for each. In this sense, ODT simulations are used to approximate the statistics of the ensemble, where the precision of the approximation depends on the number of realizations or, for statistically stationary flows, on the duration of the realization if a single realization is simulated. This raises the question of whether a formal mathematical representation of the underlying ensemble can be constructed. Such a construction is possible but has not previously been reported, so it is presented in Appendix A.

2.2. Coupling of ODT Domains

2.2.1. Lagrangian Formulation

ODT, as described in Section 2.1, differs from autonomous standalone ODT in that it includes the LES-prescribed pressure. This is one of the two types of large-scale input needed for SGS closure. The other required input is the set of LES face-normal velocities for all CVs, which are used to couple ODT domains in adjacent CVs.

This coupling is performed using a Lagrangian method that preserves the fine-scale flow features that are resolved within ODT. The approach is illustrated by comparing it to conventional Eulerian advective coupling of LES CVs involving a finite-volume treatment that does not resolve sub-structure within individual CVs.

For this purpose, consider constant-density finite-volume advection of a passive scalar field in 1D using a uniform mesh with cell width W. Continuity requires that all face velocities $\overline{v}$ are identical, where the coordinate index is omitted for this 1D illustration. The discretized scalar field has scalar value $s_j$ in cell j. For the purpose of streaming the scalar along the 1D domain with velocity $\overline{v}>0$, treat the scalar s as a function of the continuum coordinate x. Because s is uniform within each cell, it is piecewise constant in x. Streaming for a time interval $\Delta t<W/\overline{v}$ produces a shifted scalar profile $s^{*}\left(x\right)=s(x-\overline{v}\Delta t)$. Box filtering within each cell gives, after rearrangement, $(s_j^{*}-s_j)/\Delta t+\overline{v}(s_j-s_{j-1})/W=0$, which is the lowest-order discrete upwind approximation of the Eulerian 1D advection equation

$$s_{,t}+\overline{v}{s}_{,x}=0.$$

(3)

Although the spatial-continuum picture reduces to an Eulerian numerical algorithm, physically it prescribes Lagrangian time advancement as follows:

• Interpret the length-W portion of $s\left(x\right)$ within each cell j as a separate length-W 1D domain denoted $s(x;j)$.
• For each j, remove an interval $\overline{v}\Delta t$ from the downstream end of $s(x;j-1)$ and attach it to the upstream end of $s(x;j)$.
• Box filter each domain $s(x;j)$.

One consequence of the error inherent in the numerical approximation is that s is subject to numerical dissipation that is asymptotically diffusive over a long time interval with numerical diffusivity of order $W\overline{v}$. The specific cause of this dissipation is the final box-filtering step. This step is omitted when using the procedure to couple adjacent ODT domains because the intent is then to preserve rather than suppress fine-scale structure. ODT advancement processes dissipate fluctuations induced by this coupling in a manner that models the physical processes governing dissipation, so dissipation in ODT is not a numerical artifact. Nevertheless, the Lagrangian coupling introduces property discontinuities that can be partially or totally mitigated depending on the flow regime, as explained shortly.

2.2.2. Domain Coupling in the Linear-Eddy Model

Although not yet implemented numerically for ODT SGS closure, this domain-coupling procedure is commonly used for SGS closure involving the linear-eddy model (LEM), which is the antecedent of ODT. Velocity profiles are not defined or time advanced in LEM. LEM evolves scalar fields on a 1D domain using a parameterized prescription of the statistics of eddy-event sampling (where an LEM eddy event is solely a triplet map) rather than a self-contained procedure. It is used as an SGS mixing and combustion closure for LES. Features common to the LEM and ODT SGS formulations are explained here in terms of LEM, omitting details documented elsewhere [1,2] that are not relevant here. As in LEM SGS implementation, the domain-coupling procedure described in Section 2.2.1 is termed splicing.

Each exemplar object $s(x;j)$ receives inflow from the upstream side, which switches from one cell boundary to the other if the sign of $\overline{v}$ is reversed. In contrast, each LEM domain has permanently assigned inflow and outflow ends, where the inflow is always received from the whatever direction is upstream at a given instant and the outflow is always directed in the current downstream direction. In a 3D LES CV, the LEM domain is analogous to a weather vane, whose center is fixed at the center of the LES CV. All inflows through the CV faces are appended to the tail while the current contents of the LEM domain are shifted commensurately toward the head. Outflows are produced by removing fluid intervals at the head of the domain and transferring one interval through each of the faces that, based on the sign of its face-normal velocity, is an outflow face. This process is illustrated in Figure 2.

The LEM state associated with an LES CV is intended to be statistically representative of the physically correct population of SGS chemical compositions within the CV so that it provides accurate chemical closure information. For this purpose, LEM time advancement should capture the effects of all eddy motions not resolved at LES scales. This includes eddy motions as large as the size W of the CV, corresponding to eddy events of that size, so the LEM domain length D should be at least that large. This point is discussed further with reference to ODT in Section 7.1.

Another constraint is that the typical residence time of an LEM fluid parcel within an LEM domain should be roughly $W/\overline{v}$, which is the LES-prescribed flow-through time of the associated CV. This is desired because processes such as finite-rate chemistry within a CV imply a dependence of the chemical state of fluid exiting the CV on its residence time. More broadly, consistency between the SGS and LES advancement requires the SGS fluid parcels to follow trajectories through the flow domain that are consistent with the tracer statistics implied by the LES-resolved time advancement. Residence-time consistency is one facet of this requirement.

This ostensibly imposes the requirement $D=W$, but as shown in Section 7, an additional degree of freedom is available to enforce the LEM residence time $W/\overline{v}$ while allowing $D>W$ if the latter is desired, e.g., in order to reduce the statistical sample variability of outputs extracted from the LEM state.

2.2.3. Extension to ODT

The original ODT formulation [3] time advanced only one velocity component, denoted here as u. Although the SGS formulation requires the vector ODT formulation [4,5], some aspects of splicing are explained first with reference to a scalar velocity u for clarity. For this purpose, the configuration used to describe the application of splicing to the scalar field s is adopted, with u replacing s. Then on some ODT domain, whose x coordinate ranges, say, from 0 to W, splicing inserts some u profile into an interval $[0,x^{*}]$, while $[x^{*},W]$ now contains the u profile that previously occupied the interval $[0,W-x^{*}]$. Owing to time advancement on all ODT domains prior to splicing, the profile $u\left(x\right)$ after splicing is in general discontinuous at $x^{*}$.

This discontinuity is an artifact of splicing that needs to be removed because the u profile governs eddy-event sampling (Section 2.1), and in particular generates an unphysical burst of small eddy events in the vicinity of a discontinuity. Two ways of removing the discontinuity are compared.

One approach is to box filter discontinuous velocity profiles. This conserves momentum but dissipates kinetic energy, with no kinetic energy remaining to be dissipated after that, as indicated by the fact that the kinetic-energy dissipation rate is then zero. This is analogous to the lowest-order ’production equals dissipation’ level of turbulence modeling, or it alternatively can be seen as analogous to implicit LES (ILES), in which the numerical algorithm provides the closure [6].

The other approach is to use the K and J kernels (see Section 2.1) to remove the discontinuity in a globally conservative manner. First, the piecewise constant function $h=A$ in $[0,x^{*}]$, $h=B$ in $[x^{*},W]$ is added to $u\left(x\right)$, where $A-B$ is set equal to the u discontinuity $\delta$ at $x^{*}$ so as to remove the discontinuity, and the condition $A{x}^{*}+B(W-x^{*})=\delta x^{*}+BW=0$ enforces momentum conservation. This gives $B=-\delta x^{*}/W$ and $A=\delta (W-x^{*})/W$. This adjustment causes some kinetic-energy change $\Delta E$. It is counteracted to render the scheme conservative by adding $cK\left(x\right)$ to $u\left(x\right)$, where K is a kernel that is applied in this instance to $[0,W]$ and c is adjusted so that it causes the energy change $-\Delta E$. (As explained in Section 6.1, this involves solving a quadratic equation for c, where one of the two solution branches is identified as the physical solution. The discriminant can be negative, in which case energy conservation cannot be enforced, so some violations must be allowed.) The definition of $K\left(x\right)$ implies ${\int }_0^{W}K\left(x\right)dx=0$, so the kernel addition does not change the total momentum on the ODT domain. For variable-density flow, both the K and J kernels are needed, involving an additional coefficient b (see Section 6.1), to enforce conservation as well as remove the discontinuity.

A CV in a 3D Cartesian mesh is subject to one splicing operation through each of its six faces during each time-advancement cycle. Then on average there are three segments appended to the inflow end of the associated ODT domain. The required adjustments can be implemented sequentially using the method described above after each attachment, or more efficiently with a slightly more complicated algebraic system that implements one adjustment accounting for all of the newly appended segments.

To sharpen the comparison of approaches, consider the special case in which the u profile on each ODT domain is initially spatially uniform, where the spatially uniform value ${\overline{u}}_j$ is different on different ODT domains, indexed by j. This is analogous to the initial condition of the $s\left(x\right)$ example. Then splicing induces ODT-level kinetic-energy nonuniformities, implying production of subgrid-scale turbulent kinetic energy (TKE). As noted, the first approach dissipates the TKE immediately. However, the adjustment after ODT splicing is globally (within each ODT domain) conservative and therefore is operationally an advective-transport mechanism rather than a dissipative process. This is consistent with the foundational principle that the effects of viscous and any other molecular transport processes on the ODT-level system states are implemented only through time advancement of molecular-transport equations, in the present instance Equation (1).

The practical consequence is that TKE dissipation in the ODT domains is the outcome of eddy-induced cascading of velocity fluctuations to the viscous scales, leading to viscous dissipation of those fluctuations. The resulting time lag from TKE production to transport, strain-induced scale reduction, and finally dissipation of TKE results in local imbalances between production and dissipation. Thus, these imbalances arise in a manner that is more closely linked to the governing physics than is the treatment provided by, e.g., a parameterized SGS kinetic-energy equation.

Advected scalars such as species mass fractions are not amenable to adjustments that remove scalar discontinuities owing to the stricter conservation laws that must be enforced locally in order to assure chemical realizability. Therefore, scalar discontinuities, including density discontinuities, are inherent artifacts of subgrid LEM and ODT formulations that involve splicing. They can be mitigated somewhat by coarsening the CVs, which reduces the splicing frequency. Because LEM and ODT fully resolve SGS turbulent cascading within the modeling framework, the desired fidelity might be achievable with coarser meshes than are usually used.

Numerical implementation considerations such as mesh resolution are similarly consequential with regard to algorithmic accuracy and stability. In past applications of LEM SGS closure of LES, benchmark tests of the splicing algorithm, such as rigid-body translation of a cylinder at some inclination relative to the Cartesian coordinates, have been conducted albeit not reported in publications. Mesh refinement studies demonstrated satisfactory convergence to the continuum limit, but for the mesh and time-step settings of practical computations, there was considerable dispersion and distortion. As indicated by the 1D analysis in Section 2.2.1, splicing is effectively first order in space and time with no available higher-order extension, so there is no obvious way to mitigate these numerical errors. In practice, this requires limitation of the application of splicing-based closure to cases for which these numerical effects are dominated by turbulent fluctuations. As shown, e.g., in [2], this encompasses a variety of regimes of interest.

Application of splicing to ODT domains is subject to these considerations as well as additional factors. First, in AME the CV face velocities that specify the volume fluxes governing the splicing operation are based on ODT flow states, as described next in Section 2.3. This introduces stochastic variations of those fluxes in space and time, which could promote numerical instability. Second, in CVs adjacent to walls, one endpoint of the associated ODT domain is attached to the wall. Splicing into or out of that endpoint could degrade the ODT representation of the near-wall flow. Strategies that avoid this, and accordingly differ from the procedure illustrated in Figure 2, are described in Section 4.

Overall, splicing is a necessary but imperfect mechanism for preserving the physical distinction between nondissipative advection and dissipative viscous and molecular-mixing processes. However, there are limited but significant applications for which purely Lagrangian LEM or ODT closure that eschews splicing is suitable. One such application is SGS mixing within the unresolved thickness of a time-varying internal fluid interface such as a flame front, a case in which the LEM domain orientation is locally normal to the front rather than locally aligned with the flow [7].

2.3. Time-Advancement Cycle

As in Section 2.2, the notional time-advancement cycle for a simple ILES analog is outlined to serve as a reference point for description of the low-Mach-number time-advancement cycle for ODT-based SGS closure:

• Within each CV, labeled by index j, the spatially uniform velocity ${\overline{\mathbf{u}}}_j$ is subject to pressure forcing that adds $-(\Delta t/\rho ){(\nabla p)}_j$ to ${\overline{\mathbf{u}}}_j$, where the components of $\nabla p$ are the quantities $p_{,i}$ in Equation (1) and $\Delta t$ is the cycle time step.
• The next step is a pressure projection using the CV velocities ${\overline{\mathbf{u}}}_j$ as input. Based on a staggered scheme, this yields updated CV pressures and CV face-normal velocities $\overline{\mathbf{v}}$ that are in conformance with the continuity equation.
• The updated face-normal velocities are interpolated to update the quantities ${\overline{\mathbf{u}}}_j$.
• Advection of the velocities ${\overline{\mathbf{u}}}_j$ as prescribed by the CV face-normal velocities is implemented as a 3D generalization of the discretized advection of the scalar $s(x;j)$ in Section 2.2.

The numerical diffusivity associated with the advection operation, noted below Equation (3), is not necessarily an accurate representation of SGS advective transport, which is why the algorithm is termed an ILES analog rather than an ILES. In any case, this sequence of operations is LES-like advancement of the momentum equation although it is not a conventional LES procedure.

The analogous time-advancement cycle for ODT-based SGS closure is as follows:

• Based on Equation (1) or some generalization of it such as in Section 5, each ODT domain is time advanced for an elapsed time $\Delta t$. This involves sub-cycling over smaller time intervals owing to the smaller CFL time for the well-resolved ODT mesh. This advancement is subject to interruptions to implement eddy events. No-flux boundary conditions are applied at the ODT domain endpoints unless a wall boundary condition is applied at one of the endpoints. The ODT wall treatment is explained in Section 4.
• For each ODT domain, the average velocity vector $\overline{\mathbf{u}}$ is evaluated by box filtering the ODT flow state. The resulting filtered flow field is interpolated to evaluate the face-normal velocities $\overline{\mathbf{v}}$, which are the inputs to a pressure projection that enforces consistency of the face-normal velocities with the spatially discretized continuity equation, as required for consistent implementation of splicing. The momentum interpolation in conventional pressure projection schemes [8,9] needs modification owing to the nonstandard AME procedures for momentum advancement (eddy events and splicing).
• For each CV, the quantity $p_{,i}$ in Equation (1) now has a new value based on the updated CV pressures that is denoted $p_{,i}^{*}$. Accordingly, the fixed pressure-gradient value that should have been used during the time advancement of the associated ODT domain is the interpolant ${\overline{p}}_{,i}\equiv (p_{,i}+p_{,i}^{*})/2$. If the ODT domain is not wall-bounded, then $u_i\left(x\right)$ is modified for consistency with advancement during $\Delta t$ based on ${\overline{p}}_{,i}$ rather than $p_{,i}$ by adding $C_p=-(1/\rho )({\overline{p}}_{,i}-p_{,i})\Delta t$, uniformly over x, to the profile $u_i$. If the ODT domain is wall-bounded, this would violate the no-slip boundary condition at the wall, reflecting a splitting error resulting from not applying the pressure-gradient forcing in Equation (1) on a fully synchronous basis. For this case, an approach that was shown to work well in [10] is applied. Namely, instead of incrementing $u_i\left(x\right)$ uniformly over x, the correction $2(x/D)C_p$ is added to $u_i\left(x\right)$. This yields a consistent $\overline{u_i}$ value and obeys the no-slip condition at the wall location $x=0$.
• Splicing is applied as in Section 2.2.3.

A key feature of the approach is that steps 2–4 fulfill the requirements for coarse-grained advancement of the momentum equation, complementing the fine-grained advancement during step 1. Splicing requires coarse-grained information that is ultimately provided by the pressure projection, which is a correction that enforces a constraint rather than a time-stepping procedure. The correction involves the use of an auxiliary equation that is structurally a coarse-grained momentum equation but is not time advanced, as explained Appendix B.

On this basis, it is fair to say that this formulation involves no time advancement of a filtered momentum equation. In that sensem the formulation is an outgrowth of prior efforts to construct a bottom-up multiscale turbulence simulation based on ODT, notably the approach termed ODTLES [11,12,13]. The present approach is structured more like LES/LEM than like ODTLES, which involves direction-dependent hybridization of resolved and coarse-grained momentum advancement. Because the new feature is the avoidance of any time advancement of a filtered momentum equation, the approach is termed autonomous microstructure evolution (AME). This feature does not assure consistently advantageous cost/performance outcomes relative to ODTLES or conventional LES. Future studies will identify the classes of problems for which the respective approaches are advantageous.

As in ODTLES, AME distinguishes between velocity profiles $\mathbf{u}$ that represent the flow states on the ODT domains and velocities $\overline{\mathbf{v}}$ that govern coarse-grained advection. The latter are in this respect auxiliary variables. The coarse-grained flow state represented by the set of velocities $\overline{\mathbf{u}}$ is likewise auxiliary from the ODT perspective, but from the coarse-grained perspective, the $\overline{\mathbf{u}}$ velocities constitute the desired flow solution. In LES implemented on a staggered mesh, the distinction between $\overline{\mathbf{v}}$ and $\overline{\mathbf{u}}$ is purely algorithmic, while in AME there is a further distinction because $\overline{\mathbf{u}}$ is a coarse-grained approximation of the detailed physical representation of the evolving flow state produced by ODT. A nuance in this regard is that there is a fine-grained representation of the advection induced in AME by the velocities $\overline{\mathbf{v}}$, as explained in Section 3.

Finally, it is noteworthy that the velocities $\overline{\mathbf{v}}$ are adjusted to satisfy the 3D continuity equation but no such requirement is imposed on the ODT velocity profiles $u_i$. Unlike $\overline{\mathbf{v}}$, these are not advecting velocities although they influence the sampling of eddies that implement 1D advection on the ODT domain. Eddy advection is performed by triplet maps that obey the 1D analog of continuity by construction - namely, the total arc length of the mapped interval(s) is the same after a triplet map as before.

3. Reynolds Stress and Subgrid-Scale Kinetic Energy

The analogy between AME and LES that is described in Section 2.3 applies to physical interpretation as well as numerical modeling of time-advancement processes, particularly with regard to the SGS TKE budget. Features of SGS TKE production, transport, and dissipation in the ODT SGS closure are discussed conceptually in Section 2.2.3. They are now considered with reference to the formal definitions of these quantities in LES.

SGS TKE production is more precisely described as the conversion of coarse-grained kinetic energy into SGS TKE. One factor in the expression for the conversion rate is the Reynolds stress, whose evaluation in AME is considered first.

Splicing is the AME implementation of coarse-grained advection. Splicing is CV-to-CV transfer of portions of ODT velocity-component profiles, generically denoted u, where the transfer is governed by CV face velocities, generically $\overline{v}$. (Vector operations and the associated component indexing conform to standard conventions so scalar notation is used for clarity.) The u fluctuation profile on an ODT domain is denoted $u^{\prime }=u-\overline{u}$. The SGS TKE (per unit mass in what follows) is generically denoted $q_{SGS}=\frac{1}{2}\overline{u^{\prime }{}^2}$, ignoring both component indexing and summation over components. Consistent with this, a 1D array of CVs is considered (or equivalently, a multidimensional array with nonzero spatially uniform $\overline{v}$ in only one direction), as in Section 2.2.1.

On an ODT domain $[0,D]$, suppose that a sub-interval $[x^{*},x^{*}+\Delta ]$ is designated to be transferred to the ODT domain in the CV to its right, corresponding to positive $\overline{v}$. In a coarse-grained sense, the displacement of this sub-interval is $+W$, while the remainder of the u profile is unmoved by this transfer in a coarse-grained sense. (As explained in Section 7, D can exceed W while nevertheless representing the flow state in the width-W extent of the CV in direction x. Here, the default case $D=W$ is chosen for simplicity.) On this basis, the coarse-grained displacement of the u profile is zero except for the selected size-$\Delta$ interval. The advecting velocity governing the displacement of that interval is denoted $v^{*}$. Accordingly, the profile of the advecting velocity $v\left(x\right)$ in $[0,D]$ is

$$v=\left\{\begin{array}{cc}{v}^{*}\hfill & \mathrm{i}\mathrm{f}{x}^{*}\le x\le x^{*}+\Delta \hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}\right.$$

(4)

Box filtering of v over the ODT domain gives $\overline{v}=(\Delta /D)v^{*}$ and hence $v^{*}=(D/\Delta )\overline{v}$. This identifies the ODT profile of the velocity fluctuation $v^{\prime }=v-\overline{v}$ as

$$v^{\prime }=\left\{\begin{array}{cc}\left(\frac{D}{\Delta }-1\right)\overline{v}\hfill & \mathrm{i}\mathrm{f}{x}^{*}\le x\le x^{*}+\Delta \hfill \\ -\overline{v}\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}\right.$$

(5)

This shows that the face velocity can be interpreted as a box-filtered velocity, thus motivating the notation $\overline{v}$ that has been adopted.

Similarly, the fluctuation of the advected velocity is $u^{\prime }=u-\overline{u}$. Spatial averaging of the product of $u^{\prime }$ and $v^{\prime }$ then gives the contribution of the outbound transfer of the size-$\Delta$ interval to the box-filtered Reynolds stress,

$$\overline{u^{\prime }{v}^{\prime }}=({\overline{u}}^{\Delta }-\overline{u})\overline{v},$$

(6)

where ${\overline{u}}^{\Delta }$ denotes the average of u over $[x^{*},x^{*}+\Delta ]$ and $x^{*}+\Delta$ is the outflow endpoint D of the ODT domain.

The outbound fluid is replaced by inbound fluid occupying a size-$\Delta$ interval for which $x^{*}=0$ owing to its attachment to the inflow endpoint. The ODT state in the CV that supplies this fluid is denoted $u_{-}$. The filtered value of the $u_{-}$ interval that is transferred rightward, denoted ${\overline{u}}_{-}^{\Delta }$, is in general different from ${\overline{u}}^{\Delta }$, where the notation refers to the distinct size-$\Delta$ intervals containing the inbound and outbound fluid respectively. The inbound and outbound transfers convert the original profile u into a new profile U. Box filtering of U gives $\overline{U}=[(D-\Delta ){\overline{u}}^{D-\Delta }+\Delta {\overline{u}}_{-}^{\Delta }]/D$, where ${\overline{u}}^{D-\Delta }$ denotes the average of u over the size $D-\Delta$ region not removed during the outbound transfer.

The inbound fluid enters the CV through the face opposite to the outflow considered thus far. As noted, $\overline{v}$ has the same value at both faces, which is reflected by the equal sizes of the advected fluid intervals. The box-filtered Reynolds stress associated with the inbound transfer is

$$\overline{u_{-}^{\prime }{v}_{-}^{\prime }}=({\overline{u}}_{-}^{\Delta }-{\overline{u}}_{-})\overline{v},$$

(7)

where $v_{-}^{\prime }$ is specified by Equation (5) with u replaced by $u_{-}$.

The new total (filtered plus SGS) kinetic energy is $\frac{1}{2}\overline{U^2}$, which is the sum of the filtered and SGS kinetic energies, $\frac{1}{2}{\left(\overline{U}\right)}^2$ and $Q_{\mathrm{sgs}}^2=\frac{1}{2}\overline{U^{\prime }{}^2}$ respectively, where $U^{\prime }=U-\overline{U}$. Similarly, the initial SGS kinetic energy is $q_{\mathrm{sgs}}^2=\frac{1}{2}\overline{U^{\prime }{}^2}$, where $u^{\prime }=u-\overline{u}$. Then the splicing-induced SGS TKE change $\delta q_{\mathrm{sgs}}^2=Q_{\mathrm{sgs}}^2-q_{\mathrm{sgs}}^2$ is expressed as

$$\delta q_{\mathrm{sgs}}^2=\mathcal{C}+\mathcal{A}+\mathcal{T}+\mathcal{S}-\mathcal{V},$$

(8)

consisting of conversion $\mathcal{C}$ of filtered kinetic energy into SGS kinetic energy, advection $\mathcal{A}$ by the resolved flow, turbulent transport $\mathcal{T}$, subgrid transport $\mathcal{S}$, and viscous dissipation $\mathcal{V}$, respectively.

This specializes the conventional decomposition of $\delta q_{\mathrm{sgs}}^2$ [14] to the setup analyzed here. The conventional pressure and viscous transport terms are omitted because splicing involves no represention of SGS pressure fluctuations or viscous fluxes across cell faces. $\mathcal{V}$ is included to indicate that dissipation of SGS TKE occurs during ODT advancement, as explained in Section 2.2.3. Because the splicing operation is not dissipative, $\mathcal{V}$ is omitted from the application of Equation (8) to splicing that follows.

The approach is to evaluate $\mathcal{C}+\mathcal{S}=\delta q_{\mathrm{sgs}}^2-\mathcal{A}-\mathcal{T}$ using the appropriate specializations of the terms on the right-hand size and then to compare the outcome to the conventional expression for $\mathcal{C}$ and $\mathcal{S}$. Extending interval-restricted averaging to arbitrary powers p of u and $u_{-}$, the identities

$$\overline{u^p}=\left[\Delta {\overline{u^p}}^{\Delta }+(D-\Delta ){\overline{u^p}}^{D-\Delta }\right]/D$$

(9)

and

$$\overline{U^p}=\left[\Delta {\overline{u_{-}^p}}^{\Delta }+(D-\Delta ){\overline{u^p}}^{D-\Delta }\right]/D$$

(10)

will be useful in what follows. As before, each size-$\Delta$ averaging interval begins at the lower boundary $x^{*}$ of the removed or inserted interval and $D-\Delta$ refers to the complement of that interval on the ODT domain.

Evaluation of $\delta q_{\mathrm{sgs}}^2=\frac{1}{2}\left[\overline{U^2}-\overline{u^2}-{\left(\overline{U}\right)}^2+{\left(\overline{u}\right)}^2\right]$ using Equations (9) and (10) gives

$$\delta q_{\mathrm{sgs}}^2=\frac{\Delta }{2D}\left({\overline{u_{-}^2}}^{\Delta }-\Delta {\overline{u^2}}^{\Delta }\right)-\frac{{\Delta }^2}{2D^2}\left[{\left({\overline{u}}_{-}^{\Delta }\right)}^2-{\left({\overline{u}}^{\Delta }\right)}^2\right]-\left(\frac{\Delta }{D}-\frac{{\Delta }^2}{D^2}\right)({\overline{u}}_{-}^{\Delta }-{\overline{u}}^{\Delta }){\overline{u}}^{D-\Delta }.$$

(11)

In accordance with the relation $\Delta =\overline{v}\Delta t$, this change is attributed to the time-advancement step $\Delta t$. For simplicity, Equation (11) is specialized to vanishing $\Delta t$, so $\Delta /D\ll 1$. Noting that ${\overline{u}}^{D-\Delta }$ converges to $\overline{u}$ in this limit, the leading-order result

$$\delta q_{\mathrm{sgs}}^2=\frac{\Delta }{2D}\left[{\overline{u_{-}^2}}^{\Delta }-{\overline{u^2}}^{\Delta }-2({\overline{u}}_{-}^{\Delta }-{\overline{u}}^{\Delta })\overline{u}\right]$$

(12)

is obtained.

In the conventional decomposition of ${\partial }_t{q}_{\mathrm{sgs}}^2$, $\mathcal{A}=-{\partial }_x\left(q_{\mathrm{sgs}}^2\overline{v}\right)$ and $\mathcal{T}=-{\partial }_x(\overline{v{u}^2}-\overline{v}\overline{u^2})$ in the present reduced notation. In Equation (8), a difference of $q_{\mathrm{sgs}}^2$ values separated in time by $\Delta t$ is analyzed, so the conventional terms are multiplied by $\Delta t$, which will lead to cancellation of $\overline{v}$ based on $\Delta =\overline{v}\Delta t$. Additionally, each spatial derivative is replaced by $1/W$ times a difference between the CV in which the $q_{\mathrm{sgs}}^2$ is analyzed and the CV to its left, whose velocity profile has been denoted $u_{-}$. Finally D is changed to W owing to the specialization to $D=W$.

On this basis,

$$\mathcal{A}=\frac{\Delta }{2W}\left(\overline{{u^{\prime }}_{-}^2}-\overline{u^{\prime }{}^2}\right)$$

(13)

and

$$\mathcal{T}=\frac{\Delta }{2W}\left({\overline{u_{-}^2}}^{\Delta }-{\overline{u^2}}^{\Delta }-\overline{u_{-}^2}+\overline{u^2}\right)$$

(14)

are obtained, giving

$$\mathcal{T}+\mathcal{A}=\frac{\Delta }{2W}\left({\overline{u_{-}^2}}^{\Delta }-{\overline{u^2}}^{\Delta }-{\overline{u}}_{-}^2+{\overline{u}}^2\right).$$

(15)

In Equation (8), with $\mathcal{V}$ omitted, substitution of results for the terms that have been evaluated gives

$$\mathcal{C}+\mathcal{S}=\frac{\Delta }{2W}\left[{\overline{u}}_{-}^2-{\overline{u}}^2-2\left({\overline{u}}_{-}^{\Delta }-{\overline{u}}^{\Delta }\right)\overline{u}\right].$$

(16)

Rewriting this as

$$\mathcal{C}+\mathcal{S}=-\left[\overline{v}\frac{\overline{u}+{\overline{u}}_{-}}{2}\frac{\overline{u}-{\overline{u}}_{-}}{W}-\overline{v}\overline{u}\frac{{\overline{u}}^{\Delta }-{\overline{u}}_{-}^{\Delta }}{W}\right]\Delta t,$$

(17)

each term is the product of a discretized derivative of filtered u velocity and a flux-type factor involving only filtered fields. This result is compared to the conventional expressions $\mathcal{C}=-\overline{u^{\prime }{v}^{\prime }}\frac{\overline{u}-{\overline{u}}_{-}}{W}\Delta t$ and $\mathcal{S}=(1/W)(\overline{u^{\prime }{v}^{\prime }}\overline{u}-\overline{u_{-}^{\prime }{v}_{-}^{\prime }}{\overline{u}}_{-})\Delta t$, where the spatial discretization reflects the evaluation of the Reynolds stresses as domain averages rather than point values at domain endpoints. This gives $\mathcal{C}+\mathcal{S}=(1/W)(\overline{u^{\prime }{v}^{\prime }}-\overline{u_{-}^{\prime }{v}_{-}^{\prime }}){\overline{u}}_{-}$. Using Equations (6) and (7) to evaluate the Reynolds stresses,

$$\mathcal{C}+\mathcal{S}=\left[-\overline{v}{\overline{u}}_{-}\frac{\overline{u}-{\overline{u}}_{-}}{W}-\overline{v}{\overline{u}}_{-}\frac{{\overline{u}}^{\Delta }-{\overline{u}}_{-}^{\Delta }}{W}\right]\Delta t$$

(18)

is obtained.

Equation (17) reproduces the discretized derivative terms in Equation (18). To fully reproduce Equation (18), $\overline{u}$ must be replaced by ${\overline{u}}_{-}$ where $\overline{u}$ appears in each of the terms multiplying a discretized derivative. This is equivalent to adding ${\overline{u}}_{-}-\overline{u}$ to each of these occurrences of $\overline{u}$. In the limit of vanishing W, cell-to-cell variation becomes asymptotically small such that ${\overline{u}}_{-}-\overline{u}$ is asymptotically small relative to $\overline{u}$, hence negligible. This shows that the difference between Equations (17) and (18) an order-of-accuracy discretization effect rather than a difference in the physics that is represented, thereby demontrating the physical correspondence between the subgrid stress associated with splicing and the subgrid stress inferred from the conventional LES formalism.

The term $\frac{\Delta }{2W}({\overline{u}}_{-}^2-{\overline{u}}^2)$ in Equation (16) is of divergence form (first difference of a function of a filtered quantity), identifying it as the subgrid transport term $\mathcal{S}$. Hence, the remaider of Equation (16) corresponds to the conversion term $\mathcal{C}$. This completes the evaluation of the terms of the $q_{\mathrm{sgs}}^2$ budget.

The use of closure models that resolve the processes governing $q_{\mathrm{sgs}}^2$ evolution eliminates the need for modeled equations for $q_{\mathrm{sgs}}^2$ that are otherwise used for closure in various LES formulations [14]. Within the AME framework, the evaluation of $q_{\mathrm{sgs}}^2$ budget terms is for diagnostic purposes rather than for LES closure and serves to highlight the correspondence to the conventional LES formalism. Additionally, $q_{\mathrm{sgs}}^2$ budget diagnosis of AME simulation results potentially can aid in the formulation of modeled SGS equations for conventional LES, which are often elaborate, involving multiple equations [14]. This is an example of the potential use of AME to improve more economical LES formulations.

Moreover, the formal analysis of Reynolds stresses is directly involved in AME numerical implementation. In Appendix B, the Reynolds stresses given by Equations (6) and (7) are needed inputs to the evaluation of eddy viscosities in the equation solved by the pressure projection.

4. Extension to Wall-Bounded Flow

The formulation thus far addresses bulk-flow closure. Wall boundary layers require special treatment, in accordance with previous ODT near-wall treatments. Previous applications of ODT to wall boundary layers include multiphysics cases [15,16,17], and ODT has been implemented as a near-wall SGS closure in LES of channel flow [10,18]. That formulation is potentially suitable for near-wall closure within the present framework, but it involves under-resolved time advancement, so in keeping with AME concept, an alternative procedure is outlined.

A CV in a cuboidal flow domain has as many as three faces that are adjacent to walls. Such a CV is assigned as many ODT domains as the number of such faces, where each domain is associated with one of those faces and is deemed to be normal to that face. The domains associated with a given CV are treated as independent during ODT time advancement. As in [10], each domain is attached at one of its endpoints to its associated face and ODT advancement within the domain is subject to the boundary condition at that endpoint, while a no-flux boundary condition is applied at the other endpoint.

Consider splicing for the case in which each CV involved in this operation has one face adjacent to a wall. Each of the associated ODT domains has only one free endpoint available for splicing, so this becomes both the inflow and outflow point. All assigned outflows emanate sequentially from that endpoint, after which all assigned inflows are attached sequentially at that endpoint.

A limitation of this procedure is that there is no direct flow from the innermost sub-layer in a wall-bounded ODT domain to its downstream neighbor. However, there is an indirect communication path, as follows. Wall-normal turbulent transport in the ODT domain associated with the upstream CV can advect near-wall fluid close enough to the outflow endpoint of the ODT domain so that it can be spliced to the ODT domain in the downstream CV, followed by wall-normal turbulent transport down to the near-wall region of that ODT domain. Owing to the slowness of near-wall flow, this is the dominant mechanism of streamwise transport of near-wall fluid over large streamwise distances, but not necessarily between adjacent CVs. These assertions are supported by quasi-steady-quasi-homogeneous theory and its validations [19], which indicate that small-scale near-wall motions (resolved by ODT in AME) depend primarily on the larger-scale local flow state (controlled by splicing in AME).

In contrast, the method in [10,18] includes Eulerian implementation of direct advective transport, at each height above the wall, between parallel ODT domains in neighboring wall CVs. This is under-resolved time advancement of advection, hence outside the scope of the AME concept, but its adoption could be advantageous, depending on the application.

The treatment of CVs for which $f>1$ faces are adjacent to the wall is based on the interpretation of the f associated ODT domains as f instances of the representative state within the CV. Then the average over all f domains is deemed to be the average state of the CV, e.g., for the purpose of evaluating $\overline{\mathbf{u}}$. Likewise the $\overline{\mathbf{u}}$ adjustment based on the pressure projection is applied to all f domains, so the individual domains will have different filtered velocities.

For splicing in a block-cuboidal domain with cubic meshing, the possible pairs of f values in the donor and receiver CVs are $(0,0)$, $(0,1)$, $(1,0)$, $(1,1)$, $(1,2)$, $(2,1)$, $(2,2)$, $(2,3)$, $(3,2)$, and $(3,3)$. Tetrahedral edge and corner CVs in an unstructured mesh allow the additional possibilities $(0,2)$, $(2,0)$, $(0,3)$, $(3,0)$, $(1,3)$, and $(3,1)$ but exclude some of the cubic-mesh pairings. Further, for $f>1$, splicing involving a pair of ODT domains within a given CV is possible. Specification of physically consistent splicing protocols for the various cases is straightforward but tedious and will be addressed as such cases arise in future applications.

5. Extension to Variable Density

5.1. Boussinesq Approximation

The model description thus far has been restricted to constant-property flow in order to introduce the approach in the simplest possible context. The first step in the extension of the scope of the model is its generalization to variable density, which is considered next.

The minimal dynamical effect of density is obtained in the Boussinesq approximation, in which density is held fixed at a reference value except in a gravitational forcing term proportional to the small local deviation of the density from the reference value. In ODT, the representation of this forcing depends on the ODT domain orientation.

On a horizontal domain, the forcing is applied during viscous advancement by adding the gravitational forcing term to Equation (1) for the velocity component that is vertically oriented. Then relabeling the density $\rho$ in Equation (1) as ${\rho }_0$ to indicate that it is a spatially uniform reference density, the gravitational forcing $-g(\rho -{\rho }_0)/{\rho }_0$ is added to the equation for the vertical velocity component, nominally the $i=1$ equation. Owing to spatial variation of the density profile $\rho \left(x\right)$, time advancement of the modified equation tends to produce velocity fluctuations along the ODT domain, thereby contributing to eddy occurrences as governed by the sampling procedure described in Section 6.3. Variants of this approach have been applied to buoyancy-driven boundary layers along vertical walls [16,20,21].

On a vertical domain, a triplet map applied to a domain interval in which $\rho$ varies will in general change the eddy-integrated gravitational potential energy. The energy change is an input to the kernel procedure, which imposes an equal-and-opposite change of the eddy-integrated kinetic energy of the velocity profiles (see Section 6.1) and consequently influences the sampling of eddy occurrences (see Section 6.3). Details of this approach evolved in a series of publications [3,15,22,23,24,25,26,27].

Reflecting the inflow-outflow treatment of non-wall ODT domains, each of them is deemed to be aligned with the current orientation of its time-varying filtered velocity $\overline{\mathbf{u}}$. Accordingly, the gravity vector is decomposed into components perpendicular and parallel to $\overline{\mathbf{u}}$. The effect of the perpendicular (parallel) component is then represented as it would be on a horizontal (vertical) ODT domain. In this way, the formulation reduces to the implementations for the respective special cases.

Because $\overline{\mathbf{u}}$ and therefore its orientation vary during ODT advancement, the domain orientation should ideally be re-evaluated continually for accurate treatment of the gravitational forcing. Since the variation of orientation is caused primarily by filtered-scale processes (splicing and pressure correction), it is nevertheless sufficient to re-evaluate it once per time-advancement cycle in keeping with splitting errors inherent in the multiscale formulation. (From this viewpoint, the pressure forcing in Equation (1) can likewise be applied only once per time-advancement cycle with negligible overall loss of fidelity.)

For conservative tracking of gravitational potential energy changes, the ODT domain must be fully specified with reference to the fixed coordinate system. A natural convention is to locate the domain midpoint at the CV center (suitably defined for CVs in an unstructured mesh). Then an orientation change corresponds to a rotation with respect to the domain midpoint. Unless the axis of the rotation is vertical, this implies vertical displacements of fluid elements. Provided that this motion is implemented once per time-advancement cycle and specifically in conjunction with splicing, the associated change of gravitational potential energy can be compensated by a change of kinetic energy that contributes to the overall kinetic-energy change during the kernel operation associated with the removal of splicing-induced discontinuties of velocity profiles, as described in Section 2.2.3.

This mechanism can likewise subsume gravitational-potential-energy effects associated with all fluid displacements that occur during splicing, including fluid shifts along the ODT domain resulting from outbound and inbound fluid transfers at opposite ends of the domain. There is splitting and spatial-discretization error associated with the displacements that occur during this operation, but it is again within the scope of error that is inherently associated with the coarse-grained evolution.

A principal application of Boussinesq buoyant stratified flow is geophysical flow. For this application, ODT SGS resolution can extend the scale range of the simulation relative to conventional LES, but not to the finest scales of turbulent motion. In this context, ODT needs closure at its resolution scale. Standalone ODT with such a closure has been used for simulation of the atmospheric boundary layer [15,28]. ODT and LEM have also been used as SGS scale-range extensions for the purpose of moist-thermodynamics closure in LES of clouds [29,30], albeit with a different strategies for coupling ODT and LEM to LES than proposed here.

5.2. Density Variation without Dilatation

The full variable-density treatment in AME, with buoyancy effects omitted, is described next. If density is a spatially nonuniform advected quantity on the ODT domain with no density changes in the Lagrangian frame, then the density profile evolves within ODT solely by triplet mapping. (An example of this regime to which ODT has been applied is immiscible fluids advected by turbulence [31].) Density nonuniformity is nevertheless dynamically consequential in ODT because it influences eddy-event occurrences and the coefficients of the kernels that are applied during eddy events, as explained in Section 6. Further, Equation (1) is generalized to the variable-density form

$${\left(\rho u_i\right)}_{,t}={\left(\mu u_{i,x}\right)}_{,x}-p_{,i},$$

(19)

where $\mu$ is the dynamic viscosity.

These effects ultimately modify the CV face velocities and hence splicing as well as the outcome of pressure projection. The influence is two-fold. First, the evolution of the ODT velocity profiles $u_i$ is modified as indicated. Second, box filtering on the ODT domain now requires Favre (mass-weighted) averaging to obtain the mass-averaged velocities needed at the coarse-grained level.

This raises the point that coarse-grained mass-averaged velocities are based on box filtering over ODT domains, but domain intervals that are spliced across CV faces generally have different average densities than the ODT domains from which they are taken. Therefore, ODT cannot precisely match the mass fluxes implied by the box-filtered density and velocity fields. This is a feature rather than a deficiency, reflecting ODT-resolved property fluctuations, whose effects are not captured at the coarse-grained level. In any case, fluid transfers across CV faces are most consistently specified on a volumetric basis (see Section 7.1), so the mass-flux mismatch is immaterial to the implementation of splicing.

The effect of the pressure forcing in Equation (19) on momentum is spatially uniform, but with density variation this implies a source of velocity fluctuations that can contribute to the available energy that induces eddy events, as explained in Section 6. This is the ODT representation of baroclinic vorticity generation. Standalone ODT can at best (and only for compressible flow [32]) incorporate pressure gradients aligned with the domain, so baroclinic vorticity can be represented only in SGS ODT coupled to advancement of the coarse-grained pressure field.

The projection of the pressure forcing in the ODT domain direction implies a bulk acceleration of the domain that can alternatively be diagnosed from the time evolution of $\overline{\mathbf{v}}$. This acceleration is analogous to gravitational forcing. In variable-density flow it induces the Rayleigh–Taylor instability, to which ODT has been applied [27]. Thus, standalone ODT captures the instability based on prescribed gravitation or acceleration and SGS ODT can represent it as a response to the coarse-grained acceleration field.

A related mechanism involving the effect of body forces in variable-density flow is the Darrieus–Landau instability. It is specifically associated with dilatational effects, as discussed in Section 5.3 and Section 6.3, and is represented in standalone as well as SGS ODT implementation.

Unlike the ODT velocity profiles, density profiles cannot be adjusted after splicing to remove discontinuities, so these discontinuities persist and potentially adversely impact model fidelity with respect to both the eddy-event time history, as explained in Section 6, and the advancement of Equation (19). It will be important to assess the effects of this artifact on model performance. Some considerations in this regard, including possible mitigation strategies, are discussed in Section 7.1 and Section 8.

5.3. Dilatational Flow

If additionally there are Lagrangian density changes, then there is a commensurate dilatational flow along the 1D domain because it is treated as a closed system during ODT advancement. Then Equation (19) is generalized to

$${\left(\rho u_i\right)}_{,t}+{\left(\rho u_d{u}_i\right)}_{,x}={\left(\mu u_{i,x}\right)}_{,x}-p_{,i}+[{\left(\rho u_d\right)}_{,t}+{\left(\rho u_d^2\right)}_{,x}]\widehat{\mathbf{i}}·\widehat{\overline{\mathbf{u}}}.$$

(20)

Here, a dilatational velocity profile $u_d\left(x\right)$ is introduced that is based on 1D continuity in Lagrangian form, giving $u_{d,x}=-(1/\rho )(D\rho /Dt)$, where the material derivative $D\rho /Dt={\rho }_{,t}+u_d{\rho }_{,x}$ is determined by the physical processes driving the dilatation, e.g., thermal expansion. This determines $u_d$ to within an arbitrary constant, implying the freedom to advect the 1D domain by any uniform velocity. In Equation (20), $u_d$ advects the ODT momentum profiles $\rho u_i$, which are advected scalars from an operational viewpoint at the ODT level. Additionally, $u_d$ appears in a source term at the end of Equation (20) that is explained in what follows.

Dilatation changes the length of the ODT domain. As noted in Section 2.2, the model is consistent for any domain length equal to or larger than the CV scale. The handling of variable-length ODT domains during the time-advancement cycle is discussed in Section 7.1.

Recognizing that $u_d$ varies with time as well as with ODT domain location x, the dynamical implications at both the ODT-resolved and coarse-grained levels are considered. Bulk translation of the entire ODT domain is governed by CV face velocities that prescribe the splicing of fluid between ODT domains. The variable-density ODT-level processes that influence this coarse-grained transport have been noted in the discussion of dilatation-free variable-density ODT in Section 5.2. The net outcome is that global changes of the total momentum of the ODT domain are captured in Equation (20) by the pressure-gradient forcing, as prescribed by the pressure projection. Therefore, the undetermined constant in the evaluation of $u_d$ is assigned by requiring the implied domain-integrated momentum change ${\int }_0^D\rho u_{d}dx$ to be zero. This avoids any double-counting of effects captured by coarse-grained processes during the time-advancement cycle.

In the chosen frame, the spatial variation of dilatation can drive or suppress momentum fluctuations and is a source of kinetic energy that is generated by the mechanism driving the dilatation, such as exothermic chemical reactions. These consequences of dilatation are incorporated into ODT by the source term $[{\left(\rho u_d\right)}_{,t}+{\left(\rho u_d^2\right)}_{,x}]\widehat{\mathbf{i}}·\widehat{\overline{\mathbf{u}}}$ in Equation (20). This expression projects $[{\left(\rho u_d\right)}_{,t}+{\left(\rho u_d^2\right)}_{,x}]$, which is vectorially aligned with the ODT domain-orientation unit vector $\widehat{\overline{\mathbf{u}}}$, onto the unit vector $\widehat{\mathbf{i}}$ that is aligned with the coordinate direction i (see Figure 1).

As a special case, assume that the viscous and pressure terms in Equation (20) are zero, that the ODT domain is aligned with the coordinate $i=1$, and that $u_1=u_d$ and $u_2=u_3=0$ initially. For any subsequent time history of the $u_d$ profile, Equation (20) then implies that $u_1=u_d$, as required.

Because the components $u_i$ do not advect fluid along the ODT domain, there is a decoupling of the momentum and kinetic-energy changes such that the latter will not precisely reflect the conversion from thermal to kinetic energy that is induced by the dilatation. The overall self-consistency maintained by the time-advancement cycle suggests that the energy conversion process will be roughly conservative, but this remains to be demonstrated. An analogous treatment of dilatation in an application of standalone ODT to turbulent combustion [33] yielded accurate results.

Local accelerations associated with transient fluctuations of $u_d$, in conjunction with the density variations that drive them, induce the Darrieus–Landau instability. This instability mechanism is incorporated into ODT as a contribution to the frequency of eddy occurrence, as explained in Section 6.3.

6. The ODT Eddy Event

6.1. Eddy Implementation

Some details of ODT eddy implementation are explained. Except where new or modified features are introduced, physical justification and interpretation are minimal because they are provided in the cited references.

An ODT eddy event involves two operations. The first is a mapping denoted $x^{\prime }\left(x\right)$ and termed the triplet map that is applied to a selected (see Section 6.3) interval $[x_0,x_0+l]$. Within this interval, every property profile s is subject to three-fold compression and the mapped interval is filled with three copies of the compressed profile, but the middle copy is spatially flipped. Outside this interval, $s\left(x\right)$ is unchanged.

Each $x\in [x_0,x_0+l]$ is mapped to three locations $x^{\prime }$. It is sometimes convenient to refer to the single-valued inverse map $x\left(x^{\prime }\right)$ for given $x_0$ and l. Denoting the mapped property profile as $s^{\prime }$, this gives $s^{\prime }\left(x^{\prime }\right)=s\left(x\left(x^{\prime }\right)\right)$. Based on this meaning of $s^{\prime }\left(x^{\prime }\right)$, the generic argument x, and hence the notation $s^{\prime }\left(x\right)$, is used in what follows.

On a fixed uniform 1D mesh, the triplet map is approximated using a permutation of mesh cells [1]. An alternative approach that has proven to be advantageous is to use an adaptive mesh that allows implementation of the continuum definition of the triplet map as well as the continuum mathematical specfication of splicing procedures. Adaptive-mesh ODT implementation, as described in [21], is assumed here in all explanations of ODT and AME.

The second operation, which is applied only to velocity components, is the kernel addition shown in Equation (2). This operation is the model representation of scale-l pressure-fluctuation effects, reflecting the pressure-gradient term in the momentum equation. In this regard, the pressure-gradient term in Equation (1) is the filtered pressure computed during the pressure-projection step. No explicit model of ODT-resolved pressure fluctuations per se has been formulated for low-Mach-number flow.

The K and J kernels induce pointwise velocity changes and domain-integrated kinetic-energy changes. Additionally, J (always) and K (in variable-density flows) individually induce domain-integrated momentum changes, sometimes with equal and opposite effects so as to obtain zero net change. For a given velocity component, these changes can involve kinetic-energy exchange with other components as well as with momentum and kinetic-energy sources and sinks associated with coupling to processes involving, e.g., gravitational potential energy, surface-tension energy, or entrained particles or droplets. In Section 2.2.3, kernels are used to restore energy conservation after piecewise-uniform adjustments of velocity profiles to remove discontinuities.

These kernel effects and the applicable conservation laws constrain the coefficients $b_i$ and $c_i$ in Equation (2) as follows. For variable-density flow, momentum conservation implies

$$\int \rho \left(x\right)[u_i^{\prime }\left(x\right)+c_{i}K\left(x\right)+b_{i}J\left(x\right)]dV=\int \rho \left(x\right)u_i\left(x\right)dV+S_i.$$

(21)

Integration over $dV$ reflects the volume interpretation of the ODT domain based on the cross-sectional area $\gamma$ introduced in Section 7.1. The integrals extend over the entire domain, but the integrands are zero outside the eddy interval so they can be equivalently restricted to that x range. The argument x is henceforth omitted. Here and below, a previous constant-density treatment [34] that introduced momentum and energy sources (or sinks) $S_i$ and $S_E$ respectively is extended to variable-density flow by bringing $\rho$ into all integrands. This results in minor generalizations of the algebraic solutions in [34]. Details that are readily inferred from that treatment are omitted here.

Because the triplet map conserves total momentum, $u_i^{\prime }$ can be substituted for $u_i$ in the integrand on the right-hand side. This cancels a term on the left-hand side, reducing Equation (21) to

$$b_i=\frac{S_i-c_i\int \rho KdV}{\int \rho JdV}.$$

(22)

This relation can be substituted for $b_i$ in equations that follow, reducing them to relations for the unknowns $c_i$. The solution for $c_i$ then determines $b_i$ using Equation (22). Because $\int KdV=0$, Equation (22) directly evaluates $b_i$ if $\rho$ is spatially uniform.

The eddy-induced change of the domain-integrated kinetic energy of velocity component i is

$$\Delta E_i=\frac{1}{2}\int \rho {(u_i^{\prime }+c_{i}K+b_{i}J)}^{2}dV-\frac{1}{2}\int \rho u_i^{2}dV.$$

(23)

Conservation of total kinetic energy by the triplet map allows $u_i^{\prime }$ to be substituted for $u_i$ in the second integral, resulting in a cancellation of terms that gives

$$\Delta E_i=\frac{1}{2}\int \rho (2c_i{u}_i^{\prime }K+2b_i{u}_i^{\prime }J+2c_i{b}_{i}KJ+c_i^2{K}^2+b_i^2{J}^2)dV.$$

(24)

Substitution of Equation (22) for $b_i$ yields a quadratic equation for $c_i$, whose solution straightforwardly generalizes Equation (11) of [34].

$c_i$ is thus specified in terms of $\Delta E_i$ and functionals of the flow state. Energy conservation imposes the constraint ${\sum }_{i=1}^3\Delta E_i=S_E$, leaving two additional degrees of freedom governing the allotment of this net kinetic-energy change to the three components.

For this purpose, additional modeling is introduced. Although a constant-density treatment for any combination of momentum and energy sources [34] and a variable-density treatment involving only energy sources [5] are available, there is no prior variable-density treatment involving both types of sources. This is needed for the use of kernels in Section 2.2.3 to correct conservation violations caused by flow-field adjustments to remove splicing-induced discontinuities. Equation (24), in conjunction with the energy-allocation scheme introduced in Section 6.2 to specify the required energy changes $\Delta E_i$, provides the needed model extension.

Given the quantities $S_i$ and $\Delta E_i$ and the flow state, Equations (22) and (24) can be solved for $c_i$. It is possible for the discriminant of the solution to be zero or negative, corresponding to an energetically prohibited event. For this and other reasons (see Section 6.3), a sampled eddy event is not necessarily implemented.

If the discriminant is positive, then the quadratic equation has two real solutions. The solution branch is chosen on the basis of the functional dependence of $c_i$ on $\Delta E_i$. As $\Delta E_i$ monotonically approaches zero starting from its assigned value, one of the solution branches approaches $c_i=0$, corresponding to no change of the flow state in the absence of a required energy change. The solution based on the assigned $\Delta E_i$ value that corresponds to this branch is deemed to be the physically consistent solution.

6.2. Allocation of Kinetic-Energy Changes

Rotational eddy motion can redistribute kinetic energy among velocity components. Existing formulations of the kernel operation, such as in [34], take this into account. Depending on the application, the net kinetic-energy change $S_E$ might be specified on a component-wise basis, hence $S_{E,i}$ for each component i, but the rotational eddy motion similarly justifies modification of this allocation. (To preserve momentum conservation in ODT, which is a vector constraint on each individual momentum component, momentum is not redistributed among components.)

To narrow the range of options, it is noted first that there is a state-dependent bound on the amount of kinetic energy that can be extracted from a velocity component using kernels. It is determined by minimizing $\Delta E_i$ with respect to $c_i$ using Equation (24). The absolute value of this negative quantity is termed the component-i available energy $Q_i$. The variable-density expression for $Q_i$ is a straightforward generalization of Equation (12) of [34]. Through the dependence of Equation (24) on $b_i$, $Q_i$ depends on the momentum sources $S_i$.

The total available energy of the post-map velocity field is then $Q\equiv {\sum }_{i=1}^3{Q}_i$ and the net available energy, accounting also for energy sources, is $H=Q+S_E$. If $H\le 0$, then $S_E$ is an energy sink requiring a greater reduction in kinetic energy than the maximum possible reduction Q, so the eddy is deemed to be forbidden. (Criteria for eddy occurrences are explained in Section 6.3.) Therefore, nothing further needs to be said about this case, so what follows is restricted to $H>0$.

The simplest assumption is that the final eddy state has no memory of details of the initial state or external sources, so H is the sole governing quantity on which to base the enforcement of energy conservation. Then the final component-i available energies are $Q_i^{*}=\frac{H}{3}$ because there is no basis for different outcomes for different components.

This is a good baseline model and perhaps the only needed model. It implies the component-i available-energy change $\Delta E_i=Q_i^{*}-Q_i=\frac{H}{3}-Q_i=\frac{S_E}{3}+\frac{1}{3}(Q_j+Q_k-2Q_i)$, where subscripts j and k refer to the coordinate indices not equal to i. If $S_E=0$, this is the usual energy redistribution for $\alpha =\frac{2}{3}$, where $\alpha$ denotes the coefficient of $-Q_i$.

To extend this minimal formulation, a physical basis is needed for arriving at unequal values of the final component available energies $Q_i^{*}$. The initial component available energies $Q_i$ are the physically compelling basis for this because component energies are properties of the ODT physical state, involving no additional external input. This consideration, in conjunction with (1) the physical principle that an eddy region cannot supply more energy to an external sink than the eddy region contains, and (2) the additional straightforward, though not physically required, assumption that the dependence of $Q_i^{*}$ on $Q_i$ is linear, yields the generalized formulation

$$Q_i^{*}=\left[1+\chi \left(\frac{3Q_i}{Q}-1\right)\right]\frac{H}{3},$$

(25)

where $\chi$ is a ‘memory’ parameter. This recovers $Q_i^{*}=\frac{H}{3}$ for the ‘no-memory’ value $\chi =0$. Because $Q_i$ and H are non-negative and $Q_i\le Q$, $Q_i^{*}$ is also non-negative provided that $\chi \in \left[-\frac{1}{2},1\right]$. The definition of available energy precludes negative $Q_i^{*}$, so a value of $\chi$ outside this range could yield a mathematical inconsistency.

The definition of $\alpha$ as originally stated [4] does not uniquely specify the treatment of energy sources. For $S_E=0$, $\alpha$ and $\chi$ are linearly related according to $\chi =1-\frac{3}{2}\alpha$, where the $\chi$ range $\left[-\frac{1}{2},1\right]$ corresponds to the $\alpha$ range $[0,1]$. Then the equivalent result $Q_i^{*}=\left[1+\left(1-\frac{3}{2}\alpha \right)\left(\frac{3Q_i}{Q}-1\right)\right]\frac{H}{3}$ is obtained, which has the advantage that it uses the existing notational convention, but at the cost of less clarity with regard to the assumptions that extend the treatment to nonzero sources. However, the $\alpha$ notation does clarify one point. $\alpha$ is the fraction of the component-i available energy $Q_i$ that is transferred, in equal shares, to the other two components by the kernels, which requires its range to be $[0,1]$. This explains the less intuitive range $\left[-\frac{1}{2},1\right]$ of the parameter $\chi$.

Regardless of the number and types of physical processes contributing to $S_E$, this input suffices in this formulation to determine the quantities $Q_i$ and H and thus the energy changes $\Delta E_i=Q_i^{*}-Q_i$ that are used to calculate the coefficients $c_i$. On this basis, a wide range of physical processes and their couplings interact with eddy implementation (and sampling - see Section 6.3) within a unitary framework.

For spatially uniform $u_i^{\prime }$ and no energy or momentum sources, substitution of Equation (22) into Equation (24) indicates that the terms involving $u_i^{\prime }$ cancel. Then the integrand reduces to $\rho c_i^2{(K-J\int \rho KdV/\int \rho JdV)}^2$, which is non-negative. As expected for a state with zero kinetic energy, kernels can increase but not decrease the total kinetic energy, indicating that the component-i available energy is zero. Eddy selection, as described in Section 6.3, is based on available kinetic energy. In particular, positive available kinetic energy is required in a selected interval in order for an eddy to be implemented in that interval. Therefore, if all velocity components are uniform within the eddy range and there are no external sources, the eddy event is not implemented. Accordingly, the available-energy formalism satisfies the requirement that density variation in the absence of velocity fluctuations or an external source cannot induce eddy events. In buoyant stratified flow, the gravitational-potential-energy change resulting from triplet mapping a gravitationally unstable density profile is an external source that can generate eddy events in initially motionless fluid. This is the basis of ODT applications to buoyancy-driven flow that are discussed in Section 5.1.

6.3. Eddy Sampling

The physical reasoning and formal analysis on which eddy sampling is based are discussed in detail in references cited below. Here, the concepts are briefly outlined and the mathematical expression for the eddy rate distribution $\lambda (x_0,l;t)$ that is derived in those references is presented and explained.

As in mixing-length theory, $\lambda$ is expressed dimensionally in terms of the relevant length scale, in this case l, and a modeled time scale $\tau$, giving $\lambda ={\tau }^{-1}{l}^{-2}$. Unlike conventional mixing-length theory, the value of l is specific to each event rather than a fixed function of spatial location, and $\tau$ depends on the instantaneous system state in the size-l eddy interval rather than mean values of bulk metrics of the flow.

The model adopted for $\tau$ [35], re-expressed in present notation, is

$$\frac{1}{\tau }=C\sqrt{\frac{2}{{\widehat{\rho }}_{KK}{V}_{ϵ}{l}^2}\left(\frac{KK}{V_{ϵ}{l}^2}{Q}^{*}-Z{E}_{\mathrm{vp}}\right)},$$

(26)

where

• $KK\equiv \int K^{2}dV$, ${\widehat{\rho }}_{KK}\equiv \frac{1}{KK}\int \rho K^{2}dV$, and $V_{ϵ}$ is the eddy volume.
• $Q^{*}\equiv {\sum }_{i=1}^3{Q}_i^{*}$ is the net final available energy, now interpreted as the kinetic energy associated with the eddy motion.
• $E_{\mathrm{vp}}=\frac{V_{ϵ}}{2l^2}\frac{{\overline{\mu }}_{ϵ}^2}{{\overline{\rho }}_{ϵ}}$, where ${\overline{\mu }}_{ϵ}$ and ${\overline{\rho }}_{ϵ}$ are the $V_{ϵ}$-averaged dynamic viscosity and density, respectively. (${\overline{\mu }}_{ϵ}$ is a harmonic average.) $E_{\mathrm{vp}}$ is a viscous penalty that suppresses eddies that are subject to sufficient viscous damping to prevent their occurrence. (A negative argument of the square root indicates that the eddy has insufficient net final available energy to overcome viscous suppression.)
• C and Z are adjustable model parameters.
• The term $Q^{*}/\left({\widehat{\rho }}_{KK}{V}_{ϵ}{l}^2\right)$ is the underlying dimensional estimate of ${\tau }^{-2}$. Other quantities in Equation (26) are based on a variable-density formulation [5] that consistently generalizes an earlier constant-density expression for $\lambda (x_0,l;t)$, including equivalent definitions of model parameters.

With the choice $\chi =0$ in Equation (25), corresponding to the $\alpha$ value used almost all ODT studies to date, eddy implementation is effectively parameter-free and thus universal. Thus, all case-dependent behavior is embodied in the rate distribution $\lambda$ from which eddies are sampled. This is reflected in order-one variation of the best-fit parameter values for different classes of flows.

For unconfined flows that impose no intrinsic limitation on the largest eddy size $L_{\max }$, various empirical principles have been invoked to suppress unphysically large eddies [5]. In SGS implementation of ODT, C and $L_{\max }$ have been set by requiring the mesh-resolved and ODT SGS energy spectra to be complementary in wave-number range, without a gap or overlap, and to satisfy an amplitude matching condition at wave number $1/L_{\max }$ [11]. These are internal consistency requirements, involving no empirical input, so the only empirically tuned parameter is Z, which is tuned to match the empirical dissipative roll-off of the energy spectrum. (For wall-bounded ODT domains, Z is chosen to match the transition from the viscous layer to the buffer layer [10].) In AME, $L_{\max }$ will be of order W, so self-consistent determination of C and $L_{\max }$ can be performed using domain sizes $D>W$ that allow the possibility of $L_{\max }>W$. The implications of $D>W$ with regard to splicing are discussed in Section 7.1.

$\lambda (x_0,l;t)$ constrained by the assigned size bound $L_{\max }$ fully determines the statistics of eddy occurrences. Efficient sampling from this time-evolving distribution is needed for cost-effective numerical implementation. Because $\lambda$ depends only on the current system state with no explicit dependence on past history, eddy occurrences constitute a Poisson process in time. The thinning method, in which a Poisson event process is oversampled but events are selectively rejected so as to recover the correct sampling statistics, can therefore be used for eddy sampling. Details of the application of this method to sampling from $\lambda (x_0,l;t)$ are provided in the Appendix of [21].

The central role of the net available energy $Q^{*}$ in both the selection and implementation of eddy events has been highlighted. The scope of the phenomenology that this approach can capture is illustrated by ODT modeling of the Darrieus–Landau instability in premixed flames. Local accelerations induced by thermal expansion are equivalent to body forces that can result in the equivalent of triplet-map-induced changes of gravitational potential energy in variable-density flow, as explained in Section 5.3. This equivalent gravitational effect, represented in ODT as an external source contributing to the available energy so as to increase the amplitudes of the applied kernels, thereby promoting eddy occurrences, has been incorporated into combustion simulations as a representation of the Darrieus–Landau instability [16,33]. SGS ODT in AME can similarly capture this phenomenology.

7. Details of the Splicing Procedure

7.1. Enforcement of Consistent Volume Transfers

In this section, it is convenient to refer to ODT SGS closure of LES, but references to LES likewise apply to the coarse-grained fields in AME. The unclosed LES property fields that SGS ODT is intended to close are intensive properties. For this purpose, each ODT domain serves as a representative sample of flow states that can be averaged in order to obtain closure in the associated CV. This allows freedom to set the ODT domain length so that ODT provides a sufficient statistical sample of flow states for closure purposes, subject to other requirements. The relevant considerations are examined with reference to a cubic LES mesh and then generalized to an arbitrary unstructured mesh.

During an advancement step corresponding to a time increment $\Delta t$, the volume transfer across an LES CV face is $\overline{v}\Gamma \Delta t$, where $\overline{v}$ is the face-normal velocity and $\Gamma$ is the face area. In order to accommodate such LES-level transfers on a consistent volumetric basis, each ODT domain is assigned a nominal volume $\gamma D$ that is equal to the volume V of the associated CV, where D is the domain length and $\gamma$ is a nominal cross-sectional area, hence $\gamma =V/D$. To accommodate the LES-prescribed rate $\overline{v}\Gamma$ of volume inflow on the ODT domain, the nominal splicing-induced flow velocity along the ODT domain must be $\widehat{v}=\overline{v}\Gamma /\gamma$. The corresponding residence time of fluid on the ODT domain is $D/\widehat{v}=V/\left(\overline{v}\Gamma \right)$.

For the example in Section 2.2.1, the fluid residence time in a CV is $W/\overline{v}$. The cubic-mesh relations $W=V^{1/3}$ and $\Gamma =V^{2/3}$ reduce the ODT residence-time result $V/\left(\overline{v}\Gamma \right)$ to the same expression, showing that the proposed volumetric interpretation of the ODT domain is consistent in this regard.

Expressing this another way, the LES face velocity $\overline{v}$ implies a displacement of a fluid element of length $\overline{v}\Delta t$ across the face in time $\Delta t$, but for $D\ne W$, the corresponding length of the spliced ODT interval is $(D/W)\overline{v}\Delta t$. Although the volumetric picture involving $\Gamma$ and $\gamma$ is useful for motivating this splicing-length rule, numerical implementation requires only the application of this rule. In this regard, $\widehat{v}$ serves no purpose other than helping to motivate this rule and therefore has no physical significance nor does it have any role in model implementation.

As noted, $D>W$ may be desired based on sample-size considerations. Sample-size requirements related to the ODT representation of backscatter are discussed in Section 7.3. An ODT SGS closure of LES of clouds provided indications that $D=W$ produced insufficient sampling [30]. That formulation did not involve splicing, but $D=W$ was found to be sufficient in a cloud LES study using LEM SGS closure with splicing [29]. Both formulations involved dynamically active feedback of SGS moist thermodynamics to LES-resolved scales. ODTLES applications without such feedback did not indicate sample-size insufficiency, but the numerically dissipative domain coupling in ODTLES could have contributed to the suppression of fluctuations. These disparate observations do not point to any firm conclusions other than to say that the sample-size requirement might depend on both the application case and details of the model formulation.

Another consideration with regard to the domain size is that ODT eddy events should span a size range such that the combination of splicing and ODT eddy events captures the full spectral range of turbulent eddy motions, with neither double-counting nor gaps, as explained in Section 6.3. For a different form of ODT SGS closure, it was found that eddy events extending to size $L_{\max }=4W$ were required for this purpose [11]. This is qualitatively reasonable because the typical displacement of a fluid element by a triplet map is less than half the map size.

For variable-density and other variable-property flows, another advantage of $D>W$ is that it allows a significant interval of the domain starting from the inlet endpoint to be omitted from both filtering and the gathering of output data. This would reduce the impact of splicing-induced property discontinuities introduced near the inlet on the overall flow evolution. (Molecular transport mollifies these discontinuities before they are advected far from their points of origination.) For velocity profiles, the procedure for removing discontinuities is stated in Section 2.2.3 to involve kernels applied to the interval $[0,W]$, but more generally they should be applied to $[0,D]$, which includes the case $D=W$ assumed in that discussion.

Irrespective of the presence or absence of splicing-induced discontinuities, there is another potential artifact near both the inlet and outlet endpoints of the ODT domain. The no-flux boundary conditions applied at those endpoints (see Section 2.1) result in suppression of eddy-induced transport that might have occurred near those locations if the domain were instead unbounded, in effect creating a numerical boundary layer at each endpoint. Each occurrence of splicing refreshes these peripheral zones so the artifact might not be consequential. If it is nevertheless significant in some cases, then omitting regions near both endpoints from filtering and data gathering could be beneficial, subject to consideration of any conservation violations that might result from this.

In the general case of an unstructured mesh, CVs might have different volumes V and CV faces might have different areas $\Gamma$, so the volumetric picture is generalized accordingly. Subject to the considerations that have been noted, any choices of the domain lengths D of individual ODT domains based on user-specified criteria are allowed, so in general these domains will have different nominal cross-sections $\gamma =V/D$. Consider a volume transfer from CV 1 to CV 2 across a face of cross-section $\Gamma$, corresponding to a face-normal fluid-element displacement $\overline{v}\Delta t$. This implies splicing of an ODT interval of length $(\Gamma /{\gamma }_1)\overline{v}\Delta t$ out of CV 1, but volume preservation implies that it arrives in CV 2 with length $(\Gamma /{\gamma }_2)\overline{v}\Delta t$. Thus for ${\gamma }_1\ne {\gamma }_2$, the length of the spliced interval must be rescaled by the factor ${\gamma }_1/{\gamma }_2$ in order to preserve volume.

Because this rescaling is a volume-preserving change of the shape of the spliced interval, the resulting stretching or shrinking of the length of individual fluid elements does not imply any density changes. This rescaling modifies the length scales of flow structures represented by ODT property profiles, and hence is a model artifact. Provided that ${\gamma }_1/{\gamma }_2$ is neither very large nor very small, the artifact is not necessarily severe because it is akin to the physically based triplet-map-induced rescalings of intervals of the ODT domain. It can be avoided by assigning a fixed ratio $D/V$, hence fixed $\gamma$, to all CVs, but this might not be the optimal choice based on other considerations.

For dilatational flows the situation requires further explanation because the D value changes during the ODT advancement within each time-advancement cycle. This is a operator-splitting effect resulting from the fact that fluid-element density changes are accommodated during ODT advancement by varying the domain length. This is rectified during the pressure-projection step, which assigns CV face velocities so as to produce a net volume transfer into or out of a CV that is equal and opposite to the net volume change implied by the ODT advancement. The indicated adjustments involving the possibly CV-specific $\Gamma$ and $\gamma$ values remain valid, here assuming that splicing restores the D value corresponding to the CV volume.

If ODT were used as an SGS closure of unclosed terms in an LES rather than within the AME framework, then the LES advancement might produce CV face velocities $\overline{v}$ that do not precisely restore the prior D values of all ODT domains. This can be rectified by rescaling ODT domain lengths as needed to enforce the desired domain lengths and then using the relation $\gamma D=V$ to update $\gamma$. The adjustment should be small if there is reasonable consistency between ODT and LES advancement.

7.2. Ordering of Fluid Transfers

During each time-advancement cycle, a segment of an ODT domain is transferred across each CV face. For a CV with f faces, this involves some number $i\in [0,f]$ of inbound transfers and $f-i$ outbound transfers. (For non-dilatational flow, i cannot be 0 or f because these values would imply a change of the CV volume.) Figure 2 illustrates a portion of a notional 2D flow domain with Cartesian meshing and hence $f=4$, with $i=1$ on the left side and $i=2$ on the right side.

For $i>1$, the order in which inbound segments are attached to the ODT domain must specified. As in [36], residence-time considerations are invoked to address this. Suppose $i=2$, corresponding to two inbound segments, whose volumes are denoted $V_1$ and $V_2$, with $V_1>V_2$. If segment 2 is attached first, followed by segment 1, then segment 2 is immediately positioned a relatively substantial fraction of the distance from the inflow to the outflow endpoint of the domain, shortening its residence time prior to its eventual expulsion of its contents. Although the choice of ordering does not affect the mean residence time of fluid in the CV, the ordering 2 followed by 1 will tend to broaden the residence-time distribution relative to the choice 1 followed by 2. In the absence of evidence to the contrary, the default procedure should be to keep the residence-time distribution as narrow as possible. Accordingly, the attachments should be performed in decreasing order of segment volumes.

Now consider two outbound transfers requiring volumes $V_1>V_2$. Based on the same consideration, in what order should they be implemented? Detaching a segment of volume $V_1$ first would enable the second transfer to occur from a region relatively close to the inflow endpoint, corresponding to a relatively short residence time of the transferred segment. Therefore, the default ordering for outbound transfers is from smallest to largest segment volumes, where this choice is subject to comparison of the resulting residence-time distribution to any available empirical evidence.

7.3. SGS Contribution to Face Velocities

As noted in Section 3, conventional LES closures involve some implicit or explicit modeling of $q_{\mathrm{sgs}}$. This enables estimation of a subgrid velocity fluctuation $u_{\mathrm{sgs}}=\sqrt{2q_{\mathrm{sgs}}}$ that has been used to model backscatter by adding random increments of this order of magnitude to the filtered face velocities $\overline{v}$ [14].

This raises the conceptual question of whether ODT SGS closure requires some analog of this additional procedure in order to model backscatter. This question has physical implications with regard to the choice of the ODT domain size, which affects the statistical variability of filtered quantities as noted in Section 7.1. In the limit of infinite ODT domain size, filtered ODT property profiles converge to ensemble-average values of the filtered properties. Then in this limit, the filtered properties have the same meaning as in conventional LES, so as in LES, additional modeling would be needed in order to capture backscatter. The choice of a finite ODT domain size, which is a practical necessity, introduces finite-sample variability that can be viewed as either an error source relative to conventional LES or as a supplement to conventional LES in that it provides a model representation of backscatter as an inherent attribute of ODT.

Like any other model feature, the validity and fidelity of this backscatter representation depends on the physical realism of ODT flow advancement, but it also requires an objective basis for choosing the ODT domain size so as to obtain the physically correct degree of finite-sample variability. Backscatter modeling in conventional LES provides guidance in this regard. Because $q_{\mathrm{sgs}}$ is a diagnosed property of the ODT state, $u_{\mathrm{sgs}}$ can be evaluated. Then the ODT analog of conventional backscatter modeling is obtained when the domain size is chosen so that the finite-sample variability of $\overline{\mathbf{u}}$ values obtained for given flow conditions is of order $u_{\mathrm{sgs}}$ as determined from the ODT results. Relatively large ODT domains might be needed for multiphysics cases, e.g., with reaction-driven thermal expansion, owing to the intensity and sporadic occurrences of the associated backscatter effects. Based on the identified self-consistent choice of domain size, the resulting finite-sample variability of $\overline{\mathbf{u}}$ would be imprinted on the face velocities $\overline{\mathbf{v}}$ by the pressure projection. (This happens whether or not the procedure is optimized as described.)

In Section 2.3, the importance of generating realistic fluid-parcel trajectories through the flow domain was noted. This depends upon both the evaluation of the CV face velocities $\overline{\mathbf{v}}$ and the splicing protocol. For a splicing-based LEM closure of LES that lacked the SGS flow information provided by ODT, a particle-tracking approach involving an autocorrelated stochastic velocity contribution modeled SGS fluctuation effects on the splicing-induced fluid-parcel trajectories [29]. This improved the model representation of turbulent dispersion. With ODT closure, it is possible that $\overline{\mathbf{v}}$ fluctuations induced by ODT statistical variability will adequately capture this effect.

8. Discussion

LEM and ODT have been used in various ways for SGS closure of LES, such as for chemistry closure of reacting flow simulations, for near-wall momentum closure, and in ODTLES, which involves direction-dependent hybridization of resolved and coarse-grained momentum advancement. The present formulation is based on the viewpoint, within the low-Mach-number framework, that a broad range of turbulence phenomenology and associated microphysics can be well represented by modeling momentum advancement on a fully resolved basis, using only the coarse-grained pressure-projection step to capture global effects. The absence of time stepping of coarse-grained momentum or scalar transport equations motivates the description of the approach as autonomous microstructure evolution (AME).

As in LEM chemistry closures, one ODT domain is associated with each coarse-grained control volume. Relative to that and other previous formulations, AME includes the following novel features:

• The distribution of the available kinetic energy among velocity components during the kernel operation is consistently formulated for the general case of arbitrary external energy and momentum sources and sinks.
• Alignment of the orientation of each ODT domain with the evolving coarse-grained flow is enforced. Nevertheless, velocity components are defined in the fixed coordinate system of the physical domain.
• As in LEM SGS closure, the splicing operation that implements fluid transfers between ODT domains is governed by the coarse-grained velocity field but the transferred parcels maintain full spatial resolution of property profiles. A novel feature is the subsequent adjustment of velocity profiles that eliminates splicing-induced discontinuities subject to the applicable conservation constraints.
• A volumetric interpretation of splicing-induced fluid transfer provides flexibility with regard to the choice of ODT domain lengths while maintaining physically consistent fluid residence times on the ODT domains, and it straightforwardly generalizes to an unstructured mesh.
• A fine-grained interpretation of the splicing operation enables a physically consistent formulation of the Reynolds shear stress, enabling a prescription for diagnosing the terms of the SGS turbulent-kinetic-energy budget.
• The freedom to assign the ODT domain length, in combination with the availability of SGS TKE statistics, offers the possibility of a self-consistent quantitative approach to the representation of backscatter.
• A splicing protocol specialized to ODT domains attached to a wall at one endpoint incorporates the dominant fluid-transfer mechanism in boundary layers.
• For variable-density flow, the projection of the coarse-grained pressure gradient onto the ODT domain orientation identifies the bulk acceleration of the domain. On this basis, the Rayleigh–Taylor instability mechanism is embedded in the eddy selection and implementation procedures.
• For buoyant stratified flow in the Boussinesq approximation, a synthesis of standalone ODT formulations for horizontal and vertical domains respectively yields a formulation that is applicable for arbitrary domain orientation.

These features and the framework that encompasses them, embodied in the time-advancement cycle outlined in Section 2.3, offer the prospect of improved cost-performance characteristics relative to conventional LES and existing LEM and ODT closure formulations. This will be elaborated in sequels to the present contribution.

The adjustment of velocity profiles to eliminate splicing-induced discontinuities avoids a potential artifact, but it has been noted that this remedy is not applicable to density and other advected scalar fields. The alternative remedy of using long domains and excluding a limited neighborhood of the inlet endpoint from filtering operations applied to these property fields has been proposed. In this regard, several points are relevant. First, this artifact does not arise in constant-density flow or in buoyant stratified flow in the Boussinesq approximation, allowing a wide scope of AME application that is not subject to this artifact. Second, LEM SGS closure of mixing and chemical reactions in LES is subject to this artifact, yet its performance is satisfactory for a variety of combustion applications [2].

Finally there is the mitigation strategy of coarsening the mesh. Regardless of the ODT domain sizes in individual CVs, for a given mesh a fluid element must pass through the same number of CVs in order to traverse the flow and hence is spliced that number of times. A reduction in the number of splicings can only be achieved by coarsening the mesh.

This raises the question of the mesh resolution requirement in AME. AME time advancement is structured to involve only minimal large-scale influence through the pressure-projection operation, where the justification is that ODT is formulated to adequately represent flow phenomenology at all scales below those dominated by case-specific 3D large-scale forcings and constraints such as inflows, outflows, and boundary conditions. Therefore, a coarse very-large-eddy-simulation (VLES) type of mesh is not only viable but is central to the AME concept. The efficacy of coarse meshing has been shown in applications of ODTLES [11,12].

Coarse meshing not only mitigates any splicing artifacts but also reduces computational cost. For a flow that would require $N^3$ DNS CVs, AME with $M^3$ CVs where $M<N$ would require order $N/M$ fully resolved cells per ODT domain for a total of $M^{2}N$ ODT cells. Assuming that computation time scales with the total number of ODT cells, this corresponds to a cost-reduction factor ${(M/N)}^2$ relative to DNS. This cost reduction could somewhat or entirely compensate for the cost of ODT relative to other SGS closures, depending on the case. Experience with ODTLES has indicated that the parallel-processing efficiency of SGS ODT results in acceptable turnaround times on high-performance computers [13].

Apart from the discontinuity artifact, splicing is a numerical form of ‘mesomixing’ (meaning at subgrid scale, but not at the molecular-mixing microscale) in the same sense that the ODT eddy event is a physically based mesomixing operation. This aspect of splicing is analogous to numerical dissipation in LES that can dominate the physical based dissipation if it is not carefully managed. Mesh coarsening reduces the number of times that a fluid element is spliced as it traverses the flow but the number of physical mesomixing events during the traversal, whether modeled as eddy events or resolved on the 3D mesh, is unchanged. Thus, mesh coarsening has the additional advantaged of reducing numerical mesomixing relative to physically based mesomixing, thus enhancing the physical fidelity of the numerical implementation.

High fidelity with VLES-type meshing is plausible for simple flows such as the canonical free-shear flows. In the near-wall regions of boundary layers and confined flows, the demonstrated performance of ODT as a standalone and an SGS near-wall model [10,11,12,16,18,21] provides assurance that the ODT near-wall formulation for AME described in Section 4 will likewise provide suitable closure on a coarse mesh.

Unsteady boundary layers are subject to local flow separation and re-attachment. Related features such as recirculation zones arise in bulk flow. The compatibility of the present formulation with unstructured meshes provides the flexibility to refine the mesh adaptively so as to increase spatial resolution when and where it is needed. Where the local refinement leaves an insufficient range of subgrid scales to benefit from the multiscale ODT treatment, the splicing protocol can be retained but box filtering can be applied after each splicing operation so as to implement the no-model closure described in Section 2.2.1.

An approach to adaptive meshing termed the arbitrary Lagrangian Eulerian (ALE) method [37] employs operator splitting consisting of Lagrangian implementation of advection and Eulerian implementation of molecular transport. The mesh adaption advects CV faces, involving no property fluxes across faces. This overcomes the advective CFL time-step limitation although there remains a weaker time-step bound that prevents topologically forbidden mesh structures and related anomalies. After mesh adaption, the flow state is mapped onto an improved mesh structure, causing redistribution of fluid volumes among CVs. This implements the physically prescribed molecular-transport fluxes across CV faces.

In AME with ALE meshing, splicing is decomposed into two operations. First, the ODT domains are displaced in accordance with the advection of the CVs containing them. Second, the mesh remapping prescribes transfers of domain segments across CV faces. This implementation reduces both splicing frequency and volume fluxes per splicing occurrence, thus mitigating numerical artifacts as in conventional ALE.

SGS ODT within the AME framework has been emphasized because this is convenient for explaining the closure approach without reference to the particulars of LES implementation. Addressing SGS ODT within an LES framework in full detail requires consideration of the interactions between the numerical and physical modeling implementations in a particular LES code, which is beyond the present scope. Nevertheless it is pertinent to note that a basic requirement for closure of LES (and likewise AME, as shown in Appendix B) is evaluation of the filtered Reynolds stresses. In Section 3 a formalism for their evaluation based on analysis of the SGS TKE budget is developed and is shown to be broadly consistent with the Navier–Stokes equation-based LES formalism.

Another potential application of ODT SGS closure is motivated by LEM chemistry closure of LES [2]. The LES provides LEM with CV face velocities that govern splicing and a parametric representation of SGS turbulence that controls LEM SGS advection. The only LEM input to LES flow advancement is the filtered divergence, deduced from volume creation caused by thermal expansion, that is used to solve the coarse-grained pressure Poisson equation. ODT could be used in this formulation instead of LEM. ODT would receive additional input from the LES in the form of pressure-gradient forcing of the ODT momentum equation, but would provide the same output as LEM to the LES flow advancement. The ODT filtered flow field would then be at best roughly consistent with the LES flow field, but this should be sufficient for the intended chemistry closure. The potential benefit of the use of ODT instead of LEM is relaxation of the restriction to standard inertial-range similarity scalings, instead allowing dynamical effects such as SGS baroclinic vorticity, the Darrieus–Landau instability, and the kinematical effects of thermally induced viscosity fluctuations to be captured. Importantly, this substitution is not simply a matter of replacing LEM by ODT in existing formulations. As explained in Section 2.2.3, it also requires the application of a specialized splicing protocol to ODT velocity profiles.

The present restriction of AME to low-Mach-number flow is a consequence of the exclusion of coarse-grained momentum advancement. The relaxation of this restriction by extending the scope of ODT SGS implementation to LES closure leads to consideration of closure of compressible LES formulations using low-Mach-number ODT as described here or its extension to compressible flow [32]. In applications to supersonic inert mixing [38] and to detonation phenomena [39,40], LEM has been used for subgrid mixing/reaction closure in which high-speed effects are communicated by the LES thermodynamic inputs to LEM. This is a suitable strategy for some purposes, but ODT SGS closure, especially if compressible, would extend the scope of application of LES to high-speed flows. The dynamically based ODT representation of SGS kinetic energy (Section 3) would be a pertinent input to the LES energy equation.

Thus, the novel flow-advancement features that have been described have potential benefits beyond AME per se. Similarly, AME can be supplemented by coarse-grained advancement in particular situations in which it is advantageous, as illustrated by a previous approach to ODT-based near-wall closure that is mentioned in Section 4. Another example, noted in Section 5.1, is buoyant stratified flow at geophysical scales, for which SGS ODT can serve as a scale-range extension relative to conventional closures, but cannot affordably resolve the physical microscales and accordingly requires its own SGS closure. Although coarse-grained in this sense, this formulation preserves the defining AME feature that there is no time advancement of a momentum equation on the filtered 3D mesh.

Details of the numerical implementation of AME will be presented in a future publication, accompanied by applications to canonical flows for demonstration and validation purposes. The relative performance of alternative strategies for model formulation or numerical implementation that have been noted will be evaluated, involving tradeoffs with regard to physical soundness, numerical accuracy, and computational cost. For clarity of interpretation, the initial test cases will be simple and in some instances within the scope of the capabilities of conventional LES, ODTLES, or standalone ODT. These steps are targeted on the principal goal of capturing the configuration-specific coupled dynamics of microphysical processes and high-intensity turbulence that are difficult to investigate using presently available methods.