Recent Advancements in Large Eddy Simulations of Compressible Real Gas Flows

Abstract

This chapter explores some of the recent advancements in the field of computational fluid dynamics, specifically with respect to large eddy simulations. We start by introducing some fundamental concepts of scales of turbulent fluid flows. We then discuss the available simulation methods, major challenges of each method and the advantages of large eddy simulation (LES) over other methods. The focus of the chapter then shifts to real gas flows. Governing equations for LES of compressible flows and modifications for real gas flows are presented. Some of the filtering methods are discussed including physical space, Fourier space and proper orthogonal decomposition (POD) based filtering. The chapter then focuses on modeling of sub-filter / subgrid scale terms. A brief discussion of traditional modeling approaches is provided followed by a discussion of the current research and advancements. This discussion includes stochastic modeling using the filtered mass density function (FMDF) approach and machine learning based models for subgrid scales. For each method, a brief a background of the method, their applicability to different scenarios, their advantages, and disadvantages are presented.

1. Introduction

Turbulence in fluid flows is a three dimensional phenomena characterized by chaos. However, over a century of research in the field has shown that patterns exist in such chaos. Turbulence occurs at large Reynolds numbers, where the inertial forces of the fluid are larger than the viscous forces, causing vortices (eddies) to be generated due to gradients in velocity. These vortices rapidly bend, stretch and collapse into smaller eddies. As a result there exist multiple length scales of turbulent flows with vortex sizes ranging from a few kilometers (in atmospheric turbulence) to a few hundred nanometers. Further, an energy cascade phenomenon is observed where the largest scales of the flow that contain the most energy, break down into smaller scales with lesser energy. This continues until the smallest realizable scale called the Kolmogorov scale [1], beyond which the energy dissipates as heat (Figure 1). The largest flow scales are observed to be dependent on the physical constraints of the flow (initial and boundary conditions), while some degree of universality in structure and statistics is observed in the smaller scales. In computational fluid dynamics (CFD) three broad simulation techniques exist. The first is called Direct Numerical Simulation (DNS), where the computational domain is set such that the simulation captures all the length and time scales of the flow. This approach, solves the full form of Navier Stokes equations by setting the size of mesh, equivalent to the Kolmorogov scales and the number mesh points equivalent to the size of the largest eddy to capture all the length scales [2]. Needless to say, this approach is extremely computationally intensive and impractical for a majority of flows of interest. The second approach is called Reynolds Averaged Navier Stokes (RANS), where an averaging operation is applied to the flow variables in the governing equations. This results in generation of mean flow statistics where turbulence at all scales are modeled [3]. While computationally inexpensive, the effects of turbulence are oversimplified. The third approach is called Large Eddy Simulation (LES), where a filtering operation is applied to the flow variables in the governing equation, such that the largest scales of the flow are captured by the simulation, while the smaller (universal) scales are modeled. This approach results in a more accurate representation of turbulence, while being applicable to many complex flows of interest.

Figure 1.

Spectra of turbulent flow: This figure illustrates the length scales (from integral $l$ to Kolmogorov $\eta$) of turbulent flows and the energy cascade phenomena. The turbulent kinetic energy (TKE) in inertial sub-range typically exhibits a universality $E\left(k\right)\sim k^{-5/3}$. LES cutoff filter $G\left(k\right)$ separates the resolved and modeled scales.

Research in LES focuses primarily on modeling the flow physics occurring at the length scales that are not explicitly captured by the simulation. These scales are called sub-filter or subgrid scale (SGS). Unlike an averaging operation that decomposes the statistics of the flow into mean and fluctuating parts, the flow characteristics of subgrid scales are dependent entirely on the size of the filter used in the simulation. In a sense, a trade-off between the accuracy and computational expense is set by the filter size. Over many decades, a rich collection of models have been developed for the SGS fluctuations occurring in various flow scenarios including incompressible flows, compressible flows, chemically reacting flows, and boundary layer flows [4, 5]. While modeling turbulence in incompressible flows is a challenge, the difficulty increases when compressibility is introduced. From a computation point of view, there are additional equations that need to be solved to account for density variations. The boundary conditions must be designed properly to allow proper dissipation of generated acoustic waves and noise without reflection. Further, the flow physics at the subgrid scales become more localized. Modeling difficulty increases in flows with chemical reactions since the flow physics vary over a wider range of length and time scales. For instance in combustion, the flame fronts are very thin compared to the rest of the flow scales yet strong interactions exist between these scales [6], and the rates of many chemical reactions of interest are in order of a few hundred nano-seconds. These challenges are further exacerbated when the ambient pressure is increased from $1\mathit{atm}$ as many modern day combustion devices such as rocket engines, diesel engines and gas turbines operate at elevated ambient pressures $\left(10\mathit{atm}<P<100\mathit{atm}\right)$. At these pressures many of the simplifications such as the ideal gas equation of state, constant property models or simple forms of heat and mass flux models, are no longer applicable. Various regimes of flow fields are encountered such as large ranges of chemical time scales, transport rates, Reynolds and Mach numbers with weak and strong compressibility effects. Thus, the models developed must be consistent with the fundamental flow physics. This chapter provides a high level overview of a few modeling techniques in literature that are applicable to LES of compressible real gas flows. This chapter is not intended to be a detailed review of all the existing applicable methods.

2. Governing equations

In LES, filters of prescribed widths are applied to the primitive variables in the governing equations. A filtering operation of the form,

$$\overline{\mathrm{\Phi }}={\int }_{\mathrm{\Omega }}\mathrm{\Phi }\left(x^{\prime },t\right)G\left(x^{\prime },x\right)d{x}^{\prime },$$E1

is applied to the variables over the domain, $\mathrm{\Omega }$. In the above equation $G\left(x\right)$ denotes the filter function. The term $\overline{\mathrm{\Phi }}$ represents the filtered (resolved) value of variable and ${\mathrm{\Phi }}^{\prime }=\mathrm{\Phi }-\overline{\mathrm{\Phi }}$ the fluctuating subgrid component. For compressible flows, a density weighted Favre filtering operation is performed, $\tilde{\mathrm{\Phi }}=\overline{\rho \mathrm{\Phi }}/\overline{\rho }$, with ${{\mathrm{\Phi }}^{\prime }}^{\prime }=\mathrm{\Phi }-\tilde{\mathrm{\Phi }}$. For reacting flows, the resulting governing equations for mass, momentum, energy, and scalar conservation equations, respectively can be expressed as,

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial \overline{\rho }\widetilde{u_i}}{\partial x_i}=0,$$E2

