Turbulence modulation in particle-laden stationary homogeneous isotropic turbulence using one-dimensional turbulence

Turbulence modulation in particle-laden stationary homogeneous isotropic turbulence is investigated using one-dimensional turbulence (ODT), a low-dimensional stochastic flow simulation model. For this purpose, ODT is extended in two ways. First, a forcing scheme that maintains statistical stationarity is introduced. Unlike direct numerical simulation (DNS) of forced turbulence, the ODT framework accommodates forcing that is not directly coupled to the momentum equation. For given forcing the ODT energy dissipation rate is therefore the same in particle-laden cases as in the corresponding single-phase reference case. Second, previously implemented one-way-coupled particle phenomenology is extended to two-way coupling using the general ODT methodology for flow modulation through interaction with any specified energy and momentum sources and sinks. As in a DNS comparison case for Reλ = 70, turbulence modulation is diagnosed primarily on the basis of the fluid-phase kinetic-energy spectrum. Because ODT involves subprocesses with straightforward physical interpretations, the ODT mechanisms of particle-induced turbulence modulation are clearly identified and they are plausibly relevant to particleladen Navier-Stokes turbulence. ODT results for the ratio of particle-phase and fluid-phase kinetic energies as a function of particle Stokes number and mass loading are reported for the purpose of testing these predictions in the future when these quantities are evaluated experimentally or using DNS.

In this study an alternative approach, one-dimensional turbulence (ODT), is used to simulate particle-laden stationary HIT. ODT is a stochastic approach resolving the full range of length and time scales on a one-dimensional (1D) domain and is therefore a cost-efficient method to reach higher ranges of Reynolds numbers than DNS. After its first introduction [22] and extension [23], comparisons to DNS studies and experimental results for many geometrically simple flows, such as boundary layers and jets with large mean property gradients in one direction, have shown that ODT has useful predictive capabilities for such flows. In a related multiphase application Movaghar et al. [24,25] showed the ability of ODT to model primary breakup in a turbulent jet application. Schmidt et al. [26] and Sun et al. [27,28] extended the ODT model to predict fluid-particle interactions. Fistler et al. [29] introduced a new formulation of those interactions that is explained in detail here. The two latter studies showed that ODT was able to capture two-way coupling effects. Due to its cost effectiveness ODT is a promising tool to investigate turbulence modulation for a large range of parameter space.
This paper evaluates the capability of particle-laden ODT to capture turbulence modulation in particle-laden stationary HIT, a flow for which DNS results that are particularly suitable for model validation are available. The numerical methods used to capture the fluid phase, including the forcing scheme to achieve stationarity of the TKE level, are described in detail in Sec. II. In Sec. III, the governing equations of the particle phase and the modeling of its interaction with the fluid phase are discussed. In Sec. IV, results are compared with existing DNS results, and predictions for which there are not yet any comparison data are presented. Finally, in Sec. V, the results are summarized and discussed.

II. FLUID PHASE
This section describes the concept of the general ODT model and the forcing scheme. The main focus here is to provide enough details about ODT to introduce a modified version of the two-way coupling mechanism (Sec. II C) and a novel approach to large-scale forcing (Sec. II D).

A. Model formulation and evolution processes
ODT is a stochastic model to simulate turbulent flows on a one-dimensional domain that is oriented in the direction of the largest mean velocity gradients. In the present case of HIT the domain orientation is not important as no dominant direction exists. The ODT line domain can be interpreted as a line of sight through a 3D flow field, although importantly, the ODT formulation used here evolves as a closed system rather than as open system with inflows and outflows.
The basic ODT implementation used in this study was described in detail in Ref. [30]. In the following we will give a brief summary of the fundamental ideas and constructs of ODT. ODT is a numerical method to simulate realizations of turbulent flows. In ODT, time advancement consists of viscous transport to implement the dissipation of turbulent kinetic energy (TKE) and a stochastic model to mimic turbulence effects on the velocity profile along a one-dimensional line for the full range in time and space of the turbulent cascade. The latter consists of so-called eddy events in which the fluid property profile is rearranged by implementing triplet maps, defined below, in a manner consistent with turbulence scaling laws, followed by momentum changes within the individual fluid parcels comprising the eddy. The latter step represents pressure-fluctuation effects. Time advancement of viscous transport is paused at the times of occurrence of eddy events to allow eddy implementation, which is instantaneous, before viscous advancement resumes.
For the present application, a TKE forcing scheme termed a kernel event is used to achieve statistical stationarity. It implements only the second part of an eddy event. The length within the domain over which each of these events is applied is fixed, but the location is sampled uniformly within the domain and occurrence times are sampled with a fixed occurrence frequency. This is described in detail in Sec. II D. 044308-3

B. Viscous advancement
The governing equations for mass and momentum are based on the Reynolds transport theorem in incompressible form and result in the continuity and momentum transport equations, respectively. Because the mean velocity is zero in HIT, all advection is deemed to be turbulent and therefore is represented by the eddy events, so the advancement equations below do not include advective terms. In discrete form for individual control volumes with volume V c and with quantities stored at the cell centers, the continuity and momentum equations are given as and Here, u i is the ith fluid velocity component and ρ is the fluid density, which is henceforth assumed to be constant. The reason for neglecting the pressure term here is discussed in detail in Kerstein et al. [23]. To preserve the dimensionality of mass and momentum on the one-dimensional domain the cell volume is defined as the product of a cell cross-section area A c and the cell length y on the ODT line. This determination can be arbitrary, but will be important when treating extensive quantities, e.g., finite particle numbers. τ i is the ith viscous stress component evaluated on the east or west cell face and subscripted with e or w, respectively. (The kinematic viscosity ν is thus introduced into the model.) S p,i represents the momentum exchange between particle and fluid phases that accounts for two-way coupling.

