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 ﬂows 1 because the range of relevant motions is very broad and encompassing all of them is unaffordable. 2 The typical compromise is to omit the smallest motions, replacing them by estimates of their effects, 3 but this is not always adequate. Instead, it is proposed to fully resolve the ﬂow advancement along 4 selected line segments within the ﬂow volume and to couple these segments so that the resulting 5 formulation adequately represents all scales of motion. The modeling needed to implement this is 6 described and the overall approach, termed autonomous microstructure evolution, is assessed 7 conceptually and with reference to demonstrated capabilities of related methods. 8 Abstract: A multiscale modeling concept for numerical simulation of multiphysics turbulent 9 ﬂow utilizing map-based advection is described. The approach is outlined with emphasis on its 10 theoretical foundations and physical interpretations in order to establish the context for subsequent 11 presentation of the associated numerical algorithms and the results of validation studies. The 12 model formulation is a synthesis of existing methods, modiﬁed and extended in order to obtain a 13 qualitatively new capability. The salient feature of the approach is that time advancement of the 14 ﬂow is fully resolved both spatially and temporally, albeit with modeled advancement processes 15 restricted to one spatial dimension. This one-dimensional advancement is the basis of a bottom-up 16 modeling approach in which three-dimensional space is discretized into under-resolved mesh cells, 17 each of which contains an instantiation of the modeled one-dimensional advancement. Filtering 18 is done only to provide inputs to a pressure correction that enforces continuity and to obtain 19 mesh-scale-ﬁltered outputs if desired. The one-dimensional advancement, the pressure correction, 20 and coupling of one-dimensional instantiations using a Lagrangian implementation of mesh- 21 resolved volume ﬂuxes is sufﬁcient to advance the three-dimensional ﬂow without time advancing 22 coarse-grained equations, a feature that motivates the designation of the approach as autonomous 23 microscale evolution (AME). In this sense, the one-dimensional treatment is not a closure because 24 there are no unclosed terms to evaluate. However, the approach is additionally suitable for use as 25 a subgrid-scale closure of existing large-eddy-simulation methods. The potential capabilities and 26 limitations of both of these implementations of the approach are assessed conceptually and with 27 reference to demonstrated capabilities of related methods. 28


Introduction
To extend the range of applicability of numerical simulation of turbulent flows 31 beyond the limits of affordable direct numerical simulation (DNS), under-resolved 32 simulations with physically or algorithmically based subgrid-scale (SGS) closures are 33 widely used, typically within the large-eddy-simulation (LES) framework. For constant- 34 property flow, the main role of the SGS closure of LES is to dissipate kinetic energy at 35 the physically correct rate based on the flow state. A secondary process that is important 36 eterizations, but it is less costly than DNS because this fine resolution is restricted to 48 representative lines of sight within individual LES control volumes (CVs). Cost is further 49 mitigated to the extent that the fidelity of the SGS closure allows coarsening of the LES 50 mesh while maintaining the needed level of accuracy. 51 The method used for SGS closure involves one-dimensional turbulence (ODT), a 52 form of map-based stochastic advection that has been used previously for SGS closure 53 of LES as well as in various standalone applications, as cited below. Here, an ODT-based 54 SGS closure is formulated that draws upon previously reported formulations but also 55 introduces various novel elements whose details and physical interpretations are the 56 main focus of the present contribution.  The velocity profiles u i are distinct from the velocity fieldv i that governs coarse-117 grained advection in the LES sense. Here, the latter is the set of direction-i LES-prescribed 118 CV face-normal velocities, as in a typical staggered numerical scheme. The overbar 119 notation indicates box filtering. The box-filtering interpretation ofv i is explained in Sect. 120 3.