$$\begin{array}{l}\frac{\partial \overline{\rho }\widetilde{u_i}}{\partial t}+\frac{\partial \overline{\rho }\widetilde{u_i}\widetilde{u_j}}{\partial x_j}=-\frac{\partial P\left(\overline{\mathrm{\Psi }}\right)}{\partial x_i}-\frac{\partial }{\partial x_i}\left[\overline{P\left(\mathrm{\Psi }\right)}-P\left(\overline{\mathrm{\Psi }}\right)\right]+\frac{\partial {\tau }_{\mathit{ij}}\left(\overline{\mathrm{\Psi }}\right)}{\partial x_j}-\frac{\partial }{\partial x_i}\left[\overline{{\tau }_{\mathit{ij}}\left(\mathrm{\Psi }\right)}-{\tau }_{\mathit{ij}}\left(\overline{\mathrm{\Psi }}\right)\right]-\frac{\partial T_{\mathit{ij}}}{\partial x_j},\end{array}$$E3

$$\begin{array}{l}\frac{\partial \overline{\rho }\widetilde{e_t}}{\partial t}+\frac{\partial \overline{\rho }\widetilde{u_j}\widetilde{e_t}}{\partial x_j}=-\frac{\partial P\left(\overline{\mathrm{\Psi }}\right)\widetilde{u_j}}{\partial x_i}-\frac{\partial }{\partial x_i}\left[\overline{P\left(\mathrm{\Psi }\right)u_j}-P\left(\overline{\mathrm{\Psi }}\right)\widetilde{u_j}\right]+\\ \frac{\partial \widetilde{u_i}{\tau }_{\mathit{ij}}\left(\overline{\mathrm{\Psi }}\right)}{\partial x_j}-\frac{\partial }{\partial x_i}\left[\overline{u_i{\tau }_{\mathit{ij}}\left(\mathrm{\Psi }\right)}-\widetilde{u_i}{\tau }_{\mathit{ij}}\left(\overline{\mathrm{\Psi }}\right)\right]-\frac{\partial Q_j\left(\overline{\mathrm{\Psi }}\right)}{\partial x_j}-\frac{\partial }{\partial x_j}\left[\overline{Q_j\left(\mathrm{\Psi }\right)}-Q_j\left(\overline{\mathrm{\Psi }}\right)\right]+\\ \frac{\partial }{\partial x_j}\left[\sum _{\gamma =1}^{N_s}{H}_{,\gamma }\left(\overline{\mathrm{\Psi }}\right)J_{j,\gamma }\left(\overline{\mathrm{\Psi }}\right)\right]-\frac{\partial }{\partial x_j}\left[\sum _{\gamma =1}^{N_s}\overline{H_{,\gamma }\left(\mathrm{\Psi }\right)J_{j,\gamma }\left(\mathrm{\Psi }\right)}-\sum _{\gamma =1}^{N_s}{H}_{,\gamma }\left(\overline{\mathrm{\Psi }}\right)J_{j,\gamma }\left(\overline{\mathrm{\Psi }}\right)\right]-\\ \frac{\partial }{\partial x_j}\left(\overline{\rho }\widetilde{u_j{e}_t}-\overline{\rho }\widetilde{u_j}\widetilde{e_t}\right)+\overline{S_e},\end{array}$$E4

$$\begin{array}{l}\frac{\partial \overline{\rho }\widetilde{Y_{\gamma }}}{\partial t}+\frac{\partial \overline{\rho }\widetilde{Y_{\gamma }}\widetilde{u_j}}{\partial x_j}=-\frac{\partial J_{j,\gamma }\left(\overline{\mathrm{\Psi }}\right)}{\partial x_j}-\frac{\partial }{\partial x_j}\left[\overline{J_{j,\gamma }\left(\mathrm{\Psi }\right)}-J_{j,\gamma }\left(\overline{\mathrm{\Psi }}\right)\right]-\\ \frac{\partial }{\partial x_j}\left[\overline{\rho }\widetilde{Y_{\gamma }{u}_j}-\overline{\rho }\widetilde{Y_{\gamma }}\widetilde{u_j}\right]+\overline{S_{Y_{\gamma }}}.\end{array}$$E5

In these equations, $\mathrm{\Psi }$ represents the primitive variables (density $\rho$, velocities $u_i$, energy $e_i$, and mass fraction of species $Y_{\gamma }$), $\overline{\mathrm{\Psi }}$ represents the filtered values of the primitive variables and $\mathrm{\Phi }\left(\overline{\mathrm{\Psi }}\right)$ represents the variables that are calculated from the filtered primitive variables. For instance, $\overline{P\left(\mathrm{\Psi }\right)}$ represents the filtered pressure, whereas $P\left(\overline{\mathrm{\Psi }}\right)$ represents the pressure calculated from the filtered primitive variables available from LES and the equation of state. Additional closure terms are required for the source terms $\overline{S_e},\overline{S_{Y_{\gamma }}}$, the viscous shear stress tensor, the heat and the mass flux vectors. The gravity term typically appearing in the Navier Stokes equations is neglected for high pressure flows.

2.1 Modifications for real gas flows

For real gas flows at elevated pressures, a few terms in the governing equations are modified. A brief summary of the modifications is provided in this section. Complete formulation of these equations are presented in the works of Palle [7], Foster and Miller [8].

2.1.1 Equation of state

While ideal gas equation of state is predominantly used for compressible flows at atmospheric pressures, a real gas equation is required for higher pressures to account for effects of compressibility, non-equilibrium thermodynamics, variations in material properties and effects of inter-molecular forces. A cubic form of Peng Robinson equation of state is typically used for closure of pressure term in real gases due to its simplicity.

$$P\left(\overline{\mathrm{\Psi }}\right)=\frac{\mathit{RT}\left(\overline{\mathrm{\Psi }}\right)}{v\left(\overline{\mathrm{\Psi }}\right)-B_m\left(\overline{\mathrm{\Psi }}\right)}-\frac{A_m\left(\overline{\mathrm{\Psi }}\right)}{v{\left(\overline{\mathrm{\Psi }}\right)}^2+2v\left(\overline{\mathrm{\Psi }}\right)B_m\left(\overline{\mathrm{\Psi }}\right)-B_m{\left(\overline{\mathrm{\Psi }}\right)}^2},$$E6