C. Eddy events
In ODT simulations, turbulence, which can be seen as a three-dimensional vortex stretching process [31], is modeled through eddy events. These involve a map applied to the flow property profiles within a sampled eddy region, which is defined by a location y 0 and an eddy size l that specify the eddy range [y 0 , y 0 + l]. This model consists of two key components, the mapping method, which is called a triplet map, and a model to define the rate, locations, and sizes of eddy events [22]. The triplet map function compresses the original profile by a factor of three over the eddy region and three copies are filled in. To ensure continuity of the profile the second copy in the middle is spatially inverted (see Fig. 1). This process is represented mathematically as a fluid at location f (y) being mapped to location y, where f (y) is given as The triplet map function conserves all quantities, increases scalar gradients and decreases length scales, which reflects the behavior of notional turbulent eddies [27]. An essential feature of turbulence is the phenomenon of return to isotropy, which requires on the ODT modeling side a redistribution of TKE among the velocity components. Additionally, momentum and energy exchange between phases (as in the present application) and between kinetic and gravitational potential energy (as in buoyant stratified flow applications) has to be ensured. This is achieved by introducing a kernel transformation that follows the mapping operation. The complete eddy event is denoted symbolically as where u i is the velocity in the ith direction before and u TM i is the velocity after the mapping process. The kernel K (y) is defined as the fluid displacement profile y − f (y) that is induced by the triplet map and integrates to zero over the eddy region, so, for constant density, it does not induce an overall momentum change. J (y) is the absolute value of K (y) and so it does not integrate to zero over the eddy region. Thus, it forces momentum change of the profiles if its coefficient b i is nonzero. c i scales the amplitude of K (y). The effect of kernels is illustrated in Fig. 1, indicating that the eddy-integrated momentum is modified by the J kernel but not by the K kernel.
Momentum conservation requires that V e is the volume of the eddy region on the ODT line. S i represents the sum of component-i interphase momentum penalties over the particles within the eddy volume. The evaluation of S i is explained in Sec. III D. The triplet map conserves momentum, so ρ V e u i dV = ρ V e u TM i dV , and with ρ V e c i KdV = 0, Eq. (5) can be solved for b i as Similarly, the energy balance is used to evaluate c i . The map-induced change of the domainintegrated kinetic energy of velocity component i is where E i includes any sink and source terms, which in this study are the energy transfers S E ,i from the particle phase. Due to energy conservation during triplet mapping, V e u 2 i dV = V e (u TM i ) 2 dV and it follows that Reordering in powers of c i and inserting Eq. (6) for b i yields For clarity, Eq. (9) is expressed in a shorter form as The solution for c i is The choice of solution branch is explained in Ref. [23]. E i specifies the redistribution of the component energies. For this redistribution, it is useful to evaluate the maximal energy Q i that is available to subtract. This is determined by maximizing the right-hand side of Eq. (10) with respect to c i , giving Given the maximal available energy change, Q i , the energy change for a given component is computed as S E ,i represents the sum of interphase energy penalties. α is an additional model parameter with its allowed range 0 α 1. It determines the extent of the redistribution of energy among the velocity components. α is chosen to be 2/3 in accordance with a previously explained physical interpretation [23]. As a next step, it is important to define the eddy event rate λ e , which is parametrized by the eddy origin y 0 and length l, and depends on the current flow state within that interval. This rate is modeled 044308-6 by using dimensional arguments, giving where τ e is an event-specific time scale. To determine the eddy timescale τ e we use the available kinetic energy in the sampled size-l region. Based on the scaling assumption for kinetic energy E kin ∼ 1 2 mv 2 ∼ 1 2 ρV e l 2 τ 2 e , the eddy time scale is modeled as where E kin = KK V e l 2 i Q i , V e is the volume of the sampled eddy region, and KK = V e K 2 dV . The viscous penalty energy is given as E vp = V e 2l 2 μ 2 ρ and ρ and μ are eddy volume averages of density and dynamic viscosity, respectively. In this study, they will be constant. C is an adjustable eddy rate parameter and scales the overall eddy event frequency. Z is the viscous penalty parameter, which suppresses nonphysical small eddies.
The integral of λ e over y 0 and l defines the rate e of all eddies. Then the instantaneous joint probability density function (PDF) of eddy size and location is given as We assume that the occurrence of eddies follows a Poisson process in time with a mean rate e , i.e., P( t ) = e exp(− e t ). Here, P denotes the PDF of the time between successive occurrences. Technically this is solved by oversampling, i.e., generation of candidate eddies at a much higher rate than specified by the model, and thinning of the Poisson process with an acceptance-rejection method. For details we refer to Ref. [30].

D. Kernel events
The strategy for simulating stationary forced homogeneous isotropic turbulence is to keep the TKE constant in the system by performing an average of one external injection of TKE T KE per time interval 1/λ ke . λ ke is the mean occurrence rate of the TKE injection. Each injection, termed a kernel event, implements a kernel function like that used during eddy events describe above specifying so as to supply the needed TKE increment T KE, which is determined below. Based on a fixed domain length D and a fixed length l of the kernel event, chosen to match the integral length scale, the time interval 1/λ ke is specified as where T 11 is the large-eddy turnover time defined below in Eq. (33). This assures that a given location experiences TKE injection once on average per large-eddy turnover time, consistent with the phenomenology of large-scale forcing. This phenomenology is governed by the dimensional expression for the production rate P, The production rate value P, the length l, and time scale T 11 of the forcing are required as inputs akin to the DNS forcing schemes of Eswaran and Pope [19] and Mallouppas et al. [7]. Thus, C ke is uniquely defined. Using Eqs. (17) and (18) the interval can be written as The TKE increment T KE that is required so that the average rate of energy injection at any given location matches the energy production rate P is then where A c is the cross-sectional area of each mesh cell as described before. Inserting Eq. (19) in Eq. (20) it follows that Application of the kernel K (y) within the event interval l can change the TKE in that interval by any desired amount without changing the total momentum in that interval, so this procedure is suitable for implementing the external TKE injection. Accordingly, the TKE injected by a kernel applied to the ith velocity component within the kernel-event volume V k is specified as thus isotropically partitioning the injected energy among the three velocity components. The component-i kernel coefficient d i must be assigned so as to inject the specified TKE increment. The term V k K (y) 2 dV is analytically evaluated as 4 27 l 3 A c [32]. Solving Eq. (22) for d i gives where The occurrence times of kernel events are samples from a Poisson process. Therefore, as in Sec. II C the sampled time intervals t ke between randomly occurring, independent events with a mean rate of occurrence λ ke are exponentially distributed. Due to the characteristics of forced HIT the energy production rate P and dissipation rate are equal over time and in the following only mentioned as energy dissipation rate .