121
The ODT domains are Lagrangian objects whose orientations need to be defined 122 only for particular applications such as buoyant stratified flow and particle-laden flow. 123 For those cases, the domain within a given CV is assumed to be aligned with the velocity  The ODT state is subject to two advancement processes. One is viscous advancement governed by In the last term, ρ is density and p ,i is the index-i directional derivative of the LES- Equation (1) includes LES-level pressure effects but not SGS pressure fluctuations. The latter and related effects are modeled by the second eddy operation. Let K(x ) = x − x(x ) denote the map-induced displacement of the location x that is mapped to x , here expressed as the inverse-map function x(x ), and let J = |K|. 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 (x). Then the complete eddy event is the transformation where the coefficients c i and b i are chosen so as to change component momenta and 150 kinetic energies subject to applicable conservation laws. The criteria for specifying these 151 changes and the determination of c i and b i on that basis are explained in Sects. 6.1 and 152 6.2.

153
K, which is termed a kernel, is referenced here to the selected eddy interval [x 0 , x 0 + 154 l], but its mathematical definition is applicable to any x interval, which need not be an 155 eddy interval, and likewise for J. In what follows, K and J are generic tools for modifying 156 property profiles subject to specified physical requirements.

157
The main physical content of ODT is the stochastic process by which eddy intervals 158 and times of occurrence are determined. The profiles u i of the velocity components 159 within the various possible eddy intervals are the inputs to that determination, and the 160 kernels play a key role in the specification of the eddy rate distribution λ(x 0 , l; t) based on 161 that input. λ has units time −1 length −2 such that λ(x 0 , l; t) dx dl is the rate of occurrence 162 of eddies whose left boundaries are within the range [x 0 , x 0 + dx] and whose lengths 163 are within the range [l, l + dl], where the argument t reflects the dependence of λ on the 164 system state within [x 0 , x 0 + l] at time t. The specification of λ(x 0 , l; t) and the way that 165 it is used to sample eddy events are described in Sect. 6 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 velocitiesv are identical, where the coordinate index is omitted for this 1D illustration. The discretized scalar field has scalar value c j in cell j. For the purpose of streaming the scalar along the 1D domain with velocityv > 0, treat the scalar c as a function of the continuum coordinate x. Because c is uniform within each cell, it is piecewise constant in x. Streaming for a time interval dt < W/v produces a shifted scalar profile c * (x) = c(x −vdt). Box filtering within each cell gives, after rearrangement, (c * j − c j )/dt +v(c j − c j−1 )/W = 0, which is the lowest-order discrete approximation of the Eulerian 1D advection equation Although the spatial-continuum picture reduces to an Eulerian numerical algorithm, 193 physically it prescribes Lagrangian time advancement as follows:

3.
Box filter each domain c(x; j).

199
One consequence of the error inherent in the numerical approximation is that c is 200 subject to numerical dissipation that is asymptotically diffusive over a long time interval 201 with numerical diffusivity of order Wv. Owing to time advancement on all ODT domains prior to splicing, the profile u(x) after 258 splicing is in general discontinuous at x * .

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

263
One approach is to box filter discontinuous velocity profiles. This conserves mo-264 mentum but dissipates kinetic energy, with no kinetic energy remaining to be dissipated 265 after that, as indicated by the fact that the kinetic-energy dissipation rate is then zero.

266
This is analogous to the lowest order 'production equals dissipation' level of turbulence 267 modeling, or it alternatively can be seen as analogous to implicit LES (ILES), in which 268 the numerical algorithm provides the closure [6].

