An approximate inertial manifold (AIM) based closure for turbulent flows

A closure model for turbulent flows is developed based on a dynamical system theory. An appropriately discretized formulation of the governing equations is considered for this process. The key ingredient is an approximation of the system’s attractor, where all the trajectories in phase space are confined. This approximate inertial manifold based approach provides a path to track trajectories of the system in a lower-dimensional subspace. Unlike conventional coarse-graining approaches, the turbulent field is decomposed into resolved and unresolved dynamics using the properties of the governing equations. The novelty of the approach relies on the reconstruction of the unresolved field constrained by the governing equations. A posteriori tests for homogeneous isotropic turbulence and the Kuramoto–Sivashinsky equation show promising results for considerable dimension reduction with strong convergence properties. The proposed model outperforms the dynamic Smagorinsky model, and the computational overhead is competitive with similar approaches


I. INTRODUCTION
2][3] Although LES resolves only the large scales, and small scales have to be explicitly modeled, its ability to capture turbulent mixing has led to the increased predictive accuracy of reacting flows. 4However, the accuracy of LES is still dependent on small-scale modeling, especially when the dominant physics occurs at these scales. 3,5,6The effect of small scales is treated using two general approaches: explicit and implicit formulations.In implicit LES, the physics of the small scales is represented by the numerical dissipation imposed by the spatial discretization of the governing equations, assuming that the contribution of the small scales is strictly dissipative. 7,8In explicit modeling of the small scales contribution, spatial filters are applied to the governing equations to remove scales smaller than a prescribed value, introducing additional subgrid-scale (SGS) stresses and fluxes to be modeled.In this latter approach, the SGS contribution is modeled either indirectly as an algebraic function of the filtered field properties, such as eddy-viscosity models, or directly by the reconstruction of the sub-filter field. 9,10The focus here is on the explicitly modeled LES formulation.
Eddy-viscosity models represent SGS stress tensor and fluxes in terms of the filtered field variables.In particular, the Smagorinsky model 11 assumes that the SGS stress tensor is proportional to the filtered strain rate, and the proportionality factor, eddy viscosity, depends on the filter width and a constant model parameter.Despite its simplicity in implementation, in complex flows where multiple flow regimes exist, the proportionality coefficient needs to be tuned for each flow regime. 123][14] In this modeling ansatz, the local equilibrium of turbulent kinetic energy plays a key role, which implies a forward cascade of energy between resolved and modeled scales.These models have been found to introduce errors in near-wall and transitional regions. 14,15 different approach is to approximate the unfiltered field from the information at the resolved scales without invoking the universality of the subgrid field.These models are developed mostly on the scale similarity premise, i.e., the contribution of the small resolved scales to the large resolved scales is similar to the contribution of the unresolved scales to the small resolved scales.Scale similarity model (SSM) 16,17 approximates the SGS stress tensor by applying a second spatial filter like in the dynamic Smagorinsky model.There, the proportionality constant is determined to reproduce the exact average SGS kinetic energy.SSM gives the correct rate of energy flux to the subgrid scales and can predict backscatter reasonably.A generalized form of the scale similarity model with repeated filtering is the approximate deconvolution model (ADM). 9,10In ADM, the unfiltered field is recovered by applying the truncated series expansion of the inverse filter operator to the filtered field; hence, nonlinear terms in the governing equation of the resolved scales are computed directly.Unlike eddy viscosity models, SSM and ADM cannot predict adequate SGS dissipation, and an additional relaxation regularization is required.Bardina et al. proposed a mixed-model of scale-similarity and eddy-viscosity models to account for the twofold SGS contribution: (1) energy transfer from large scales and (2) dissipation of energy contained in the SGS.This model has a superior performance in transitional flows. 17However, the ratio of each component of the SGS model needs to be determined.
There exist other LES modeling approaches that do not rely on traditional spatial filters for the separation of scales.Such models use projection-based decomposition for scale separation.For instance, variational multi-scale LES 18,19 uses a variational projection to decompose the range of scales in groups of large resolved scales, small resolved scales, and unresolved scales.In this approach, the direct contribution of the SGS physics is confined to the small resolved scales, and large resolved scales are solved directly (i.e., without any modeling) but influenced indirectly by the subgrid-scale model due to the inherent coupling of all scales. 18,20Projectionbased scale separation provides a theoretical framework for model reduction that is largely used in dynamical systems.For instance, this approach has been successful in weather prediction. 21A fluid system described by a set of partial differential equations (PDEs) governing nv variables can be treated as a finite-dimensional dynamical system after appropriate spatial discretization using ng grid points.In this setting, the dynamics of the system reside in an N-dimensional state-space defined by nv × ng degrees of freedom.The spatial and temporal evolution of the turbulent flow can then be expressed as a trajectory in this state-space.In many systems dominated by coherent structures, the long-time behavior of the system is known to be confined to a low-dimensional subspace of the full N-dimensional state-space.3][24] Constantin et al. 23 showed that the dimension of the attractor scales nonlinearly with the Reynolds number of the flow.Direct estimations of this attractor dimension for turbulent flows using the Kaplan-Yorke conjecture 25 showed that attractor dimensions are orders of magnitude lower than the number of degrees of freedom required by DNS. 26,27n certain dissipative dynamical systems, an invariant subset of the state-space attracts all trajectories of the system exponentially.The long-time dynamical behavior of the system can be studied in this low-dimensional subspace, known as the inertial manifold (IM). 24These manifolds, when they exist, contain the global attractor of the system.The dynamics of the inertial manifold can be described by a finite-dimensional system of ordinary differential equations (ODE), called the inertial form, which completely describes the long-time dynamical behavior of the original infinitedimensional system.9][30] However, these theorems cannot explicitly determine the topology of the inertial manifold, and they can only provide an upper bound for its dimension.Inertial manifold theories require strong restrictions, which are not satisfied in practical systems governed by the Navier-Stokes equations. 24,28Consequently, an approximation of the inertial manifold (AIM) is necessary to describe the system in its inertial form.4][35][36][37] The main assumption is that the dynamics in the complement space, between the AIM and the full state-space, equilibrate to the dynamics on the AIM.In other words, the motion in the complement space responds instantaneously to changes in the trajectory on the AIM.While this assumption is justified by the theoretical studies discussed above, its validity needs to be scrutinized more rigorously for different systems.
Recently, the use of AIM for turbulence simulations has been pursued. 37,38In one study, 37 a priori analyses of the AIM formulation for the one-dimensional Kuramato-Sivashinsky equation (KSE) and the three-dimensional Navier-Stokes equations were carried out.It was demonstrated that the AIM approach provides a viable pathway for modeling the unresolved scales.In a related study, 38 the AIM approach was used to model non-premixed turbulent combustion.In this particular case, the chemical reactions are known to occupy a lower-dimensional manifold.The proposed AIM correctly identified this manifold, without being prescribed a priori.The a priori study examines the suitability of the inertial manifold theory in modeling turbulent flows.The proposed AIM is investigated for determining the dimension of the inertial manifold or the attractor of the system and reconstructing the unresolved variables based on the information of the exact resolved field. 37Given these prior results, the focus here is on a posteriori validation, where the AIM approach is used to obtain the small-scale contribution to the evolution of the large resolved scales and close the turbulence modeling.As opposed to the a priori studies, this work examines the accuracy of the modeled trajectory of the system in the approximate inertial manifold.Moreover, a computationally efficient approach is provided, and possible simplifications to accelerate the model will be studied.The inertial manifold theory and the AIM methodology have been explained in Sec.II, followed by results in Sec.III.Finally, concluding remarks and future paths are discussed in Sec.IV.