III. PARTICLE PHASE
The particle phase treatment described in the following is an extension of previously discussed and verified models by Schmidt et al. [26] and Sun et al. [27,28].

A. Drag law
Particles are assumed to be spherical (Clift et al. [33]) and reasonably modeled in a Lagrangian way following Newton's second law of motion [34]. The main assumption for the point-source approximation is that the particle size is small relative to the Kolmogorov length scale of the fluid phase. For low particle Reynolds number and high density ratios between the fluid and the particle phase, the governing equations for each individual particle, here omitting gravity and adopting the point-particle approximation, are: The subscripts p and f represent the particle and the fluid phase, respectively, and the convention is adopted that the coordinate i = 2 is aligned with the ODT coordinate direction y. Here, v f ,i represents the undisturbed fluid velocity. x p,2 is the particle spatial coordinate on the ODT line and v is used instead of u to denote velocities because the velocities used to advance these equations are not necessarily the same as the ODT fluid-phase velocities u i . The velocities v are accordingly termed eddy velocities. Where the component index i is omitted, the component i = 2 aligned with 044308-8 the ODT domain is implied, hence y p is shorthand for x p,2 . ODT does not implement particle displacements in directions i = 2 because ODT captures only the domain-aligned displacements, which is sufficient for modeling representative behavior in isotropic turbulence. Nevertheless, there is an intermediate step in the advancement process that requires evaluation of particle displacements ) based on Stokes flow, is expressed here in terms of the diameter d p , the particle-to-fluid density ratio ρ p /ρ f , and the fluid kinematic viscosity ν. (Where convenient, the dynamic viscosity μ = ρ f ν is introduced.) An empirical correction factor f is specified as based on studies of Schiller and Naumann [35] for nonslip Reynolds numbers Re p smaller than 800 to capture the standard drag curve [33]. This correction is used in the DNS comparison cases, so its use here maintains consistency of the ODT formulation with the DNS cases. Here, Re p is evaluated with the interacting gas velocities, which differ for particle time advancement and particle-eddy interaction, as explained in what follows. The drag law [Eq. (24)] is solved by a first-order Euler method. Re p statistics that are presented in Sec. IV C for representative cases verify that Re p is well within the range of validity of the correction.

B. Particle time advancement
Much as the fluid-phase treatment consists of time advancement of viscous transport punctuated by eddy events, the particle treatment consists of time advancement of the drag law, Eq. (24), punctuated by interaction of particles with the fluid phase during the eddy events, which involves a different application of the drag law, as described in Sec. III C. Particle advancement concurrent with viscous transport is described first.
One might suppose that the velocity components v f ,i in Eq. (24) can be set equal to the ODT fluid velocities u i , but this is not the case because the advection of fluid elements occurs only by triplet mapping. One implication for particles is that, in the zero-inertia limit of vanishing τ p , particles cannot deviate from the trajectory of the fluid elements that contain them. During viscous transport, fluid elements are motionless in the ODT domain direction. Then zero-inertia particles must likewise be motionless. The only value of v f ,2 in Eq. (24) that is consistent with this requirement is zero. Considering the case of vanishing τ p in Eq. (24), v p,2 is forced to the value of v f ,2 and hence to zero. These results imply that there is no i = 2 contribution to the computation of Re p for the time advancement. In the other coordinate directions i = 2, v f ,i has no bearing on the displacement of particles relative to fluid parcels, so v f ,i is taken to be u i for those components, where u i is subject to viscous-transport advancement concurrently with particle time advancement.
Another limiting case that must be treated consistently is the infinite-inertia ballistic limit. Equation (24) gives constant-velocity motion in this limit for any time history of the fluid velocity components, so the treatment of this limiting case is trivially correct. In contrast, the formulation of the particle-eddy interaction requires some care in order to enforce the correct behavior in this limit, as described in Sec. III C.
As particle momentum changes during the integration of Eq. (24), overall momentum conservation is maintained by distributing an equal and opposite momentum change uniformly to the fluid within the mesh cell that contains it, thus implementing one of the particle-fluid two-way coupling mechanisms in the model. (The other is described in Sec. III D.) This feedback to the fluid phase during particle time advancement is indicated by the particle contribution S p,i in Eq. (2), which is evaluated as  where N c is the number of particles in the mesh cell (but particle and mesh cell indices are omitted in the summand for clarity). This treatment of two-way coupling is standard in DNS with point particles, and in that context is known to dissipate total kinetic energy, which is unphysical because only the fluid phase can dissipate kinetic energy [20]. This is intended as a compensation due to the fact that not all true fluid dissipation is resolved on the grid. In any case, the key behavior of this formulation with regard to its role in the outcomes described in Sec. IV C is the net loss of particle kinetic energy. This is largely due to the condition v f ,2 = 0, which assures that drag-induced momentum change reduces the i = 2 component of particle kinetic energy. Because mean particle kinetic energy, like any global quantity, must be constant in the quasisteady state, there must be a commensurate source of particle kinetic energy. This source is the particle-eddy interaction that is described next.