269
The other approach is to use the K and J kernels (see Sect. nonuniformities, implying production of subgrid-scale turbulent kinetic energy (TKE).

295
As noted, the first approach dissipates the TKE immediately. However, the adjustment 296 after ODT splicing is globally (within each ODT domain) conservative and therefore is 297 operationally an advective-transport mechanism rather than a dissipative process. This 298 is consistent with the foundational principle that the effects of viscous and any other 299 molecular transport processes on the ODT-level system states are implemented only 300 through time advancement of molecular-transport equations, in the present instance Eq.

302
The practical consequence is that TKE dissipation in the ODT domains is the local imbalances between production and dissipation. Thus, these imbalances arise in 307 a manner that is more closely linked to the governing physics than is the treatment 308 provided by, e.g., a parameterized SGS kinetic-energy equation.

309
Advected scalars such as species mass fractions are not amenable to adjustments 310 that remove scalar discontinuities owing to the stricter conservation laws that must be 311 enforced locally in order to assure chemical realizability. Therefore scalar discontinu-

3.
The updated face-normal velocities are interpolated to update the quantitiesū j .

4.
Advection of the velocitiesū j as prescribed by the CV face-normal velocities is 329 implemented as a 3D generalization of the discretized advection of the scalar c(x; j) 330 in Sect. 2.2.

331
The numerical diffusivity associated with the advection operation, noted below Eq. (3),

332
is not necessarily an accurate representation of SGS advective transport, which is why 333 the algorithm is termed an ILES analog rather than an ILES. In any case, this sequence 334 of operations is LES-like advancement of the momentum equation although it is not a 335 conventional LES procedure.

336
The analogous time-advancement cycle for ODT-based SGS closure is as follows:  an approach that was shown to work well in [9] is applied. Namely, instead of 362 incrementing u i (x) uniformly over x, the correction 2(x/D)C p is added to u i (x).

363
This yields a consistentū i value and obeys the no-slip condition at the wall location 364 x = 0.

366
A key feature of the approach is that steps 2-4 fulfill the requirements for coarse-  This is a correction that enforces a constraint rather than a time-stepping procedure, so it 372 is fair to say that this formulation involves no time advancement of a filtered momentum identify the classes of problems for which the respective approaches are advantageous.
As in ODTLES, AME distinguishes between velocity profiles u that represent the 383 flow states on the ODT domains and velocitiesv that govern coarse-grained advection.

384
The latter are in this respect auxiliary variables. The coarse-grained flow state represented 385 by the set of velocitiesū is likewise auxiliary from the ODT perspective, but from the 386 coarse-grained perspective, theū velocities constitute the desired flow solution. In LES 387 implemented on a staggered mesh, the distinction betweenv andū is purely algorithmic, 388 while in AME there is a further distinction becauseū is a coarse-grained approximation 389 of the detailed physical representation of the evolving flow state produced by ODT.

390
A nuance in this regard is that there is a fine-grained representation of the advection 391 induced in AME by the velocitiesv, as explained in Sect. 3. The analogy between AME and LES that is described in Sect.

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

401
Splicing is the AME implementation of coarse-grained advection. Splicing is CV-

410
On an ODT domain [0, D], suppose that a sub-interval [x * , x * + ∆] is designated to be transferred to the ODT domain in the CV to its right, corresponding to positivē 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 Sect. 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-∆ interval. The advecting velocity governing the displacement of that interval is denoted v * . Accordingly, the profile of the advecting velocity This shows that the face velocity can be interpreted as a box-filtered velocity, thus 411 motivating the notationv that has been adopted.

412
Similarly, the fluctuation of the advected velocity is u = u −ū. Spatial averaging of the product of u and v then gives the contribution of the outbound transfer of the size-∆ interval to the box-filtered Reynolds stress, whereū ∆ denotes the average of u over [x * , x * + ∆] and x * + ∆ is the outflow endpoint 413 D of the ODT domain.

414
The outbound fluid is replaced by inbound fluid occupying a size-∆ interval for 415 which x * = 0 owing to its attachment to the inflow endpoint. The ODT state in the CV 416 that supplies this fluid is denoted u − . The filtered value of the u − interval that is trans-417 ferred rightward, denotedū ∆ − , is in general different fromū ∆ , where the notation refers 418 to the distinct size-∆ intervals containing the inbound and outbound fluid respectively.

419
The inbound and outbound transfers convert the original profile u into a new profile U.  The inbound fluid enters the CV through the face opposite to the outflow considered thus far. As noted,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 where v − is specified by Eq. (5) with u replaced by u − .

423
The new total (filtered plus SGS) kinetic energy is 1 2 U 2 , which is the sum of the filtered and SGS kinetic energies, 1 2 (Ū) 2 and Q 2 sgs = 1 2 U 2 respectively, where U = U −Ū. Similarly, the initial SGS kinetic energy is q 2 sgs = 1 2 u 2 , where u = u −ū. Then the splicing-induced SGS TKE change δq 2 sgs = Q 2 sgs − q 2 sgs is expressed as The approach is to evaluate C + S = δq 2 sgs − A − T using the appropriate specializations of the terms on the right-hand size and then to compare the outcome to the conventional expression for C and S. Extending interval-restricted averaging to arbitrary powers p of u and u − , the identities and will be useful in what follows. As before, each size-∆ averaging interval begins at  Evaluation of δq 2 sgs = 1 2 U 2 − u 2 − (U) 2 + (u) 2 using Eqs. (9) and (10) gives In accordance with the relation ∆ =v dt, this change is attributed to the time-advancement step dt. For simplicity, Eq. (11) is specialized to vanishing dt, so ∆/D 1. Noting that u D−∆ converges toū in this limit, the leading-order result is obtained.

436
In the conventional decomposition of ∂ t q 2 u 2 ) in the present reduced notation. In Eq. (8), a difference of q 2 sgs values separated in 438 time by dt is analyzed, so the conventional terms are multiplied by dt, which will lead to 439 cancellation ofv based on ∆ =v dt. Additionally, each spatial derivative is replaced by 440 1/W times a difference between the CV in which the q 2 sgs is analyzed and the CV to its 441 left, whose velocity profile has been denoted u − . Finally D is changed to W owing to 442 the specialization to D = W.