where, $A_m\left(\overline{\mathrm{\Psi }}\right),B_m\left(\overline{\mathrm{\Psi }}\right)$ represent resolved mixture parameters, $R$ the molar gas constant and $v\left(\overline{\mathrm{\Psi }}\right),T\left(\overline{\mathrm{\Psi }}\right)$ the molar volume and temperature calculated from the resolved primitive variables.

2.1.2 Property models

Flow properties at atmospheric pressures are typically assumed as constant or simple functions of temperature. However, for real gas flows with the aforementioned equation of state, realistic property models related to the equation of state through the Gibbs function must be derived. These property models include departure functions that represent the thermodynamic deviations of the properties from those at low pressures. For instance, modified the molar heat capacity term can be expressed as,

$$C_p\left(\overline{\mathrm{\Psi }}\right)=C_p^0\left(\overline{\mathrm{\Psi }}\right)+f\left(\mathrm{\Phi }\left(\overline{\mathrm{\Psi }}\right)\right),$$E7

where $C_p^0\left(\overline{\mathrm{\Psi }}\right)$ represents the low pressure heat capacity and $f\left(\mathrm{\Phi }\left(\overline{\mathrm{\Psi }}\right)\right)$ the deviation term modeled as a function of resolved variables. Similarly, property models for molar specific enthalpy $H\left(\overline{\mathrm{\Psi }}\right)$, partial molar volume $v\left(\overline{\mathrm{\Psi }}\right)$, partial molar enthalpy $H_{,\gamma }\left(\overline{\mathrm{\Psi }}\right)$, and the mixture molar internal energy $u\left(\overline{\mathrm{\Psi }}\right)=H\left(\overline{\mathrm{\Psi }}\right)-P\left(\overline{\mathrm{\Psi }}\right)v\left(\overline{\mathrm{\Psi }}\right)$ must be derived. The mass diffusion coefficients, mixture viscosity and thermal conductivity are typically obtained from look up tables.

2.1.3 Diffusion terms

In flows at atmospheric pressures, the Fourier form of heat diffusion $\left(Q_j=-\alpha \partial T/\partial x_j\right)$ and Fickian form of mass diffusion $\left(J_j=-D_m\partial Y_{\gamma }/\partial x_j\right)$ models are generally applicable. However, for flows at elevated pressures, generalized multi-component forms of the heat and mass flux vectors derived from non-equilibrium thermodynamics and fluctuation theory [9, 10, 11] are required. The modifications must include the effects of Soret, Dufour and cross diffusion for higher pressures. The general form for heat diffusion can be expressed as,

$$\begin{array}{l}{Q}_j\left(\overline{\mathrm{\Psi }}\right)=-f\left(\alpha ,D_m,\mathrm{\Phi }\left(\overline{\mathrm{\Psi }}\right)\right)\frac{\partial T\left(\overline{\mathrm{\Psi }}\right)}{\partial x_j}\\ -f\left(\alpha ,D_m,\mathrm{\Phi }\left(\overline{\mathrm{\Psi }}\right)\right)\frac{\partial P\left(\overline{\mathrm{\Psi }}\right)}{\partial x_j}\\ -f\left(\alpha ,D_m,\mathrm{\Phi }\left(\overline{\mathrm{\Psi }}\right)\right)\frac{\partial X_{\gamma }\left(\overline{\mathrm{\Psi }}\right)}{\partial x_j},\end{array}$$E8

where, $\alpha ,D_m$ represent the thermal and mass diffusivities, $\mathrm{\Phi }\left(\overline{\mathrm{\Psi }}\right)$ the filtered variables including molar fluxes $\left(X_{\gamma }\right)$, molecular weights and mixture parameters. The cross diffusion effects are incorporated using non equal values of species mass diffusion coefficients. The Soret effect is associated with terms proportional to $\partial T\left(\overline{\mathrm{\Psi }}\right)/\partial x_j$ and $\partial X_{\gamma }\left(\overline{\mathrm{\Psi }}\right)/\partial x_j$, while the Dufour effect is associated with terms proportional to $\partial P\left(\overline{\mathrm{\Psi }}\right)/\partial x_j$ and $\partial X_{\gamma }\left(\overline{\mathrm{\Psi }}\right)/\partial x_j$. A similar form is used for the mass diffusion term. These terms produce similar effects as the Fourier and Fickian models at low pressures.

2.1.4 Reaction terms

Chemical kinetics of varying level of complexities can be included in simulation of reacting flows. However, to strike a balance between computational efficiency and accuracy, reduced form of pressure dependent kinetics must be used for reacting flows at elevated pressures. Pressure dependent mechanisms for high pressure reacting flows are proposed in the studies of Sohn et al. [12] and Chung et al. [13].

2.2 Unclosed terms

The governing equations (Eqs. (3)-(5)) contain many unclosed subgrid terms that require modeling.

• $T_{\mathit{ij}}=\overline{\rho }\widetilde{u_i{u}_j}-\overline{\rho }\widetilde{u_i}\widetilde{u_j}$—subgrid turbulent stress

• $\overline{S_e},\overline{S_{Y_{\gamma }}}$—filtered reaction source

• $\left[\overline{{\tau }_{\mathit{ij}}\left(\mathrm{\Psi }\right)}-{\tau }_{\mathit{ij}}\left(\overline{\mathrm{\Psi }}\right)\right]$—subgrid viscous stress

• $\left[\overline{P\left(\mathrm{\Psi }\right)}-P\left(\overline{\mathrm{\Psi }}\right)\right]$—subgrid pressure

• $\left[\overline{Q_j\left(\mathrm{\Psi }\right)}-Q_j\left(\overline{\mathrm{\Psi }}\right)\right]$—subgrid heat flux vector

• $\left[\overline{P\left(\mathrm{\Psi }\right)u_j}-P\left(\overline{\mathrm{\Psi }}\right)\widetilde{u_j}\right]$—subgrid pressure flux

• $\left[\overline{u_i{\tau }_{\mathit{ij}}\left(\mathrm{\Psi }\right)}-\widetilde{u_i}\overline{{\tau }_{\mathit{ij}}}\left(\overline{\mathrm{\Psi }}\right)\right]$—subgrid viscous diffusion

• $\left[\overline{\rho }\widetilde{u_j{e}_t}-\overline{\rho }\widetilde{u_j}\widetilde{e_t}\right]$—subgrid convective energy