II. MATHEMATICAL FORMULATION OF AIM
In this section, a reduced-order description of dynamical systems based on the IM theory is discussed.Any fluid system governed by a set of partial differential equations can be cast as a dynamical system after spatial discretization such that the discrete vector of variables of interest, e.g., momentum or energy in the case of nonreacting flows, is given by the set v = {v 1 , v 2 , . . ., vn v }, and the governing differential equations are written as where  is a discretized linear operator, which is taken to be positive and self-adjoint.Thus, the set of eigenvectors of  forms an orthonormal eigenbasis for the Hilbert space ℋ , in which the dynamics reside.The nonlinear term, ℱ , induces a computational challenge as it couples all scales of the solution.
The goal is to describe the dynamical features of the flow in a lower-dimensional manifold described by the dynamics of a subset of the variables of interest.An orthogonal projection operator P is defined to decompose the vector of variables (v) into the resolved u and unresolved w subsets.Applying the projection operator to the discrete governing equations, the evolution equations for the resolved and unresolved fields can be obtained as where Q = I − P is the complement operator of P, and it maps its operand to the null-space of the projection operator, P. Here, the projection P is taken to be onto the space formed by the first m eigenfunctions of the linear operator  .The main challenge in describing the large-scale dynamics only by using the resolved scales is the projected nonlinear term Pℱ (v), which cannot be computed directly from u.This is the typical closure challenge in turbulence modeling.The goal is to estimate w using only information of u and compute the projected nonlinear term.The approach here is to leverage the unresolved-scales governing equation by utilizing the IM approximation: the dynamics of u directly determine the dynamics of w.In other words, the components of w adjust to changes in u instantaneously.With the approximation dw/dt = 0, 39 Eq. (3) results in The above nonlinear equation can be solved by an iterative method to obtain a converged solution for w.With this approximation of the unresolved dynamics, the nonlinear term Pℱ (u, w) and the governing equations of the resolved modes u are closed.Before any spatial discretization, the unresolved subspace is infinite-dimensional.After spatial discretization, the unresolved subspace is the subspace between approximated IM and the entire state space.To develop a low-dimensional AIM, the unresolved subspace becomes considerably higher dimensional such that solving Eq. ( 4) can be adversely expensive in terms of cost and memory.This limitation can be lifted by considering only part of the unresolved dynamics that resides in a close neighborhood of the resolved subspace.Such simplification is in agreement with the exponential tracking property of dissipative dynamical systems, 29 and the resultant lower grid resolution reduces the cost of the approach.As aforementioned, the projection operator is positive and self-adjoint, with an ascending set of eigenvalues.Moving farther from the AIM, the timescales of the unresolved dynamics become smaller.By removing these exponentially decaying dynamics, the stiffness of the problem reduces, and a larger time step can be used.In the a priori analysis, 37 this modeling ansatz has been evaluated for the KSE and homogeneous isotropic turbulence (HIT).Direct numerical simulations are used to validate the accuracy of the IM approximation and study the convergence properties of the AIM approach.The purpose of this a posteriori study is to assess the AIM closure as a reducedorder model for the prediction of turbulent flows.Here, only the resolved scales will be evolved, and at each time step, the small scales will be reconstructed using the AIM approach.