443
On this basis, and are obtained, giving In Eq. (8) with V omitted, substitution of results for the terms that have been evaluated gives Rewriting this as each term is the product of a discretized derivative of filtered u velocity and a fluxtype factor involving only filtered fields. This result is compared to the conventional where the spatial discretization reflects the evaluation of the Reynolds stresses as domain averages rather than point values at domain endpoints. This gives C Using Eqs. (6) and (7) to evaluate the Reynolds stresses, is obtained.

444
Equation (17)  The use of closure models that resolve the processes governing q 2 sgs evolution 458 eliminates the need for modeled equations for q 2 sgs that are otherwise used for closure 459 in various LES formulations [13]. Within the AME framework, the evaluation of q 2 sgs 460 budget terms is for diagnostic purposes rather than for LES closure and serves to  each domain is associated with one of those faces and is deemed to be normal to that face.

477
The domains associated with a given CV are treated as independent during ODT time 478 advancement. As in [9], each domain is attached at one of its endpoints to its associated by quasi-steady-quasi-homogeneous theory and its validations [18], which indicate 496 that small-scale near-wall motions (resolved by ODT in AME) depend primarily on the 497 larger-scale local flow state (controlled by splicing in AME).

498
In contrast, the method in [9,17] includes Eulerian implementation of direct ad-499 vective transport, at each height above the wall, between parallel ODT domains in 500 neighboring wall CVs. This is under-resolved time advancement of advection, hence out-501 side the scope of the AME concept, but its adoption could be advantageous, depending 502 on the application.

Density variation without dilatation 583
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 [30].) 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 Sect. 6. Also, Eq. (1) is generalized to the variable-density form where µ is the dynamic viscosity.
Here, a dilatational velocity profile u d (x) is introduced that is based on 1D continuity generic argument x, and hence the notation c (x), is used in what follows.

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

688
The second operation, which is applied only to velocity components, is the kernel piecewise-uniform adjustments of velocity profiles to remove discontinuities.