C. Particle-eddy interaction
Due to the instantaneous character of eddy events it is necessary to model the interaction between a particle and an eddy event, which we refer in the following to as particle-eddy interaction (PEI). Here, it is important to have in mind the distinction between physical vortex/eddy structures and eddy events in ODT. Schmidt et al. [26] developed the so-called instantaneous particle-eddy model (noted as type-I), which governs the domain-aligned displacement of particles due to an eddy event. This PEI model is used to capture particle-eddy interaction for each particle that is located in the sampled eddy region. The main model assumption is that the eddy lifetime equals the eddy time scale τ e [Eq. (15)] times a coefficient β p that is a model parameter. The PEI can persist no longer than the eddy lifetime but can be shorter in duration if the particle exits the volume occupied by the eddy, defined below, before the PEI terminates at the end of the eddy lifetime. That means the integration of the motion of a particle in the ODT line direction based on Eq. (24) must take into account the time until a particle exits the region occupied by the eddy if that time is less than the eddy lifetime.
Before the latter point is addressed, the treatment of the drag law during the PEI is explained. It will be shown that modifications of the drag law are needed in order to satisfy consistency conditions. It is useful for this purpose to simplify the drag law by specifying that v f ,i and τ p / f are constant in time during the PEI. For i = 2, the constant value is chosen to be the value of u i at the particle initial location y p0 . The specification of v f ,2 is explained shortly. The auxiliary variables x p,i for i = 2 are initialized to nominal values x p0,i = 0.
On this basis, the analytical result of integrating the drag law over an arbitrary time interval t is yielding the particle trajectory In Sec. III B, v f ,2 was set equal to zero in order to enforce consistency between particle and gas motion in the zero-inertia limit. The PEI likewise requires modification, in this instance to enforce consistency in the infinite-inertia (ballistic) limit. The reason is that the PEI advances particles during a time interval t pei starting at the instant of an eddy occurrence, but the resulting modified values of particle locations and velocities are applied at that instant rather than at a time t pei after eddy occurrence. Upon completion of the PEI, the particle time advancement described in Sec. III B resumes, resulting in a second implementation of the advancement during the ensuing time interval t pei , albeit using a different algorithm. This double counting of particle advancement has no consequence for zero-inertia particles because such particles are motionless during the second advancement interval as a result of the consistency enforcement described in Sec. III B. All finite-inertia particles are subject to this artifact, but there is an exact remedy only for ballistic particles. As noted in Sec. III B, the advancement 044308-10 Eddy velocity v f ,2 is defined as the massless particle displacement by the triplet map divided by the particle-eddy interaction time t pei . For the displacement one (black circle) of three possible positions (gray circles) is chosen randomly.
after eddy occurrence correctly produces the trajectory of a ballistic particle. So the artifact can be avoided by modifying the drag law solution, Eqs. (27) and (28), so that its specialization to the ballistic limit produces no change in y p or v p,2 when the drag law is applied during the PEI. This gives the differences in position and velocity The corresponding post-PEI particle location and velocity are y p = y p0 + y p and v p,2 = v p0,2 + v p,2 , respectively, where v p0,2 appears in the latter equation because it has been zeroed out only in the expressions for the changes.
In the ballistic limit τ p → ∞, Eq. (29) yields no change in y p or v p,2 , as desired. Equation (29) is likewise consistent with the zero-inertia limit, as explained below after the quantities v f ,2 and t pei are defined.
For finite St, the modifications that yield Eq. (29) and the resulting expressions for y p and v p,2 do not constitute a physically consistent drag-law formulation. Nevertheless, this formulation smoothly interpolates between limits in which the formulation is physically correct. From this viewpoint the approximation is no more severe than idealizations inherent in the overall modeling framework. Now it is necessary to define the eddy velocity v f ,2 and to evaluate the interaction time t pei . Determination of v f ,2 is based on the displacement of a massless (zero-inertia) particle by a triplet map as specified by Eq. (3). The triplet map provides three possible massless particle positions and a unique position is sampled randomly with a uniform distribution from those three possible ones. (The mathematical justification for this is presented in Appendix A.) The chosen displacement, see Fig. 2, divided by t pei defines the fluid velocity v f ,2 during the PEI. Based on the zero-inertia limit of Eq. (29), this enforces the zero-inertia consistency condition.
As a next step the interaction time scale t pei has to be determined and therefore a so-called eddy box is introduced with the dimensions [l × l × l], corresponding to the eddy interval l in ODT line direction and intervals [−l/2, l/2] in each of the two coordinate directions normal to the ODT line direction. If the particle exits the eddy box through any face during the eddy lifetime, then the duration of the PEI is deemed to be the time until the exit of the box. Otherwise it is the eddy lifetime.
The eddy box reflects the consideration that the PEI can terminate either because the eddy motion has ended or because the particle leaves the volume containing the eddy before the end of the eddy motion. Although the model lives on the 1D domain, it evolves all three components 044308-11 of fluid and particle velocity and therefore can evolve the particle trajectory in three dimensions under simplifying assumptions such as constancy of domain-normal fluid velocity components in the domain-normal directions. The eddy box then defines the boundaries of the eddy volume, allowing the particle exit time to be evaluated accordingly.
For directions i = 2, the time until exit (which need not exist because the particle could come to rest within the eddy box) is based on the direction-i trajectory given by Eq. (28). For i = 2 it is given by the expression given below Eq. (29), where the time argument in Eq. (28) is now the generic time variable t. The first crossing of any face of the eddy box determines the exit time.
v f ,2 is defined in terms of t pei , but now there is the possibility that t pei depends on v f ,2 because the exit time might be determined by the particle exiting the eddy box in the y direction. This possibility would exist if the drag law were implemented without modification. However, the initial particle velocity in direction i = 2 has been set equal to zero to obtain Eq. (29), which assures that the new particle position is in the range between initial particle position and the massless particle displacement for St > 0 and monotonically approaches the initial particle position as St increases (by construction, for consistency with the ballistic limit). Therefore particles cannot exit the eddy box in the y direction, so only directions i = 2 need to be checked. This is why it is possible to evaluate t pei before evaluating v f ,2 .