• $\left[{\sum }_{\gamma }\overline{H_{,\gamma }\left(\mathrm{\Psi }\right)J_{j,\gamma }\left(\mathrm{\Psi }\right)}-{\sum }_{\gamma }{H}_{,\gamma }\left(\overline{\mathrm{\Psi }}\right)J_{j,\gamma }\left(\overline{\mathrm{\Psi }}\right)\right]$—subgrid molecular enthalpy flux

• $\left[\overline{\rho }\widetilde{Y_{\gamma }{u}_j}-\overline{\rho }\widetilde{Y_{\gamma }}\widetilde{u_j}\right]$—subgrid scalar flux

• $\left[\overline{J_{j,\gamma }\left(\mathrm{\Psi }\right)}-J_{j,\gamma }\left(\overline{\mathrm{\Psi }}\right)\right]$—subgrid mass flux vector

Of these, the subgrid turbulent stress, and filtered reaction terms are usually given the most importance in LES modeling. The rest of the chapter discusses some of the modeling techniques in literature that are used to close these terms. The effects of subgrid pressure, subgrid molecular enthalpy flux, the subgrid shear-stress tensor, the subgrid heat and mass flux vectors are typically negligible for flows at atmospheric pressures. However, for reacting flows at elevated pressures the subgrid pressure, heat and mass flux vectors are not negligible as shown in the studies of Ma et al. [14], Foster and Miller [15], due to the presence of large scalar gradients. Closure for these terms are not presented here, but some of the methods described in this chapter may be used to model these terms as well.

2.3 Numerical approach

LES of compressible flows result in generation of acoustic waves and large instabilities in both low and high Mach number regimes. Choice of appropriate numerical technique, low diffusivity schemes for discretization of the equations, appropriate filters, boundary conditions for proper inflow and outflow of acoustic waves, entropy waves, and numerical instabilities with minimal reflection for high quality simulations are critical. Several methods to address the numerical implementations for compressible LES are reported in the studies of Jiang and Lai [16].

3. Filtering techniques

Filtering in LES may be performed either in a physical space or in a spectral space as illustrated in Figure 1. The filter functions are generally required to satisfy a few mathematical properties including normalization, linearity and commutativity [17].

3.1 Physical space filters

The simplest, physical space filter is a grid based filter called top hat / box filter that can be represented as, $G\left(x\right)=\frac{1}{\mathrm{\Delta }}H\left(\frac{1}{2}\mathrm{\Delta }-|x|\right),$ where $H$ is the Heavyside function, and $G$ a filter kernel that is spatially and temporally invariant. This is an implicit filter, where the grid / mesh size is used as the control. The flow characteristics larger than the grid are resolved, while those smaller than the grid require modeling.

3.2 Spectral space/Fourier mode based filters

In spectral space, a cut off wave number (or frequency) can be chosen and a filter kernel at that wave number can applied. Applying this low pass filter to the variable removes the small scales associated with high frequencies. To perform the filtering operation, Fourier transforms of the variable $\left(\widehat{\mathrm{\Phi }}\right)$ and the filter kernel $\left(\widehat{G}\right)$ are first obtained. The filtering is then performed as a convolution operation, $\overline{\widehat{\mathrm{\Phi }}}\left(k,\omega \right)=\widehat{G}\left(k,\omega \right)\ast \widehat{\mathrm{\Phi }}\left(k,\omega \right)$.

3.3 Proper orthogonal decomposition (POD) assisted filters

POD is a method primarily applied to development of reduced order models in CFD [18, 19]. For a given data, POD extracts the modes that contain the most dominant characteristics of the data. POD of a turbulent flow results in scale separation. The POD modes represent all the characteristic scales of motion. In this method, a singular value decomposition (SVD) of the turbulent flow field data $\mathrm{\Phi }\left(x,t\right)$ is performed. This results in a decomposition of the form, $\mathit{SVD}\left(\mathrm{\Phi }\right)={\mathrm{\Theta \Sigma \Xi }}^{\ast }$, where $\mathrm{\Theta }$ represents the spatial modes (orthonormal basis functions) of the flow field data, $\mathrm{\Xi }$ the corresponding temporal coefficients and $\mathrm{\Sigma }$ the magnitude of significance of the modes. The instantaneous value of the field variable can be expressed as,

$$\mathrm{\Phi }\left(x,t\right)=\sum _{n=1}^{\infty }{\theta }_n\left(x\right){\xi }_n\left(t\right)$$E9

It is noted that not all $\infty$ modes are actually required to reconstruct the data. A finite number of modes $\left(n<N_c\right)$ typically contain upto $80\%$ of the energy of the flow field data. This is analogous to other LES filtering methods, where a set filter size captures the most energetic scales of flow. Here, the filter criteria is the number of modes chosen to reconstruct the data. The filtering operation may be defined as follows,

$$\mathrm{\Phi }\left(x,t\right)=\overline{\mathrm{\Phi }}\left(x,t\right)+{\mathrm{\Phi }}^{\prime }\left(x,t\right)=\sum _{n=1}^{N_c}{\theta }_n\left(x\right){\xi }_n\left(t\right)+\sum _{n=N_c+1}^{\infty }{\theta }_n\left(x\right){\xi }_n\left(t\right).$$E10

In this equation, $N_c$ represents the number of modes used to represent the resolved (most energetic) scales of the data. The sub-filter scale modeling for POD assisted filters is discussed in Section 4.1.2.

4. Subgrid scale modeling

The unresolved terms identified in Section 2.2 must be closed using modeling techniques such that forward scatter (energy transfer from large scales to small scales) and possible back scatter (energy transfer from small scales to large scales) are captured. The general approach to modeling is to calculate the subgrid terms as a function of the appropriate resolved variables. This section presents some of the most predominantly used modeling approaches to close the subgrid turbulent stress and the reaction source terms. The general principles outlined here can further be used to model various other unclosed terms discussed in the preceding section.

4.1 Closure for turbulent stress

4.1.1 Eddy viscosity based models

Turbulent stress appears as a result of filtering applied to the momentum equation. While the large scale turbulent features cause mixing and transport the energy, the smaller scales are generally associated with diffusion. Thus, the subgrid turbulent stress is generally modeled similar to viscous diffusion with a special diffusivity factor. This form of model is called an eddy viscosity model or Smagorinsky model. The unclosed turbulent stress is expressed as a function of the resolved strain energy and a modeled diffusivity factor.

$$T_{\mathit{ij}}=-2{\nu }_{\mathit{SGS}}\widetilde{S_{\mathit{ij}}}=C_s{\mathrm{\Delta }}^2\mid \tilde{S}\mid \widetilde{S_{\mathit{ij}}},$$E11