III. AIM-BASED ROM FOR CANONICAL PROBLEMS
In this section, the system's dynamics are tracked in a lowdimensional AIM for the KSE and HIT.By tracking the dynamics of the system in a lower-dimensional space, convergence of the AIM model to the full-dimensional solution is shown for the KSE that possesses an IM.A computationally efficient framework of AIM is investigated for HIT because it is a more realistic problem of interest.A correction model is proposed and assessed in the forecast of statistical properties, dynamics of spatial statistics, and energy transfer between the resolved and unresolved subsets.Finally, the reducedorder model performances are compared against other prevalent turbulence models.

A. Kuramoto-Sivashinsky equation-based spatiotemporal chaos
The KSE is known to have a low-dimensional IM and has traditionally served as a surrogate problem for studying IM properties. 40,41The evolution of the KSE in spectral space is governed by where t is the time, q k = 2πk L , k ∈ Z, L is the spatial period, and μ is the viscosity.
Equation ( 5) can be arranged as Eq. ( 1) with the linear operator  = μq 4 k .The KSE has two linear operators, but only the diffusion operator (μq 4 k ) satisfies properties required by the theory of inertial manifolds (i.e., linear, unbounded, and self-adjoint).The full-dimensional system can be solved by truncating the discretized KSE at a sufficiently large wavenumber called nF.Here, the goal is to predict the dynamics of the system by evolving only the first m Fourier coefficients such that m ≪ nF, while the nonlinear contribution of the higher modes is modeled by the AIM approach.
After decomposition by projections P and Q, the dynamics of the system are predicted by time evolution of the resolved variables alone, u = (v −m/2 , . . ., v m/2 ), and the unresolved variables (w) are reconstructed by the AIM graph [Eq.( 4)], to close the governing equation of the resolved variables [Eq.( 2)].
Here, j denotes the iteration index in the approximation, and in the following analysis, j = 1, unless otherwise mentioned.The accuracy and convergence of Eq. ( 6) were assessed a priori in Ref. 37. The focus here is to examine the accuracy of AIM-based ROM in the prediction of the resolved dynamics.

Numerical results
The KSE exhibits spatiotemporal chaos, where infinitesimal perturbations can lead to exponential energy accumulation.The quadratic nonlinear term transfers energy from the low linearly unstable modes to the high modes with rapid exponential decay.Therefore, insufficient spatial resolution cannot capture energy dissipation and leads to numerical instability.The chaotic dynamics of the system are controlled by a bifurcation parameter defined as Re = L 2π √ ν .To identify a range of chaotic structures, the viscosity is kept constant at μ = 0.001, and the length of the domain is varied to have 158 ≤ Re ≤ 1000.As the range of scales increases with the Re number, the number of Fourier modes needed to resolve the dynamics also changes.At each Re number, there is a minimum resolution required for the stability of the solution (N min ), while a considerably higher resolution is needed to capture the spatiotemporal chaos of the KSE NDNS. 37The objective of the proposed reduced-order model is to evolve the dynamics of the system on an AIM spanned by m ≪ N min .When discussing the AIM dimension (m in Sec.II), the minimum resolution will be used as a reference.The initial condition of the dynamical system in physical space is chosen to be g(x) = sin(x) (1 + cos(x)) for the reference DNS simulation, and for the AIM prediction, the projected initial condition Pg(x) is used.For both full-dimensional and reducedorder simulations, the exponential time difference fourth-order Runge-Kutta method (ETDRK4) 42,43 with standard 3/2 de-aliasing is implemented.
Figure 1 (left) shows an example solution of the full-dimensional system for the Reynolds number 316.23 using NDNS = 4096 Fourier modes.At t > 0.15, the dynamics enter the chaotic regime, where a truncated system with 1024 Fourier modes for the same Reynolds number becomes unstable and blows up.However, including the AIM subgrid model in the evolution of the first 1024 Fourier modes stabilizes the solution and predicts transition to turbulence accurately (Fig. 1, right).
The accuracy of the model prediction depends on the size of the AIM, m, and the accuracy of the approximation of the unresolved dynamics [Eq.( 6)].The KSE is known to have a relatively low-dimensional inertial manifold that scales with Re. 36,40 Here, the model accuracy is assessed for a range of resolutions, m, and Fig. 2 shows the root mean square of the error between the AIM prediction and the full-dimensional system solution in physical space computed in the chaotic regime (t > 0.15), over a range of Re numbers and AIM dimensions.The AIM resolution is normalized by the integer part of the Reynolds number (Re), which is the number of linearly unstable modes.The AIM prediction converges uniformly to the exact solution when the approximate IM is large enough to contain all of the unstable dynamics.These results suggest that strong convergence properties can be obtained for resolutions exceeding this point, which is in agreement with prior works. 36,40n the AIM model, the ground assumption is that the unresolved variables respond instantly to the resolved dynamics, i.e., dw/dt = 0.The validity of this assumption and the rate of convergence to the fixed-point solution of Eq. ( 6) have been assessed in the a priori study. 37The optimum number of iterations depends on the Re number and the resolution of AIM (m).Seeking the fixed-point solution can be computationally expensive, and its feasibility should be judged based on the improved accuracy of the resolved dynamics prediction.Figure 3 (left) shows the energy spectrum of Fourier modes in the resolved and unresolved subspaces.The modeled resolved spectrum generally follows the DNS spectrum, but there are discrepancies at the largest resolved scales.The reconstructed unresolved spectrum can be improved by implementing more iterations in Eq. ( 6).The first-order approximation ( j = 1) considers only the nonlinear interaction between the resolved scales for the transfer of energy to the unresolved scales.This approximation is improved by seeking the fixed-point solution of Eq. ( 6) with more iterations to reconstruct the unresolved modes.However, unlike the unresolved dynamics, this higher-order approximation does not improve the resolved dynamics spectrum significantly.To assess this higher-order approximation more precisely, the probability density function (PDF) of the predicted resolved modes is compared for various closures obtained by different numbers of iterations ( j) in Eq. ( 6). Figure 3 (right) compares the PDF of the real part of the Fourier mode at the cut-off wavenumber for Re = 158.11and m = 256 predicted by AIM and DNS.It is shown that a higher-order approximation of the unresolved dynamics gives a better prediction of the resolved scales throughout the distribution.More details on the convergence properties of Eq. ( 6) are discussed in Ref. 37.