703
These kernel effects and the applicable conservation laws constrain the coefficients a i and b i in Eq. (2) as follows. For variable-density flow, momentum conservation implies Integration over dV reflects the volume interpretation of the ODT domain based on 704 the cross-sectional area γ introduced in Sect. 7.1. The integrals extend over the entire 705 domain, but the integrands are zero outside the eddy interval so they can be equivalently 706 restricted to that x range. The argument x is henceforth omitted. Here and below, a 707 previous constant-density treatment [33] that introduced momentum and energy sources 708 (or sinks) S i and S E respectively is extended to variable-density flow by bringing ρ into 709 all integrands. This results in minor generalizations of the algebraic solutions in [33].

710
Details that are readily inferred from that treatment are omitted here.

711
Because the triplet map conserves total momentum, u i can be substituted for u i in the integrand on the right-hand side. This cancels a term on the left-hand side, reducing Eq. (21) to This relation can be substituted for b i in equations that follow, reducing them to relations The eddy-induced change of the domain-integrated kinetic energy of velocity component i is Conservation of total kinetic energy by the triplet map allows u i to be substituted for u i in the second integral, resulting in a cancellation of terms that gives Substitution of Eq. (22) for b i yields a quadratic equation for c i whose solution straight-715 forwardly generalizes Eq. (11) of [33].  If the discriminant is positive, then the quadratic equation has two real solutions.

732
The solution branch is chosen on the basis of the functional dependence of c i on ∆E i .

733
As ∆E i monotonically approaches zero starting from its assigned value, one of the Rotational eddy motion can redistribute kinetic energy among velocity components.

739
Existing formulations of the kernel operation, such as in [33], take this into account.

740
Depending on the application, the net kinetic-energy change S E might be specified on a where subscripts j and k refer to the coordinate indices not equal to i. If S E = 0, 765 this is the usual energy redistribution for α = 2 3 , where α denotes the coefficient of −Q i .

766
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 where χ is a 'memory' parameter. This recovers Q * i = H 3 for the 'no-memory' value 767 χ = 0. Because Q i and H are non-negative and Q i ≤ Q, Q * i is also non-negative provided 768 that χ ∈ − 1 2 , 1 . The definition of available energy precludes negative Q * i , so a value of 769 χ outside this range could yield a mathematical inconsistency.

770
The definition of α as originally stated [4] does not uniquely specify the treatment The model adopted for τ [34], re-expressed in present notation, is where 810 • KK ≡ K 2 dV,ρ KK ≡ 1 KK ρK 2 dV, and V is the eddy volume. to the available energy so as to increase the amplitudes of the applied kernels, thereby 856 promoting eddy occurrences, has been incorporated into combustion simulations as a 857 representation of the Darrieus-Landau instability [15,32]. SGS ODT in AME can similarly 858 capture this phenomenology.  insufficient sampling [29]. That formulation did not involve splicing, but D = W was 891 found to be sufficient in a cloud LES study using LEM SGS closure with splicing [28]. found that eddy events extending to size L max = 4W were required for this purpose [10].

903
This is qualitatively reasonable because the typical displacement of a fluid element by a 904 triplet map is less than half the map size.   that has been used to model backscatter by adding random increments of this order of 975 magnitude to the filtered face velocitiesv [13]. here or its extension to compressible flow [31]. In applications to supersonic inert 1143 mixing [37] and to detonation phenomena [38,39], LEM has been used for subgrid 1144 mixing/reaction closure in which high-speed effects are communicated by the LES 1145 thermodynamic inputs to LEM. This is a suitable strategy for some purposes, but ODT On the right-hand side, the first term represents Liouville evolution of the PDF due prime denotes x-derivative and a functional derivative is applied to the quantity in 1172 braces. The second term represents the loss in probability for the configuration u due 1173 to eddy events, which occur at the rate R u and lead to configurations different from u.

1174
The third term represents the gain in probability for the configuration u due to eddy