D. Interphase coupling terms
Eddy events displace particles and change their momentum in the domain-aligned direction. As in the treatment of particle time advancement in Sec. III B, this implies momentum exchange and thus energy exchange between the fluid and particles within the eddy volume, implemented in this case as the final step of the eddy event. S i and S E ,i in Eqs. (5) and (13), respectively, enforce fluid momentum and energy changes that are equal and opposite to those of the particle phase, where and These are summations over the N e particles within the eddy volume, with the summation indices suppressed in the summands.

IV. RESULTS
As a summary, using ODT to study particle-laden, stationary forced HIT requires the following steps to set it up: 1. Set single phase turbulence parameters (C, Z); 2. Set forcing parameters for TKE injection (P, T 11 , l); 3. Set particle-eddy interaction parameter (β p ), where Steps 1 and 2 tune ODT to match the single-phase DNS and only Step 3 involves tuning to match particle-laden DNS data. Steps 1 and 3 are validated, here, through comparison with DNS, but once found they are independent of the turbulence level and can be used globally without new validation. In Step 2 only the length parameter l of the kernel event has to be validated, as the kernel event geometry is not directly reproducing an equal size turbulent length scale. However, the production rate P and time scale T 11 can be directly applied as in DNS.

A. Validation of unladen case
The ODT data is validated against results from direct numerical simulations (DNS) studied by Abdelsamie and Lee [6] with Taylor microscale Reynolds number Re λ = 2k √ 5/(3ν ) of 70 for the 044308-12  Table I including macro and micro time and length scales. The microscales τ η and η and the Taylor microscale λ are based on combinations of the kinematic viscosity ν, the TKE k, and its dissipation rate as follows: The value for k in ODT is implicitly set through the choice of eddy frequency parameter C [Eq. (15)] to match the DNS k value and thus the inferred quantities Re λ and λ. The integral length and time scales, L 11 and T 11 , are given as where E (κ ) is the kinetic-energy spectrum. u is the average fluctuation given as u = √ 2 3 k. L 11 and T 11 are predicted in contrast to the outputs that precede them in Table I.
The ODT domain is set up as follows. Periodic boundary conditions are used. To minimize artifacts of periodicity in simulations of homogeneous turbulence, such as truncation of the range of eddy sizes, the ODT domain length D is set to D = 14, which is on the order of ten times the integral length scale. The largest occurring eddy event is limited to 90% of the domain size. This is an input parameter of the ODT simulation and it corresponds statistically to a rare large eddy event. It has turned out to be necessary in ODT to capture the frequency with the most energy in the nondimensionalized kinetic-energy spectrum.
The initial velocity profile for all test cases is a fully developed turbulent profile, which was generated by a single-phase ODT simulation until it reached stationary behavior for TKE. ODT is solved on an adaptive mesh, but to gather the statistics of production (kernel events) and rate of dissipation of TKE during the simulation, the current system state is duplicated, with interpolation, on a uniform grid with 500 grid cells that is used for further data reduction. Additionally, the simulation creates output files of the current state at uniform time intervals. These are used to obtain the kinetic-energy spectrum. Based on a parameter study the ODT model parameters are set to C = 6.0 and Z = 600.

B. Validation of particle-laden case
The particle phase is completely defined for a two-way coupling simulation by the three quantities, Stokes number St, mass loading φ, and density ratio ρ p ρ f . For one-way coupling between particle and fluid phase the two latter quantities are not required. Using these three quantities and fluid properties in Table I, the particle relaxation time τ p , the particle diameter d p , the mass of a TABLE II. ODT particle characteristics grouped as inputs that define the cases and quantities inferred from the particle and fluid (Table I) inputs. The table identifies quantities needed for one-way coupling and others that are additionally needed or convenient for simulating two-way coupling. single particle m p , and the particle number density n p can be calculated by, Here, ρ f is assigned the nominal value 1.0 kg/m 3 for the purpose of evaluating μ and each particle is assumed to be spherical to calculate their volume. After particles are introduced into an equilibrated single-phase ODT state, the flow is time advanced for sufficient additional time for the two-phase flow to reach statistical stationarity before the gathering of output statistics begins. The case-specific values for the comparisons of ODT to DNS results of Abdelsamie and Lee [6] are shown in Table II. Consistent with the DNS results, St is not reevaluated for case-specific values of τ η . Instead, the single-phase value is used in order to simplify the interpretation in the following sections. All cases are monodisperse because the available DNS comparison cases are monodisperse. However, the model straightforwardly generalizes to polydispersions. All presented cases result in particle volume loadings, which are in the two-way coupling regime without significant collision effects (which are therefore neglected here) [36]. As mentioned before, ODT is solved on a nonuniform adaptive mesh with a minimum size of 0.001 of the domain length, which gives d p (St = 1)/ min = 0.25 and d p (St = 5)/ min = 0.58. However, the simulation almost never reached that point.
The PEI model parameter β p is set to 0.8 based on a parameter study, which can be found in Appendix B. That study shows that the sensitivity of the kinetic-energy spectrum to the value of β p is negligible.
Equation (25) is valid provided that the nonslip velocity Reynolds number Re p is smaller than 800. Table III displays the average Re p for both St, showing that the values for both cases are orders of magnitude lower than the limit. The results show higher Re p for higher St, reflecting the increase of particle slip with increasing inertia.