B. Homogeneous isotropic turbulence
In this section, the evolution of a turbulent flow in a domain of 2π × 2π × 2π m with periodic boundary conditions is modeled by AIM.The flow field is governed by the three-dimensional incompressible Navier-Stokes equations where ξ i is the velocity component in the ith direction, p is the hydrodynamic pressure, μ is the kinematic viscosity, and ρ is the density.Large-scale motions are forced by a constant linear forcing coefficient B to compensate for the viscous dissipation and reach statistical stationarity. 44,45By expanding the solution of Eq. ( 7) in Fourier space: ξi = ∑⃗ k vi( ⃗ k, t)e ⃗ k.⃗ x , a Galerkin projection of the governing equations leads to a system of ODEs that govern the evolution of the Fourier coefficients vi( ⃗ k, t), Equation ( 8) can be rearranged as Eq. ( 1) using the linear operator  = μ|k| 2 , with ℱ (v) containing all other terms.A sharp spectral projection operator with a cut-off wavenumber kc separates the resolved and unresolved subspaces such that all the modes with wavenumbers √ k 2 x + k 2 y + k 2 z ≤ kc are included in the resolved space.The number of modes satisfying this requirement is the dimension of AIM, m.With the IM assumption, the unresolved variables with wavenumber | ⃗ k| > kc can be approximated as where u is the vector of the resolved variables.It can be shown that the velocity field reconstructed by the AIM satisfies continuity. 37ith this approximation of the unresolved dynamics, the nonlinear term (Pℱ (u, w)) can be computed directly.Hence, the governing  6): red dashed-dotted line, and AIM with j = 3 in Eq. ( 6): green dashed line.
equation of the resolved dynamics [Eq.( 2)] is closed, and the resolved space dynamics can be evolved in time independently.
Unlike the KSE, the rate of convergence in Eq. ( 9) is slower than quadratic, 37 and more iterations are needed to approximate higher modes in the unresolved subspace, which can make the AIM approach inefficient.Removing the smallest unresolved scales from the computational grid makes the AIM approach more efficient in two ways: (1) by reducing the size of the domain, the cost of computing the nonlinear term ℱ (u, w) decreases and (2) small-scales require higher-order approximation of Eq. ( 9), and by removing them, fewer iterations are needed to approximate all unresolved modes.Therefore, the unresolved subspace can be decomposed into the unresolved but represented scales approximated by Eq. ( 9) and the unresolved and unrepresented scales with wavenumbers larger than the Nyquist wavenumber (kn g ) of the computational grid.The decomposition of the computational domain is shown in Fig. 4. The AIM dimension refers to m, which determines the degrees of freedom of the resolved subspace, and AIM resolution is the ng of an AIM simulation that includes both the resolved and represented unresolved subspaces.The effect of the unresolved and unrepresented scales on the dynamics of the system needs to be modeled.