where, ${\nu }_{\mathit{SGS}}$ represents the subgrid scale diffusivity, $\mathrm{\Delta }$ the size of filter, $\widetilde{S_{\mathit{ij}}}=\frac{1}{2}\left(\frac{\partial \widetilde{u_i}}{\partial x_j}+\frac{\partial \widetilde{u_j}}{\partial x_i}\right)$ the resolved strain rate tensor and $C_s$ a model constant typically derived via dimensional analysis. While the form of the model is simple and consistent with the idea of energy cascade, a constant value of $C_s$ makes it too restrictive. Further this model does not predict backscatter of energy. A number of improvements have been proposed for this model over the years. The most popular method is the usage of a dynamic model constant $C_D\left(x,t\right)$ that varies in space and time. The idea here is to perform LES at two different grid (filter) sizes (original grid size ${\mathrm{\Delta }}_0$, and test grid size ${\mathrm{\Delta }}_T\ge 2{\mathrm{\Delta }}_O$). The subgrid scales of the test simulation contain many of the resolved scales of the original simulation. Here it is assumed that the interaction between the turbulent eddies at the length scales between ${\mathrm{\Delta }}_T$ and ${\mathrm{\Delta }}_O$ must be similar to the interaction between the eddies at ${\mathrm{\Delta }}_O$ and its subgrid scales. It is assumed that the same eddy viscosity type model is applicable in both levels. While the subgrid turbulent stress terms $\left(T_{{\mathit{ij}}_O},T_{{\mathit{ij}}_T}\right)$ at each of these filter levels cannot be directly computed, a relationship exists between the two. This is expressed by the Germano identity [20] as, $L_{\mathit{ij}}=T_{{\mathit{ij}}_T}-\tilde{T_{{\mathit{ij}}_0}}=\widetilde{\tilde{u_i}\tilde{u_j}}-\widetilde{\tilde{u_i}}\widetilde{\tilde{u_j}}$. The term $L_{\mathit{ij}}$ is called Leonard’s stress tensor that represents the subgrid turbulent stress at length scales between test and the original grid filters. This can be explicitly calculated from the resolved velocities as a function of space-time and can be modeled similar to the Smagorinsky model as,

$$L_{\mathit{ij}}=C_D\left({\mathrm{\Delta }}_T^2|\tilde{\tilde{S}}|\tilde{\widetilde{S_{\mathit{ij}}}}-{\mathrm{\Delta }}_O^2\tilde{\mid \tilde{S}\mid \widetilde{S_{\mathit{ij}}}}\right).$$E12

The dynamic model constant $C_D$ is then obtained via least square minimization. This model is self contained and requires no prior knowledge of the model constant. This form of model is still predominantly used in LES of compressible real gas flows. However, the variable values of eddy viscosity may cause numerical instabilities and care must be taken when applying this model to highly localized flows.

4.1.2 Scale similarity based models

Scale similarity model proposed by Bardina et al. [21] assume similarities in characteristics between the eddies that are slightly larger than the grid size (smallest resolved scale) and the eddies are slightly smaller than the grid size (largest unresolved scale). With this hypothesis, the subgrid turbulent stress is expressed in terms of the smallest resolved velocities. This approach does not use an eddy viscosity assumption. This approach is also used to model the sub-fliter terms generated by POD based filtering presented in Section 3.3. Applying the scale similarity assumption between the smallest resolved mode ($n=N_c$) and the largest subfilter POD modes ($n>N_c+1$), the drain of modal energy at the cut-off mode can be expressed in terms of a ratio of the singular values $\sigma$ (magnitude of significance of the modes), ${\beta }_i={\sigma }_{N_C}/{\sigma }_{N_C-1}$, which is assumed to be constant for all the sub-filter scales. The energy distribution of higher energy modes can be approximated in terms of coefficients ${\alpha }_i={\left(1+{\sum }_{n=N_c+1}^{\infty }{\beta }_i^{n-N_c}\right)}^{1/2}$. The subgrid turbulent stress is then computed as a function of the energy coefficient ${\alpha }_i$ and the resolved modes $\theta \left(x\right),\xi \left(t\right)$. A complete derivation of the turbulent stress for POD assisted filter is presented in the works of Shinde [22].

4.2 Closure for chemical source terms

Reaction source terms are highly non-linear and computationally stiff. Some of the popular models for closure of the reaction source terms in LES include flamelet libraries and linear eddy models. Each has its own advantages and limitations. However, both are approximations. Major improvements in closure of this term were made after introduction of the Probability Density Function (PDF) method by Pope [23] and Filtered Mass Density Function (FMDF) method by Jaberi et al. [24]. In these methods, the complete statistical information of the flow is provided by a filtered joint PDF for velocity and scalars of the flow. The moments of statistics of the flow are obtained from the solutions of a transport equation derived for the PDF. The primary advantage is that the chemical source term is closed and does not require any modeling. Thus, for combustion at high pressures this approach may generate more accurate results. For a scalar $\varphi \left(x,t\right)$, the filtered mass density function $F_L$ can be defined as,

$$F_L\left(\psi ;x,t\right)={\int }_{-\infty }^{\infty }\rho \left(x^{\prime },t\right)\zeta \left[\psi ,\varphi \left(x^{\prime },t\right)\right]G\left(x^{\prime }\right)d{x}^{\prime },$$E13

where, $\zeta \left[\psi ,\varphi \left(x,t\right)\right]=\delta \left[\psi -\varphi \left(x,t\right)\right]={\mathrm{\Pi }}_{\alpha =1}^{\sigma }\delta \left[\psi -\varphi \left(x,t\right)\right]$ is the fine-grained density with $\delta$ the delta function and $\psi$ the composition domain. The integral property of FMDF results in filtered density, ${\int }_{-\infty }^{\infty }{F}_L\left(\psi ;x,t\right)\mathit{d\psi }={\int }_{-\infty }^{\infty }\rho \left(x^{\prime },t\right)G\left(x^{\prime }\right)d{x}^{\prime }=\overline{\rho }\left(x,t\right)$. From the properties of FMDF, the filtered value of a scalar is obtained by integration over the composition space, as $\widetilde{Q\left(x_i,t\right)}={\int }_{-\infty }^{\infty }\rho \left(x_{i^{\prime }},t\right)Q\left(x_{i^{\prime }},t\right)F_L\left(\psi ;x_i,t\right)\mathit{d\psi }$. Any term that is a function of the variables within FMDF can be directly calculated. The transport equation for FMDF can be expressed as,