C. ODT vs DNS
ODT and DNS have significant differences regarding the forcing scheme that are important to have in mind when comparing the results. To force the stationary turbulence, the DNS uses source terms in the large-scale modes of the spectral momentum equation [19] that interact directly with the source term caused by the fluid-particle interaction [1], leading to a modified rate of TKE production relative to an unladen case with the same source terms. Therefore , and hence quantities inferred from , are DNS outputs rather than inputs owing to coupling between the turbulence forcing mechanism and the fluid-particle interaction.
In contrast, the ODT forcing scheme acts independently of the viscous advancement such that the TKE production rate resulting from forcing by kernel events is unaffected by particle loading. However, by adjusting the production of TKE for each test case to match the DNS value of , ODT is capable of capturing the DNS results as seen in Figs. 3 and 4. Table IV shows the effect of particles on k and , where the latter was externally enforced in ODT. However, because this is a driven, dissipative system, k predicted by ODT need not match the DNS k values, where the difference reflects model error, which are explained below. The results show that ODT qualitatively captures the decay of k with some underprediction towards St = 5. The reason for it can be seen in Fig. 3, which displays the nondimensionalized kinetic-energy spectrum. Due to the different forcing schemes in ODT and the DNS the particles are interacting with different turbulence scales. The spectrum for St = 5 in ODT shows roughly uniform multiplicative reduction  of spectral amplitude over all wavelengths and, unlike the DNS, does not cross the single-phase spectrum above the Kolmogorov wave number. The spectrum for St = 1 follows a similar trend, but with less amplitude reduction than for St = 5. However, the different adjusted dissipation rates can contribute to the difference between the St = 1 and St = 5 spectra. Figure 4 shows this behavior by displaying the ratio of TKE with and without particles. The most significant difference is in the high-wave-number region, in which the DNS spectrum is strongly enhanced but the ODT spectrum is not. Nevertheless it is seen that ODT captures the overall tendency of the DNS results. Because  [7] found in their study with the parameters Re λ = 58, St = 9.2, and φ = 0.11 a value for k/k φ=0 = 1.06 and so an increase of TKE.
To the authors' knowledge the ratio of particle and fluid phase mass averaged kinetic energies has not been reported in this or other DNS studies. Figure 6 provides it for the same parameter range. For St = 5, k p /k values are relatively constant with variations between 0.25 and 0.3, whereas for St = 1, k p /k has a larger spread with variations between 0.35 and 0.5 for corresponding φ values from 0.2 to 1. To remove the complicating influence of case-specific adjustment of the dissipation rate, the same test cases were rerun using the dissipation rate for the single-phase case, as described in Sec. IV D.

D. ODT with constant dissipation
In the following the results of particle-laden ODT simulations with the same dissipation rate as the single phase, given in Table I,  for each Stokes number are evaluated. Figures 7(a) and 7(b) show the nondimensionalized kineticenergy spectra with different mass loadings for St = 1 and St = 5, respectively. The dependence on Stokes number is not as significant as it was with the adjusted dissipation. The differences are highlighted by Figs. 8(a) and 8(b), which present again the quotients of particle-laden and singlephase test cases.
For φ = 0.4, the effects of St = 1 and St = 5 particles on the kinetic-energy spectrum are quite similar. Both decrease the TKE over the whole spectrum and meet the single-phase spectrum at the largest wave number, kη > 1. The S-shape form suggests strongest attenuation in the wave number range 0.2 < kη < 0.8, with stronger attenuation for St = 1. Higher mass loadings are seen to accentuate these behaviors for both St values. On each plot, the curves for all mass loadings cross near kη = 0.8.
The significant difference between St = 1 and St = 5 is that for the latter case, the quotient curves appear to level off close to the single-phase spectrum value at kη > 1, but for St = 1 the curves exceed that value at the highest wave numbers. These tendencies are consistent with the fact that the interphase coupling influence is weakest at the critical point 1/η for St near unity [3], and they also reflect the increase of particle influence on the fluid phase as mass loading increases.
Holding the dissipation fixed is thus seen to remove some of the apparent case sensitivity seen when the dissipation varies. As noted, loading dependence of the DNS dissipation is a consequence of coupling of the forcing and the particle-flow interaction, so it is appropriate to base comparisons on constant dissipation. A related consideration is that, even with constant dissipation, St varies with loading because TKE decreases with increasing φ, as shown below in Fig. 9, implying an increase of τ η . In principle this could be compensated by evaluating τ η on a case-specific basis so as to obtain case-specific St values. It would be interesting to ascertain how much of the apparent φ dependence would thus be absorbed in case parametrization based on a more relevant evaluation of St. This level of analysis would rely for validation on reanalysis of DNS as well as ODT outputs. For this and other reasons, it is beyond the present scope and therefore is deferred for future investigation. For now, this consideration is recognized as a caveat to the further interpretations that are proposed below.
Notwithstanding the noted complications, it is instructive to examine the probability density function (PDF) of the time scale τ e of implemented eddy events, shown in Fig. 10. τ e depends on the eddy length and its available kinetic energy [see Eq. (15)]. The overall effect of increasing St or particle loading is to shift the PDF to the right, that is, towards eddies with larger time scales. This is reflected by the reduced area under the PDF curve at small τ e /τ η , which by PDF normalization raises the curve elsewhere, albeit by a small apparent amount in the log-log coordinates. The eddy time scale increases with decreasing available energy. Increase of particle energy during the PEI results in such a decrease of the eddy available energy. This indicates that particles are net sinks of the available energy of the frequent small eddies. Owing to statistical stationarity, there must be an equal and opposite flow of kinetic energy back to the fluid. As noted in Sec. III B, particle time advancement between eddy occurrences has this tendency. Moreover, the point-particle two-way-coupling algorithm in the viscous advancement [Eq. (2)] in ODT transfers the energy to the smallest available fluid scale, namely the mesh cell containing the particle.
Overall, this picture is consistent with the shapes of the DNS and ODT curves in Fig. 4. Namely, the dips reflect the depletion of fluid TKE implied by Fig. 10 and the upswings at the highest wave numbers reflect the redeposition of that energy into the fluid at the smallest scales. Because the latter is algorithmically controlled in both DNS and ODT, the difference between their outputs at the highest wave numbers cannot be interpreted unambiguously in the absence of a comparative study of the associated numerical artifacts of the two simulation methods.
It is additionally feasible to interpret details of the parameter dependences seen in Fig. 10. there is no particle-fluid interaction expected for St = ∞. Figure 10(a) is consistent with maximum particle influence for St near unity. Variation of mass loading might scale the frequency and strength of particle-turbulence interactions that act as both sources and sinks of TKE without affecting the net distribution of τ e values, consistent with Fig. 10(b). However, there are several associated considerations such as possible effects of mass loading on the effective values of Re λ and St, so further investigation would be needed in order to arrive at a definitive interpretation of Fig. 10(b).
The trend seen in Fig. 9 up to φ = 1.0 resembles the trend in Fig. 5, but for constant dissipation there is a more gradual decline of k with increasing φ. Figure 9 shows additional results for φ = 2.0 that indicate a reversal of the St dependence. This could reflect a variety of possible influences, ranging from the case parametrization issue that has been mentioned to increasing importance of particle clustering at high mass loading. (Triplet mapping of inertial particles is known to induce clustering that is consistent with physical and mathematical requirements [37,38]) In view of the limited evidence and the potentially complicated phenomenology at high mass loading, no interpretation of the φ = 2.0 results is attempted, but they serve as additional predictions that are subject to possible future validation by DNS or other means. Figure 11 shows k p /k, similarly to Fig. 6. Here, Fig. 11 illustrates that apparent complexity can partly reflect the choice of case parametrization. For constant dissipation, there is an overall 044308-20 downward trend with increasing φ, in contrast to the complicated dependence seen in Fig. 6. This follows the trend of the fluid phase TKE, which is the source of the particle phase TKE, so the latter varies in tandem with the former.