A modified approximate inertial manifold
The a priori analysis has demonstrated that approximating only the larger unresolved scales reconstructs the interaction between FIG. 4. Representation of the resolved and unresolved subspaces in the wavenumber space.A circle with a radius km encloses the resolved modes.The light gray shaded area denotes the unresolved and unrepresented modes, and the dark gray area represents the unresolved but represented modes.kmax is the highest wavenumber in DNS calculations, and k N is the highest wavenumber in AIM-ROM calculations.The dashed red rectangle denotes the computational domain of AIM-ROM.
the resolved and unresolved dynamics with sufficient accuracy.However, it cannot provide adequate dissipation in the system. 37his behavior is similar to other approaches, which reconstruct the unresolved dynamics. 10,17Here, a dissipative modeling component similar to the eddy-viscosity approach is added to address this shortcoming.Most eddy-viscosity subfilter models in LES assume that the SGS contribution to the filtered field is dissipative and thus cannot predict the transfer of energy to the large scales (backscatter).An improvement has been made by adding additional non-dissipative terms to these models, such as the mixed model and the scalesimilarity model. 1,17Dynamic subfilter modeling can also account for backscatter in transitional flows if locally negative eddy-viscosity is allowed. 13On the contrary, AIM recovers the nonlinear interaction between the resolved and unresolved scales by reconstructing the subfilter field without assuming a forward cascade of energy.The recovered unresolved energy spectrum can be used to determine the rate of backward/forward transfer of energy between the resolved and unresolved subsets.The energy stored in the unresolved scales either transfers backward to the resolved scales or transfers forward to the smaller unresolved scales where it finally gets dissipated by molecular viscosity.If only a subset of the unresolved space is reconstructed by the AIM approach, the energy of the unrepresented scales should be drained from the system.The energy of unrepresented scales can be approximated by the forward scatter of energy of the reconstructed scales as this energy will eventually transfer to smaller unrepresented scales.Therefore, energy virtually transferred to the unresolved and unrepresented scales is dissipated.First, the AIM approximation is studied to determine whether backward/forward energy transfer is captured accurately.The subgrid-scale dissipation, εSGS, is defined as where τij = P(vivj) − uiuj is the SGS stress, and Sij = ( ∂u i ∂x j + ∂u j ∂x i ) is the resolved-scales strain rate.When the unresolved scales remove energy from the resolved ones (forward scatter), εSGS is negative; and if SGS transfers energy to the resolved scales (backscatter), SGS dissipation is positive.Therefore, forward and backward energy transfer can be defined as 46 ε− = 0.5(εSGS − |εSGS|), ε+ = 0.5(εSGS + |εSGS|). ( The energy from the forward scatter will eventually dissipate at the smallest scales.If the smallest unresolved scales are discarded, this energy needs to be removed to avoid energy accumulation beyond the cut-off wavenumber.The energy spectrum reconstructed by AIM can provide the rate of energy dissipation at the unresolved scales.A dynamic spatially varying viscosity can be determined to dissipate the forward cascade of energy beyond the cut-off wavenumber, where ε− is the Fourier transfer of the forward cascade of energy, and E(k) is the energy spectrum that represents the contribution to the turbulent kinetic energy 1 2 ⟨vivi⟩ from all modes with | ⃗ k| in the 1: At t = tn: 2: μ T ( ⃗ k) = 0 3: w 0 = 0 4: for j = 1 : l do 5: Compute ε− by Eq. ( 11) 9: Compute μ T (k) by Eq. ( 12) 10: end for 11: Compute Pℱ (un, wn) 12: Advance resolved scales by Eq. ( 2) Here, the energy spectrum is computed for the resolved and reconstructed unresolved modes.With this turbulent viscosity, the effective viscosity at the unresolved scales is , and the unresolved dynamics are approximated by It should be noted that μ eff changes at each iteration.For j = 1, μ eff = μ because the unresolved subspace is not reconstructed yet.Algorithm 1 summarizes the modified AIM-ROM approach that advances from time step t = tn to t n+1 .By implementing more iterations (increasing l in Algorithm 1), the modified viscosity is updated with the reconstructed unresolved scales.Backward and forward scattering of energy occurs at all scales, and turbulent viscosity obtained from Eq. ( 12) can be defined at both the resolved and unresolved scales.Three different approaches have been considered: (1) modifying viscosity only at the resolved subspace similar to the optimal LES model, 47 (2) modifying viscosity at both the resolved and unresolved subspaces, and (3) modifying viscosity only at the unresolved subspace.Here, the effective viscosity is modified only at the represented unresolved scales as the forward energy transfer from the resolved scales will eventually dissipate at the represented unresolved scales.Hence, modifying the effective viscosity at the resolved scales based on ε− introduces too much dissipation in the system.It should be mentioned that while modifying the effective viscosity at the unresolved scales changes the linear operator upon which separation of scales is based, it is aligned with the inertial manifold theory requirements as it increases the spectral gap of the linear operator of the Navier-Stokes equations.As unresolved scales become more dissipative, the separation of scales between the resolved and unresolved scales is more prominent.In turn, the assumption that unresolved scales equilibrate to the AIM dynamics is more justified.

Numerical results
Direct numerical simulations of HIT at two spatial resolutions are used to investigate the accuracy of the AIM prediction.For DNS and AIM-ROM simulations, a pseudo-spectral method with dealiasing is used for the non-linear term.Exact time integration is used for the linear viscous term, and second-order Runge-Kutta and the flow statistics are monitored for several eddy turnover times (τ) to ensure that the forcing method does not lead to instability and energy pile-up at small scales.
Current theories cannot prove the existence of an inertial manifold for the Navier-Stokes equations. 249][50] For instance, the attractor dimension is lower than the number of degrees of freedom required by DNS for turbulent flows. 27,51More specifically, in forced HIT, it is shown that the attractor dimension scales with ( L η ) , 27 where L is the domain length.For the turbulent fields considered here, this estimated attractor dimension (DKY ) is used as a reference for the assessment of AIM accuracy over a range of AIM dimensions.It should be mentioned that when an inertial manifold exists, it contains the attractor of the system.Hence, the AIM should be larger than the attractor of the system regardless of the existence of an inertial manifold.