$$\begin{array}{l}\frac{\partial F_L\left(\psi ;x,t\right)}{\partial t}+\frac{\partial }{\partial x_i}\left[\widetilde{\left(u_i|\psi \right)}\right]F_L\left(\psi ;x,t\right)]=\\ \frac{\partial }{\partial \psi }\left[⟨\frac{1}{\rho }\frac{\partial J_{j,\alpha }}{\partial x_j}|\psi ⟩F_L\left(\psi ;x,t\right)\right]-\frac{\partial }{\partial \psi }\left[S_{\alpha }\left(\psi \right)F_L\left(\psi ;x,t\right)/\overline{\rho }\right].\end{array}$$E14

The chemical reaction source term in the above equation is closed but two new terms appear that remain unclosed. The first unclosed term $\widetilde{u_i\mid \psi }$, is called SGS scalar convection and represents the motion of scalars at the subgrid scale due to the flow velocity. This term may be decomposed as, $\widetilde{u_i\mid \psi }{F}_L=\widetilde{u_i}{F}_L+\left[\widetilde{u_i\mid \psi }-\widetilde{u_i}\right]F_L$, where the second term is the effect of SGS convective flux. This is typically modeled using an eddy viscosity type model. The second unclosed term $⟨\frac{1}{\rho }\frac{\partial J_{j,\alpha }}{\partial x_j}|\psi ⟩$, is the conditional scalar diffusion that represents the effects of molecular transport (diffusion) in physical space and molecular mixing in composition space. The resulting transport equation is then written as,

$$\frac{\partial F_L}{\partial t}+\frac{\partial \widetilde{u_i}{F}_L}{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(\mu +{\mu }_t\right)\frac{\partial F_L/\overline{\rho }}{\partial x_i}\right]+\frac{\partial }{\partial \psi }\left[⟨\frac{1}{\rho }\frac{\partial J_{j,\alpha }}{\partial x_j}|\psi ⟩-S_{\alpha }/\overline{\rho }\right]F_L.$$E15

This equation can be further integrated to obtain the transport equations for the SGS moments (mean, variance). The spatial transport of the FMDF is represented by a diffusion process governed by a stochastic differential equation

$${\mathit{dX}}_i\left(t\right)=\widetilde{u_i}+\left(1/\overline{\rho }\right)\partial \left(\mu +{\mu }_t\right)/\partial x_i+\sqrt{2\left(\mu +{\mu }_t\right)/\overline{\rho }},$$E16

and the composition evolution is governed by an equation,

$$\frac{d{\varphi }_{\alpha }}{\mathit{dt}}={\mathrm{\Theta }}_{\alpha }-S_{\alpha }=⟨\frac{1}{\rho }\frac{\partial J_{j,\alpha }}{\partial x_j}|\psi ⟩-S_{\alpha }.$$E17

This equation is typically solved via Lagrangian Monte Carlo approach. Here, the term ${\mathrm{\Theta }}_{\alpha }$ is the conditional scalar diffusion (CSD) which remains as the last unknown term in this transport eq. A few mixing models that stem from different physical concepts have been proposed for modeling the CSD. The most popular choice is the Interaction by exchange of mean (IEM) method that is linear and causes the composition of a Lagrangian particle within a filter volume to relax to a mean value. However, it is local only in the physical space and not in the composition space. Improvements to this model include interaction by exchange of conditional mean (IECM) and IEM plus mean drift but they also lack the locality in composition space property. Additional models include Mapping Closure [25] that maps the composition space of a single scalar to a Gaussian reference field and then solves for it evolution. This preserves the properties of mean, variance decay and is local in composition space but is applicable only to a single scalar. An improvement to this model is the Euclidean Minimum Spanning Tree (EMST) [26] that forms a minimum spanning tree based on composition of particles in an ensemble. The lengths of the vertices of the tree are dependent on the weights, and age properties of particles and enables treatment of multiple scalars. The mixing is performed for particles lying in the minimum spanning tree. This preserves the locality property but does not satisfy the linearity and independence properties. Evaluation of various mixing models are presented in the studies of Shetty et al. [27], Chen et al. [28] and Padmanabhan [29]. Most molecular mixing models however, are developed for pure mixing flows. Their application to reacting flows is based on an assumption that the reaction and diffusion are indirectly coupled as the molecular diffusion smooths toward the scalar means in state space. This assumption is not true for reacting flows and the models that worked for pure mixing flows may not work for reacting flows. Further, a common factor between all the models is the mixing frequency ${\mathrm{\Omega }}_m=f\left(\overline{\rho },{\mathrm{\Delta }}_G,\widetilde{u_i},\mu ,D_m\right)$ that is determined from the simulation parameters that represents the rate of change of composition of particles within a filter volume. The standard form of this term is only applicable to flows with constant diffusivity at atmospheric pressure. For flows at elevated pressures with generalized diffusion models, pressure dependent diffusivity corrections must be included in the mixing frequency term as shown in the study by Padmanabhan and Miller [30].

5. Current research in SGS modeling—Machine learning based models

Application of machine learning (ML) to turbulence modeling and control has been a topic of research since early 2000s [31, 32], but over the last several years major strides have been made in terms of development of new ML models, available libraries and collaborative open source codes that enable development of better turbulence closure models for LES. This section briefly covers a few popular supervised and semi-supervised techniques that have been successfully applied to modeling of SGS terms in LES. The SGS closures currently in literature typically pertain to turbulent stress and chemical source terms, but the general approach can be applied to model any of the unclosed terms identified in Section 2.2. Details of implementations of algorithms are not covered in this chapter, as it is beyond the scope of this work.

5.1 Artificial neural network (ANN)/deep learning (DL) based models

An ANN is a network of interconnected nodes in a layered structure (typically, an input, an output and one or more hidden layers). The nodes in the input layer hold numerical values (activations) that represent the features generated from a training dataset and the nodes in the output layer hold values that represent the generated prediction. The nodes in the hidden layers contain activations that map the input to the output. The core of ANN lies in automatic adjustments of these activations in a way that the machine is able to identify a correct mapping between a given input and output, not only for the set of data on which it is trained but also for a new set of data it has never encountered before. The structure of ANN is based on the Perceptron [33], where the activation of each node in a particular layer (except input layer) is a weighted sum of the activations of the nodes in the previous layer. This is expressed as,