V. DISCUSSION
The ODT strategy for economical turbulence simulation using a representation of flow evolution on a 1D domain was previously extended to particle-laden flow and then to two-way fluid-particle coupling. Here, a version of the latter combined with a novel approach to large-scale forcing has enabled an ODT study of turbulence modulation by particles in stationary HIT, highlighted by validation against several sets of DNS results for this flow.
The fluid-particle coupling in ODT requires significant modifications of the drag law treatment relative to the drag law implementation in DNS. The guiding principles that preserve the key physics of this coupling are that particle motion in the zero-inertia limit must exactly conform to the local fluid motion and that infinite-inertia particles must travel ballistically, i.e., at constant velocity. The formulation smoothly interpolates between these limiting behaviors, yielding qualitatively correct phenomenology and a useful degree of quantitative accuracy.
In particular, ODT reproduces the spectral profile of turbulence modulation that is induced by the two-way coupling in DNS. An important feature of this modulation is the shift of spectral intensity from scales somewhat larger than the Kolmogorov microscale to scales in the vicinity of the Kolmogorov microscale. Owing to the transparent phenomenology implied by the model formulation, this behavior has been readily diagnosed. The dissipative nature of point-particle time advancement in DNS as well as ODT contributes to the high-wave-number upturn. This raises the question of the extent to which this behavior is physical rather than a numerical artifact of both methods. For ODT, there is a nondissipative alternative that has not been described here. Its formulation might raise other concerns about physical fidelity, but it might at least be informative in terms of sensitivities and therefore will be addressed in future work. It is also noted that an extension based on Horwitz and Mani [20] is planned to improve the computation of the undisturbed fluid velocity in the drag law [Eq. (24)].
These points illustrate that model features enable physics investigation presently beyond the scope of DNS or any other existing methodology. The ODT turbulence forcing mechanism is nondissipative for particle-laden as well as single-phase flow, so parameter studies are possible with the dissipation rate held fixed while varying any combination of Stokes number, particle mass loading, and fluid-particle density ratio. It has been shown that this can yield more consistent trends than DNS for which the dissipation rate is affected by details of the fluid-particle coupling. Additionally, the run time of the simulation is much lower and needs for example for the test case Re λ = 70, St = 5, and φ = 2.0 only 9 h on a single CPU (3.1GHz). The low computational cost of ODT creates opportunities for extensive future parameter studies that can potentially explore, e.g., alternatives to defining Stokes number based on the unladen flow, such as inferring the case-specific turbulence microscale from the flow statistics in order to obtain a physically more meaningful Stokes number. Irrespective of the accuracy of point predictions, this might identify an approach to rationalizing the parameter dependences of the two-way coupling regime of particle-laden turbulence.
The ratio of particle-phase and fluid-phase TKE is a fundamental property of particle-laden turbulence whose evaluation based on DNS of stationary HIT does not appear to have been reported in the literature. ODT results for this ratio are reported for various combinations of St and φ, and can readily be generated over a wider parameter space encompassing density ratio and Re λ , including values of the latter that are unattainable using DNS. Validation of these predictions by DNS or other means would be of interest.
Although stationary HIT is a useful test case in terms of relative simplicity of interpretation and availability of validation data, it is not the most relevant case for generating information useful for SGS closure of LES. SGS turbulence in LES is governed mainly by the mesh-resolved shear, hence the velocity differences across the coarse-grained control volumes. The natural paradigm for this configuration is homogeneous sheared turbulence (HST). The periodicity of an HST simulation then corresponds to the control-volume size. Accordingly, particle-laden periodic HST will be another focus of two-way-coupled particle-laden ODT research in the future.

ACKNOWLEDGMENT
This project has received funding from the European Union Horizon-2020 Research and Innovation Program. Grant Agreement No. 675676.