FIG. 6.
Top: subgrid-scale dissipation normalized by the total resolved dissipation, middle: SGS backscatter normalized by the total resolved dissipation, and bottom: fraction of points with backscatter of the energy in the computational domain.DNS of 256 3 : blue solid line with circle, AIM modeling for 256 3 DNS case: red dotted line with circle.DNS of 512 3 : green solid line with square, and AIM modeling for 512 3 DNS case: purple dotted line with square.The horizontal axis is the resolved subspace dimension normalized by the full-dimensional system dimension (m/ng).
To evaluate the proposed AIM model against common turbulence modeling approaches in terms of accuracy and efficiency, large-eddy simulations with the dynamic Smagorinsky subfilter model 11 have been conducted.For the sake of comparison, the sharp spectral filter is used for the LES to make the resolved subspace of AIM and LES almost identical.However, the resulting variable separations are not identical as the resolved subspace of AIM in the wavenumber space is a sphere with radius kc, while the same cut-off wavenumber in LES resolves all the wavenumbers enclosed in a cube of side length 2kc.Following the conventional LES practices, the subfilter field representation is implicit such that the subfilter scales are not represented on the computational grid, and their contribution is modeled by the dynamic eddy viscosity.On the other hand, AIM reconstructs the subfilter field either entirely or just a subspace of it.Therefore, for the same cut-off wavenumber, the computational grid is larger in AIM compared to the LES.
First, the AIM approximation [Eq.( 9)] has been examined to see if the reconstructed turbulent field captures the forward and backward scatter of energy between the resolved and unresolved scales accurately.Figure 5 shows the rate of energy transfer over the range of scales for the 256 3 field with kc = 16.As expected, subgridscale energy transfer is dominated closer to the cut-off wavenumber and at larger unresolved scales.This behavior confirms that there is no need to recover all of the unresolved subspace, and approximating only the largest unresolved scales is sufficient to capture subgrid-scale effects on the resolved dynamics.It can be seen that AIM captures energy transfer in both directions.However, it overestimates the backward scatter at the unresolved scales.This may be due to the limited approximation of the unresolved scales, and by implementing more iterations in Eq. ( 9) and recovering smaller unresolved scales, this approximation can improve. 37tatistics describing the energy transfer between the resolved and unresolved subspaces are provided in Fig. 6, where the (SGS) dissipation (top) and energy backscatter of the subfilter field (middle) as a function of the normalized AIM dimension (m/ng) are shown.Here, the SGS dissipation and backscatter of energy to the resolved scales are computed from turbulence fields modeled by the AIM over a range of resolutions (m).Also, DNS fields are filtered, and the exact values of these quantities are computed at different filtering widths for comparison with AIM.At each cut-off wavenumber, the SGS dissipation (⟨εSGS⟩) and the energy backscatter (⟨ε+⟩) are normalized by the total resolved dissipation (⟨ε⟩).It can be seen that by increasing the filter width, i.e., by using a lower AIM resolution, the amount of SGS dissipation increases.Accordingly, the amount of backscatter of energy to the resolved scales increases, because the cut-off wavenumber is farther away from the rapidly dissipative scales, and a larger part of the inertial range is in the unresolved subspace.The number of points in the physical domain  experiencing energy backscatter is almost independent of the cutoff wavenumber (Fig. 6, bottom), which shows that even when the amount of energy backscatter is not considerable compared to the total dissipation, subgrid energy backscatter occurs between the smallest scales at the dissipation range.AIM predicts the same characteristics but more locations in the domain experience backscatter.This is not surprising as the AIM approach [Eq.( 9)] models nonlinear interaction between the scales but does not have a completely dissipative component.Overall, these results show that AIM can capture energy transfer between the resolved and unresolved subspaces accurately, and this property can be used to implement a dynamic dissipative component to account for the unresolved and unrepresented scales [Eq.(13)].Turbulent statistics predicted by original and modified AIM models are compared against the DNS calculation in Fig. 7, where turbulent kinetic energy is computed based on the resolved field and the total dissipation rate is computed based on the effective viscosity in each of the modeling approaches including LES.It can be observed that the original AIM approximation alone cannot predict enough dissipation in the system, but adding a dissipative component solves this problem, and dissipation is almost accurately predicted by the modified AIM model.This improvement in AIM prediction is not dominant in the resolved turbulent kinetic energy as the modified viscosity removes energy only from the unresolved scales.These statistics show that AIM models outperform the dynamic Smagorinsky approach.However, it should be noted that here all of the unresolved subspace is approximated by AIM [Eqs.(9) and ( 13  To assess the modified AIM approach in modeling the effect of the unresolved and unrepresented scales, only a part of the unresolved subspace is kept in the computational domain and approximated by AIM [Eq.( 13)]. Figure 8 shows the statistical properties for the 256 3 field and kc = 16, where the unresolved subspace dimension in AIM is reduced by discarding higher wavenumbers and using a lower grid resolution.For the full-dimensional system of this case, the highest wavenumber in the computations is 120.7, and the full-dimensional unresolved subspace contains all modes with 16 < | ⃗ k| < 120.7.By reducing the AIM resolution to N g = 64 3 and N g = 128 3 , the highest unresolved modes considered are roughly 30 and 60, respectively.While reducing the unresolved subspace dimension does not affect the turbulent kinetic energy of the resolved subspace substantially, it underestimated total dissipation in the system considerably.In the lowest-resolution AIM for kc = 16, the largest wavenumber reconstructed by AIM is 30.2, which gives | ⃗ k|/kc < 2. Figure 5 shows that a considerable amount of forwarding cascade of energy is discarded at this AIM resolution.In this case, the dimension of AIM, m, is only 0.75 of the estimated dimension of the system's attractor (DKY in Table I).By increasing the AIM dimension, the approximation of the unresolved dynamics and AIM prediction improve considerably.Considering approximate inertial manifolds for the 512 3 case with two different dimensions obtained from kc = 64 and kc = 128 results in m/DKY ≈ 0.5 and m/DKY ≈ 4. Figure 9 compares the time evolution of turbulent kinetic energy and dissipation rate predicted by AIM against the DNS data.By increasing the AIM dimension, both approximations have improved especially earlier in the prediction time.
The resolved energy spectrum of the 512 3 field predicted by the DNS, AIM, and LES models for projection wavenumber kc = 64 is compared in Fig. 10 (left).At this resolution, the AIM dimension, m, is almost half of the estimated dimension of the system's attractor (DKY in Table I).The energy spectrum predicted by AIM models is more accurate than the dynamic Smagorinsky model.The original AIM model overestimates the energy of the smallest resolved scales.This issue has been alleviated in the modified AIM model.The spectrum predicted by the dynamic Smagorinsky model is quite different from the exact spectrum, which can be due to the limited forcing of the large scales.The dynamic Smagorinsky approach relies on the scale-similarity between larger resolved scales, smaller resolved scales, and unresolved scales.Limited forcing at the larger resolved scales can falsely impose a higher rate of energy transfer at the smaller scales and lead to overshooting of energy at the smallest resolved scales.To test this explanation, the same initial condition for the 512 3 field in Table I is used for decaying HIT where there is no force.Figure 10 (right) shows the energy spectrum predicted by the DNS, modified AIM, and LES at t/τ ≈ 1, and it can be seen that the LES spectrum follows the exact spectrum even though it is more dissipative at the smallest resolved scales.This behavior of the LES and AIM models is consistent over the range of cut-off wavenumbers considered.However, the prediction of the dynamic Smagorinsky model improves when the cut-off wavenumber is high enough to be in the dissipation range of the energy spectrum.
Finally, the contours of the velocity field predicted by the DNS, modified AIM, LES models are compared in Fig. 11.Here, the decaying case of the 512 3 field has been chosen for comparison.The decaying turbulent field allows for a fixed time step for all simulations.Velocity fields are compared at t/τ 0 ≈ 1.Both modeled fields look quite similar to the DNS field, but it can be seen that AIM preserves more details of the smaller structures.It should be mentioned that in decaying HIT, as turbulent energy dissipates, the size of the attractor changes and shrinks.Hence, the approximation of AIM becomes more accurate for longer prediction times.The computational cost of the AIM and LES models are compared over a range of cut-off wavenumbers in Fig. 12, where computational costs of the AIM and LES models are normalized by the cost of the corresponding DNS.The comparison is based on the reduced grid resolution of the models (N/N g ), where N is the grid resolution of AIM or LES simulations, and N g is the grid resolution of the corresponding DNS.It is shown that AIM is more expensive than the constant Smagorinsky model, but it is more efficient compared to the dynamic LES modeling.It should be mentioned that for the same grid resolution of AIM and LES, the resolved space in AIM is lower-dimensional than in LES, and it is not possible to compare the computational cost of these models with the same accuracy.

IV. CONCLUSIONS
A reduced-order model of turbulent flows has been developed based on the inertial manifold (IM) theory.Casting the discretized governing equations as a dynamical system provides a path for the decomposition of variables without relying on traditional scale-separation methods, such as spatial filtering.Here, governing equations of the system have been leveraged to define resolved variables and recover unresolved variables to directly compute the nonlinear term.The proposed model has been examined on two canonical flows: the one-dimensional KSE and the three-dimensional HIT.
The existence of an inertial manifold is not yet proven for turbulent flows.Nevertheless, the construction of an AIM for the Navier-Stokes equations shows promising results.The AIM-based reduction requires that the dimension of the reduced-order model should be higher than the dimension of the attractor.Since exactly obtaining the attractor dimension is not feasible for most practical problems, the proposed AIM is examined over a range of dimensions to understand the validity of this approximation.Convergence properties of the AIM conform with direct estimations of the size of the attractor for these systems, proving that the proposed AIM can approximate the dynamics of the attractor.Studying the attractor of chaotic systems provides new paths for the development of reduced-order models to predict and control complex systems.Direct methods for finding the topology of the attractor are prohibitively expensive, but strong convergence properties observed for approximate inertial manifolds over a range of problems considered in this study show the potential of this approach in locating the attractor of more practical systems.The next steps will involve extensions to non-homogeneous systems such as wall-bounded flows, for which the AIM approach should be cast in physical space. 52or a given resolved field, AIM reconstructs a single realization of the unresolved dynamics.This contribution of an AIM to the resolved dynamics can also be seen as a subgrid-scale model.In all configurations, for a sufficiently large dimension of the AIM, the unresolved dynamics were found to respond to the dynamics of the AIM instantaneously.However, smaller scales in the unresolved dynamics are less responsive to the dynamics of the IM, and there is a time delay in their response.A higher-order estimation of the unresolved dynamics, where the interactions between the resolved and unresolved dynamics are included, improves the AIM estimation of the unresolved dynamics.The rate of convergence is controlled by the nonlinear interaction between the resolved and unresolved scales.However, turbulence is broadband, and the approximation of the unresolved dynamics farther from the approximate inertial manifold can be cost-prohibitive.It is shown that reconstruction of the entire unresolved subspace is not necessary, and recovering the unresolved dynamics in the vicinity of the AIM captures the nonlinear interaction sufficiently.The information recovered by AIM is used to model the effect of the dynamics far from the AIM.The modified AIM approach is robust, efficient, and more accurate in the prediction of statistical properties of the system.The modified model shows the capacity of the AIM approach for an adaptive modeling framework.

FIG. 1 .FIG. 2 .
FIG. 1. Solution of the KSE for Re = 316.23.Left: DNS with N DNS = 4096 and right: AIM with m = 1024.Only part of the computational domain is shown.

FIG. 5 .
FIG. 5. Spectrum of the backward (blue solid line) and forward (black solid line) SGS dissipation rate for ng = 256 3 and kc = 16.DNS: solid lines and AIM: dashed lines.

FIG. 7 .
FIG. 7. Time evolution of the resolved turbulent kinetic energy (left) and the dissipation rate (right) of the 256 3 field with kc = 16.DNS: black solid line, AIM: red dashed line, modified AIM: green dashed-dotted line, and LES: violet dotted line.

FIG. 8 .
FIG. 8. Time evolution of the resolved turbulent kinetic energy (left) and the total dissipation rate (right) of the 256 3 field with kc = 16.DNS: black solid line, modified AIM with Ng = 256 3 : green dashed-dotted line, modified AIM with Ng = 128 3 : red dashed-dotted line, modified AIM with Ng = 64 3 : blue dashed-dotted line, and LES: violet dotted line.

FIG. 9 .
FIG. 9. Time evolution of the resolved turbulent kinetic energy (left) and the total dissipation rate (right) of 5123 field with kc = 64 (blue lines) and kc = 128 (black lines).The solid lines are obtained from DNS, and the dashed lines are predicted by a modified AIM model.Since the total dissipation does not depend on kc, only one line is shown here for the dissipation evolution of the DNS field (right).

FIG. 10 .
FIG. 10.Left: the resolved kinetic energy spectrum of a forced 512 3 field for kc = 64, DNS: black solid line, AIM: red dashed-dotted line, modified AIM: green dashed-dotted line, and LES: violet dotted line.Right: the resolved kinetic energy spectrum of the decaying 512 3 field for kc = 64, DNS: black solid line, modified AIM: red dashed-dotted line, and violet dotted line.

FIG. 12 .
FIG.12.Reduced computational cost of DNS for Ng = 256 3 with AIM: blue filled circle, dynamic Smagorinsky: blue filled up-pointing triangle, and constant Smagorinsky: blue filled square models.Reduced cost of DNS for Ng = 512 3 with AIM: red filled circle, dynamic Smagorinsky: red filled up-pointing triangle, and constant Smagorinsky: red filled square models.

TABLE I .
Numerical setup for HIT.RK2) is used for the other terms.TableIpresents the numerical setup for DNS cases.The Taylor microscale Reynolds number Re λ and the Kolmogorov length scale η (Δx is grid spacing in each direction) are monitored over the initialization time to make sure the turbulent field is fully developed and resolved.The forcing of the velocity field has been limited to the large energy-containing scales [B( ⃗ (