$$a_j^L=\sigma \left(\sum _{k=1}^N{w}_{\mathit{jk}}^L{a}_k^{L-1}+b_k^L\right),$$E18

where $N,a_k^L$ represent the number of nodes and activation of node $k$ in layer $L$, $w_{\mathit{jk}}^L$ the weight connecting node $k$ in layer $L$ and node $j$ in layer $L-1$, $b_k^L$ a bias term typically added to rescale this weighted sum to a positive value and $\sigma$ a non-linear activation function (ReLU or Sigmoid function) that is used to force the output between 0 and 1. For CFD simulations, data or statistics computed at various mesh locations within the domain may be used as the input features. The activations of the input nodes are propagated forward using Eq. 18 until the last layer is reached. This is called forward propagation. A loss function is then defined to evaluate the error between the prediction in the output layer and the actual output $y$ determined from a training dataset,

$$\mathcal{L}=\frac{1}{N}\sum _k^N{\left(a_k^L-y_k\right)}^2$$E19

The weights and biases of the nodes in all the layers are then adjusted (from the last layer to first layer) using an optimization technique like gradient descent such that the loss function is minimized. This process is called back propagation. The updates to the weights and biases are applied as, $V=V-{\alpha }_R\mathrm{\partial ℒ}/\partial V$, where $V=\left[w,b\right]$ denotes the weight, bias terms, ${\alpha }_R$ denotes the rate of descent and the derivatives $\mathrm{\partial ℒ}/\partial V$ are determined using the functional dependencies between each layer and chain rule. This form of network can be trained to identify mapping between any input / output pairs. Identifying and generating the right set of inputs / outputs to train a model for optimal performance is the biggest challenge. Modeling approach for SGS terms using ANN can be classified into two categories:

1. Direct modeling of SGS terms using the resolved variables from LES as inputs. A DNS dataset can be explicitly filtered to compute the resolved and SGS terms. The SGS terms from this filtered DNS (FDNS) form the ground truth. A training dataset can then be created using the resolved FDNS variables, filter parameters and measured SGS terms. The resulting network can create a mapping between the resolved variables and the SGS terms. Implementations similar to this approach are presented in the studies of Gamhara and Hattori [34], Wang et al. [35], Park and Choi [36]. A similar approach for closure of reaction source term was proposed by Ranjan et al. [37] where DL was used to compute reaction rates from the filtered species mass fractions, temperatures and pressures. The training data were generated using a variety of methods including 1D linear eddy mixing (LEM), flamelet models and appropriately sliced DNS datasets.

2. Tuning of the model parameters of existing SGS models. Xie et al. [38], Meng et al. [39] generated coefficients used in dynamic SGS modeling, with the filtered flow variables and their gradients computed from LES as inputs, and the coefficients determined from filtered DNS as the ground truth. An illustration of this approach is presented in Figure 2.

Figure 2.

An illustration of deep learning based model via tuning of existing SGS model is presented. The network inputs are the resolved variables measured from LES and the output is a model constant for an applicable eddy viscosity model.

These two approaches are often very cumbersome requiring large training datasets with equivalent DNS and LES states. Equivalency in filter kernel used in LES and explicit filter performed on DNS is not always guaranteed [40]. Time dependence is not typically captured in the network and training is typically performed with a snapshot data. This is acceptable with statistically stationary dataset but is error prone for instantaneously evolving flows. This approach has not shown to outperform the traditional models in terms of computational costs [41]. Further, the predictions by the models are not universal, but specific to the flow conditions used in the training dataset.

In addition to traditional form of NN, Physics Informed Neural Networks (PINN) can be applied to SGS modeling. In this approach, in addition to the traditional loss function, physics constrained loss functions can be imposed on the initial, boundary conditions and the governing equation. As a result the updated network parameters are constrained to satisfy the underlying governing equations. This enables the NN to generate a robust model with less data, and may enable faster computation for a well defined problem. An implementation of PINN for SGS modeling in lagrangian LES is presented in the study by Tian et al. [42].

Another powerful addition to NN, typically applied in image processing, is the idea of Convolutional Neural Network (CNN). In this approach, a convolutional filter is applied to the input, such that a feature map (activations) can be created. From these activations, certain known features of the data can be extracted. In image processing, this translates to extraction of known patterns within the image. These features are then pooled together and then passed to a fully connected neural network (FCNN), which then categorizes the input. In LES, CNN can be used to extract known flow patterns, SGS statistics or other quantities of interest. An implementation of this approach is presented in the study of Xing and Lapeyre [43] who applied CNN to predict the total flame surface density (unresolved due to thin flamefront) from the resolved progress variable $c$. This was then used to predict SGS flame wrinkling.

For high pressure real gas flows, there is a lack of modeling effort and studies evaluating the applicability of existing SGS models for the unclosed terms. Hence, the approach of direct modeling of SGS terms from resolved LES variables is more applicable to closure of the SGS pressure, molecular enthalpy, heat and mass flux vectors terms. While the concerns of equivalency in filter kernels and non-universality of the developed models exist, this approach does not disregard the effects of these terms. This would however require development of a DNS database for high pressure combustion [7, 44].

5.2 Reinforcement learning (RL)/deep reinforcement learning

Reinforcement learning is a ML technique where a model is not explicitly trained by a training set, instead it learns to improve itself by interacting with a dynamic environment. A brief overview of the method and a few core terminologies are provided in this section. An RL framework consists of an agent within an environment at a state $s$. The agent can take an action $a\left(s\right)$ to receive a reward $R\left(s\right)$. The rewards are cumulative in nature. The environment is typically stated as a Markov Decision Process (MDP). For deterministic systems, a transition model $T\left(s^{\prime }|s,a\right)$ can be developed that describes the new state $s^{\prime }$ the agent will be in, given it is at a particular state and takes an action. A policy for the agent can then be defined as a probability ${\pi }_{\theta }\left(s\right)=P\left(a|s\right)$, which determines the actions it must take, given it is in a particular state. Further, a value function can be defined for the agent which quantifies the expected rewards it can receive, for a given policy $V\left(s\right)=\mathbb{E}\left({\sum }_{i}R\left(s_i|s,{\pi }_{\theta }\right)\right)$. The value function therefore is dependent on the path the agent takes. When defining a value function a discount factor $\gamma$ is typically added to the rewards to prioritize the immediate rewards over the time delayed rewards. Policy and value functions can have any functional form including an ANN (in deep reinforcement learning).