APPENDIX A: MATHEMATICAL BASIS OF THE RANDOM CHOICE AMONG MASSLESS-PARTICLE DISPLACEMENTS BY THE TRIPLET MAP
As illustrated in Fig. 2, triplet mapping of a massless particle identifies three possible particle displacements. A unique displacement is obtained by randomly choosing one of the displacements, all of which are equally probable. This procedure was previously presented as a model assumption with no statement of an underlying physical or mathematical basis [37]. Here, the mathematical basis of this procedure is explained.
Application of the triplet map to a function of y corresponds operationally to a three-step transformation. First, the portion of the function that is within the mapped interval is compressed by a factor of three. Second, the original function h(y) is replaced within the interval by three copies of the compressed function so as to fill the interval without gaps or overlaps. Third, the middle copy is flipped (reflected with respect to its midpoint). The resulting mapped function h (y) obeys the model analog of the solenoidal condition, namely y : a<h(y)<b dy = y : a<h (y)<b dy for any a and b > a. Choosing h(y) to be any two-valued (0 or 1) indicator function clarifies that this relation is equivalent to the requirement that the measure of any subregion of the y domain before a map is the same as the measure of the subregion onto which it is mapped. This is the 1D analog of volume preservation in 3D solenoidal flow, as illustrated from another perspective in Fig. 12. However, this conservation property becomes ill defined for a generalized function such as h(y) = i δ(y − y i ), i.e., a sum of δ functions that might in the present context represent a set of particle locations y i . The underlying issue is that the effect of the map on the dependent variable h is uniquely defined but its effect on the independent variable y is nonunique in that it specifies three distinct results y of the map y → y . An alternative map representation that is commonly used in numerical implementations [39] resolves this ambiguity. A discrete triplet map is defined by partitioning the mapped interval into 3k equal-sized segments, where k > 1. The discrete map applied to these segments, sequentially labeled 1, 2,. . ., 3k is a permutation that yields the segment sequence 1, 4, 7,. . ., 3k − 8, 3k − 5, 3k − 2, 3k − 1, 3k − 4, 3k − 7, . . ., 8, 5, 2, 3, 6, 9,. . ., 3k − 6, 3k − 3, 3k. For example, an interval consisting of six segments sequentially labeled 1, 2, 3, 4, 5, 6 is rearranged to obtain 1, 4, 5, 2, 3, 6. The k = 3 discretization is shown in Fig. 12.
Indexing these maps by k, the resulting transformation h(y) → h k (y) of the discretized dependent variable h converges with increasing k to the continuum map h(y) → h (y) for suitably regular functions h(y). On this basis, the continuum triplet map is redefined as the large-k limit of the sequence of discrete maps. Consider the subset of discrete maps specified by k = 2 × 3 j−1 for j 1. The maps, now indexed by j, can be viewed as a sequence of successive spatial refinements, where a refinement partitions each segment into three identical smaller ones. Denote the y range of the mapped interval as [0,2], which has the convenient property that 3 j y ranges from 0 to 3k. Let I j (y) denote the integer part of 3 j y. Then for y in segment n, I j (y) = n − 1. Label the three compressed copies produced by a map by the indices 0, 1, and 2 going from left to right. Then by inspection of the permutation rule, for every j and every positive n 3k, segment n is mapped to copy (n − 1) mod 3, and thus to I j (y) mod 3. Next, y and multiples of y are expressed in base 3, so I j (y) is the integer part of 10 j y. The factor 10 j promotes j numerals of the fractional part of y to the integer part. Then evaluation of the integer part modulo 3 yields the jth numeral in the fractional part of y.
Combining results, that numeral selects the mapped copy to which y is displaced, so the ternary representation of the real number y specifies the sequence of copy choices for any map refinement j. It is now asserted that y is almost surely simply normal in base 3, meaning that in the limit J → ∞, the first J numerals in the ternary representation of y have the value 0, 1, or 2 with equal frequency of each. This holds because every simply normal number in any base is a normal number, and according to the normal number theorem [40], almost all real numbers are normal. Then for sufficiently large j, the likelihoods of the three copy choices are equal to within an arbitrarily small tolerance, suggesting that in the continuum limit j → ∞ the mathematically correct procedure is to randomly sample the copy to which location y is displaced. The discrete maps indexed by j are a subset of those indexed by k. The remaining values of k are treated as follows. Identify the smallest k that does not obey k = 2 × 3 j−1 for some positive integer j. This value k * is the first member of a subset k = k * × 3 h−1 for h 1 to which the same reasoning is applied, with the y range of the mapped interval specified for this purpose as [0, k * ]. Then the smallest value k * * not included in either subset similarly becomes the first member of another such subset. By induction, every k value is eventually included in at least one subset. Thus it is concluded that to within any specified tolerance, however small, there is a k value above which it is assured that the likelihoods of the three copy choices are equal, thereby confirming the correctness of the random sampling procedure used in the continuum formulation. This required augmentation of the continuum map specification confirms that δ-function singularities of h(y) are not sufficiently regular for the convergence property of the sequence of discrete maps to be satisfied.

APPENDIX B: PARAMETER STUDY FOR PEI PARAMETER β p
The PEI model parameter β p is adjusted to obtain agreement with DNS results for particle dispersion. In Fig. 13, mean-square particle displacement σ 2 is compared for a range of β p values to the DNS results of Abdelsamie and Lee [41]. The reported DNS data was limited to a short time interval during which the ODT results look linear rather than quadratic as indicated by the DNS. However, over a longer time interval ODT shows a quadratic behavior, too. It can be seen that β p = 0.8 agrees best with the DNS results and therefore is used in the present study. This value differs from the value used by Schmidt et al. [26] in their study because the present PEI treatment differs in some respects from their formulation. Figure 14 shows that the kinetic-energy spectrum quotient for the Table IV case settings has negligible sensitivity to β p . This implies that the tuning of β p is not an important contribution to the degree of agreement between ODT and DNS results that has been obtained.