The goal of an agent in RL, is ultimately to determine the best policy that maximizes the value function. When a transition function for a system is available, the policy design for an agent is straightforward and can be achieved via value iteration (updating the value function iteratively for every possible sequence to determine the best value) or via policy iteration (for each policy, determine the best value and choose the best policy that maximizes the value). This is called model based learning. For non-deterministic systems, it is not practical to derive a transition function because of the large action space. Instead, an optimal policy can be directly derived from the agent’s interaction with its environment. This approach is called model free learning. One such method is $Q$ learning where a function $Q\left(s,a\right)$ can be learnt for a state-action pair. This function represents the value (quality) of an action at a state. Once a $Q$ function is learnt, the optimal value function and policy can be derived as the ones that return the highest $Q$ value.

However, for continuous action spaces like a CFD simulation, $Q$ learning is not always feasible because a $Q$ value is required for every action and the method becomes unstable. In this case, policy gradient method is more suitable. The premise of policy gradient method is to directly optimize the policy instead of learning a $Q$ function first. One of the most efficient policy gradient methods is proximal policy optimization (PPO). In PPO, the policy and value functions are each represented as an ANN. For each time iteration, the policy network takes states $s$ as input and generates the probabilities $P\left(a|s\right)$. This information is stored as a tuple $\left(s,a,R,P\right)$ in a dataset. The value function network then takes $\left(s,a\right)$ pairs as input and generates $Q\left(s,a\right)$ that quantify the goodness of the decision. Future expected rewards for the decisions at each time step are then computed. The difference between the computed $Q$ values and the total future rewards (for a certain set of actions), generates an advantage function $\widehat{A}\left(s,a\right)=Q\left(s,a\right)-V_{\pi }\left(s\right)$. This advantage function is used to compute a loss, which is then back-propagated through the value function to improve it. A similar loss function is computed using the advantage function $\widehat{A}$ and the stored probabilities $P\left(a|s\right)$, and back-propagated to improve the policy network. A clipped loss function used for the policy network increases the stability of training by constraining the maximum change of policy at each time step. Effectively the policy and value functions are trained and improved together.

An implementation of deep reinforcement learning to model SGS terms is presented in a study by Kurz et al. [45]. In this study, an agent interacts with an environment of LES of homogeneous isotropic turbulence (HIT), and adapts the Smagorinsky model coefficients $C_s\left(x_i,t\right)$ for each element at each time, in the domain. The framework for the model is stated as follows: the agent’s states $s$ are defined by the current flow state of the LES, the agent’s actions space are defined as the element wise Smagorinsky model constants $C_s$, the policy ${\pi }_{\theta }\left(a|s_t\right)$ is defined by a CNN that takes resolved flow variables (momentum field, invariants of velocity) from LES as inputs and predicts applicable $C_s$, and the rewards $R\left(s\right)$ are defined by closeness of the turbulent energy spectra computed for the LES using the predicted $C_s$ to the spectra measured from a filtered DNS solution, $R\left(s\right)=f\left(E_{\mathit{\text{FDNS}}}\left(k\right)-E_{\mathit{LES}}\left(k\right)\right)$. The form of reward used in this case can also be extended to other flow configurations where the energy spectrum follows $K^{-5/3}$ trend. An illustration of this approach is presented in Figure 3. This approach is more efficient than a supervised learning model because the rewards are independent of the filters used. In another study by Novati et al. [46], a multi-agent reinforcement learning (MARL) approach is adopted for SGS turbulent stress tensor modeling. In this work a PPO method is adopted for policy improvement. $N_{\mathit{\text{agents}}}$ are distributed across the computational domain each of which collect the state, action, reward and statistics information. These agents then send the information to a master at the end of each simulation and receive updated policy parameters. The policy network takes invariants of gradients and Hessians of velocity and generates a Smagorinsky constant as output. Two types of rewards are used to control the model. One reward function ensures the sum of resolved and modeled contribution to the SGS stress is independent of the LES resolution (filter size), while the other reward function ensures the statistics between LES and DNS data are comparable. The MARL approach exploits the potential of cooperating agents resulting in an improved computation efficiency.

Figure 3.

At each time step, the agent computes $C_s$ and passes to the environment. The environment then computes the SGS terms, the energy spectra and the corresponding reward. This information is then passed to the agent, which then adjusts its actions ($C_s$) based on the rewards it received.

5.3 Super resolution

With the continuous development of newer ML methods and libraries, creative ML based modeling techniques are being applied in the field of turbulence modeling. One such approach is super resolution (SR) that aims at scaling up the data from a low resolution simulation to produce high resolution data. A trivial method to accomplish this would be to perform interpolation of data from the neighbors of a point. However, interpolation typically produces an average of the data which is not indicative of high resolution data. Instead, an ANN can be trained to learn this scaling based on prior information it collected from large sets of data. The ANN performs upscaling in a manner that fills the missing information. This approach does not violate the data processing inequality theory because the missing information is inherently present in the large set of input data. One approach to perform a direct mapping between low resolution LES and high resolution DNS data would be to use a super resolution convolution neural network (SRCNN), which uses a CNN for pattern recognition (flow structure / statistics in LES) and a SR model to fill the information gap. Some implementations of SR for modeling in LES include a Generative Adversarial Netowrk (GAN) based model by Nista et al. [47] that has a generator network to upscale the LES information and a discriminator to evaluate the quality of scaling, a physics guided super resolution network (PGSRN) by Chen et al. [48] that combines the idea of PINN and SR, and a meta learning deep CNN by Zhao et al. [49] that generates high resolution data without manually inputting filtered DNS information.

6. Future research challenges and concluding remarks

Research and progress in ML is seeing an unparalleled growth because of its applications in every field of science, technology, and business. While ML has tremendous potential, research focus must also be directed toward incorporation of fundamental physics and mathematical attributes of turbulence in the models. This is especially true for problems like combustion because of the complex flow physics and two-way coupling between turbulence and flames. Further, the vast majority of modeling efforts in literature only cover the SGS turbulent stress and chemical reaction terms. However, for flows at high pressures, additional unclosed SGS terms are observed to be significant and future modeling efforts in LES of real gas flows must take these terms into account. Many of the ML based models presented in this chapter use simple flow configurations and are intended to be a proof of concept. Application of these models to more complex flow scenarios must be evaluated. The ultimate goal of LES modeling is to develop accurate models with high computational efficiency. Thus, along with exploration of new innovative methods, improving the computational efficiency of the current methods must also be explored.