Implicit Large Eddy Simulations (ILES)

• ILES is a turbulence modeling method where discretization intrinsically filters unresolved scales, eliminating the need for explicit subgrid-scale models.
• ILES relies on spectral selectivity, using upwinding and stabilization operators to minimally affect large scales while damping only high-wavenumber modes.
• ILES implementations in high-order DG, flux reconstruction, and stabilized finite elements have been shown to accurately reproduce energy budgets and turbulent flow characteristics.

Implicit Large Eddy Simulations (ILES) denote a class of large-eddy methodologies in which the effect of unresolved scales is supplied by the discretization itself rather than by an explicit subgrid-scale closure. In this sense, ILES resolves the energy-containing motions while relying on scheme-dependent numerical dissipation—arising from upwinding, monotone fluxes, flux reconstruction filters, stabilization operators, or related mechanisms—to remove energy or enstrophy near the grid cutoff. Across high-order discontinuous Galerkin, flux reconstruction, stabilized finite elements, spectral vanishing viscosity, kinetic Boltzmann formulations, and astrophysical shock-capturing schemes, the central requirement is not merely dissipation, but dissipation with the correct spectral selectivity: resolved large scales should remain minimally contaminated while the highest resolvable modes are damped in a controlled way (Ata et al., 2 Jan 2026).

1. Definition and conceptual scope

ILES is distinguished from explicit LES by the absence of an added SGS model in the governing equations. In explicit LES, filtered equations are advanced together with a closure such as an eddy viscosity, for example Smagorinsky-type, dynamic, WALE, or related models; in ILES, the discrete approximation of transport terms itself provides the regularization associated with unresolved scales (Maulik et al., 2019). In the formulation quoted in the wingtip-vortex study, and attributed there to Sagaut, ILES assumes that “the action of subgrid scales on the resolved scales is equivalent to a strictly dissipative action” (Lombard et al., 2015).

The filtered-equation perspective remains standard. For incompressible flow, the SGS stress can be written as

$\tau_{ij}=\overline{u_i u_j}-\bar{u}_i \bar{u}_j,$

whereas in compressible settings Favre filtering and its associated SGS stresses and fluxes provide the corresponding decomposition (Yuan et al., 2020). ILES does not eliminate these unresolved interactions physically; rather, it represents their net effect through the discretization. This suggests that the decisive question is not whether a scheme is formally “implicit” or “explicit,” but whether its numerical transfer of energy across the cutoff reproduces acceptable cascade behavior.

The modern literature represented here spans several distinct realizations of that idea. High-order DG-Boltzmann ILES uses a nodal discontinuous Galerkin discretization of the Boltzmann equation with BGK relaxation to recover nearly incompressible Navier–Stokes behavior while relying on interelement jump dissipation and physical relaxation instead of explicit SGS terms (Ata et al., 2 Jan 2026). Nonlinearly stable flux reconstruction treats truncation error, face upwinding, and correction operators as an implicit filter for subsonic viscous turbulence (Brillon et al., 2024). Stabilized finite elements interpret residual-based SUPG/PSPG terms as the implicit SGS mechanism in separated incompressible flow (Saavedra et al., 2022). Spectral vanishing viscosity supplies high-wavenumber dissipation in spectral/$hp$ element ILES of vortex-dominated flows (Lombard et al., 2015). In astrophysical and cosmological practice, ILES often denotes shock-capturing finite-volume simulations in which unresolved turbulence is dissipated by numerical diffusion alone (Schmidt-Brückner, 8 Sep 2025).

2. Spectral and energetic principles

A recurrent principle across the cited work is that successful ILES requires narrow-band dissipation. The DG-Boltzmann study states this explicitly: dissipation should be concentrated in a narrow band of high wavenumbers so that the energy-containing and inertial-range scales remain minimally affected; otherwise nonlinear energy transfer is suppressed, the $k^{-5/3}$ spectrum is distorted, and cascade physics is undermined (Ata et al., 2 Jan 2026). High-order DG is attractive precisely because compact stencils, high spectral resolution, and suitable numerical fluxes can keep dispersion and dissipation low over most of the spectrum and allow damping to rise only near the cutoff.

This spectral selectivity is commonly diagnosed through energy budgets. For the Taylor–Green vortex in the DG-Boltzmann formulation, the kinetic energy, enstrophy, and dissipation measures are

$$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$

$$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$

In that setting, high-order DG on a moderate grid produced small $\epsilon_{\mathrm{num}}$ at $N=5$–6 and spectra close to the inertial-range $k^{-5/3}$ slope, while the local Lax–Friedrichs flux increased numerical dissipation and was therefore disfavored for ILES (Ata et al., 2 Jan 2026).

Energy-stability arguments make the same point from the discrete side. In the DG-Boltzmann formulation, upwind numerical fluxes generate a non-positive semidefinite quadratic form in interelement jumps, so the linear DG operator is energy stable and dissipates only through those jumps; spectrally, this damps the highest resolvable wavenumbers (Ata et al., 2 Jan 2026). In nonlinearly stable flux reconstruction, the correction operator $K(c)$ modifies the mass matrix and acts as a linear filter on the DG residual, while face upwinding provides targeted dissipation; the resulting scheme realizes a fully implicit LES without an explicit eddy viscosity (Brillon et al., 2024).

The consequences of insufficient selectivity can be severe. In weakly compressible, anisotropic turbulence relevant to core-collapse supernovae, ILES based on high-resolution shock-capturing methods was found to capture the large-scale dissipation rate reasonably well at modest resolution, yet intermediate scales were dominated by a strong bottleneck, and an inertial range was not recovered until roughly $512^3$ resolution for the least dissipative schemes (Radice et al., 2015). This is not a contradiction of ILES as such; it is evidence that the numerical SGS action can be acceptable at large scales while remaining overly intrusive at intermediate ones.

3. Numerical realizations

The numerical mechanisms that realize ILES vary substantially, but several families recur.

High-order DG-Boltzmann ILES solves the continuous Boltzmann equation with a BGK collision operator,

$hp$0

and reduces velocity space through a second-order trivariate Hermite–Galerkin expansion to a 10-component first-order PDE system. The first four moments carry mass and momentum, the last six the deviatoric stress, and the nearly incompressible limit recovers Navier–Stokes with $hp$1. Spatial discretization uses nodal DG on hexahedra with Gauss–Legendre–Lobatto points; upwind fluxes provide interface dissipation, the nonlinear BGK term is over-integrated to suppress aliasing, and a third-order semi-analytic Runge–Kutta scheme integrates the stiff relaxation exactly while preserving the advective CFL limit (Ata et al., 2 Jan 2026).

Flux reconstruction and correction-procedure approaches implement ILES differently. In the nonlinearly stable flux-reconstruction work, entropy-conservative two-point volume fluxes are combined with optional Roe-type face upwinding, and the correction parameter $hp$2 damps the highest modal residual content while allowing larger stable explicit time steps. The paper reports that explicit eddy-viscosity SGS models did not improve Taylor–Green or decaying homogeneous isotropic turbulence in the tested subsonic viscous free-shear flows; NSFR’s built-in dissipation was sufficient, though care with face terms and oversampling was required for reliable pressure-dilatation diagnostics and spectra (Brillon et al., 2024). A separate FR study for turbomachinery combined HLLC or Roe internal-face fluxes with a local modal filter and showed that HLLC was more robust, while an entropy-based positivity-preserving filter stabilized the solution at the cost of greater dissipation (Wang, 2024). A dynamically load-balanced $hp$3-adaptive FR method extended the same ILES logic by using Roe fluxes, BR2 viscous treatment, ESDIRK time integration, and spectral-decay smoothness indicators to refine or coarsen order in under-resolved transitional flows (Wang et al., 2019).

Projection-based upwind formulations reinterpret ILES as a variational multiscale selection of which jump content should be damped. In the mass-conserving mixed stress framework, exact pointwise divergence-free velocity fields are constructed in $hp$4-conforming spaces, so the baseline discretization is intrinsically low-dissipative. Standard DG upwind stabilization then acts as the implicit SGS model, while the high-order projection-based upwind method removes the low-order component of tangential velocity jumps from the penalty and applies dissipation only to the high-order remainder. In turbulent channel flow and periodic hills, this reduced excessive near-wall damping and improved wall-bounded predictions relative to standard ILES (Lederer et al., 2023). The later HOPU formulation applied the same principle to Eppler 387 airfoil flows and showed improved turbulent boundary-layer statistics without materially changing pressure-driven quantities (Lederer et al., 2024).

Other realizations are more explicitly regularizing. Spectral vanishing viscosity applies

$hp$5

with a modal kernel $hp$6 that vanishes below a cutoff $hp$7 and damps only higher modes. In the NACA 0012 wingtip vortex study, this SVV-iLES approach stabilized an under-resolved spectral/$hp$8 simulation at $hp$9 while preserving low-wavenumber accuracy (Lombard et al., 2015). Stabilized finite elements, by contrast, realize ILES through residual-based SUPG/PSPG terms whose stabilization parameter depends on $k^{-5/3}$0, convection, and diffusion scales; in periodic hills, this numerical dissipation played the SGS role and interacted strongly with mesh and time step (Saavedra et al., 2022). Interior Embedded DG provided yet another variant, where Roe or Lax–Friedrichs stabilization on face traces, together with a hybridizable first-order formulation and Newton–GMRES solution strategy, supported high-order ILES of transitional compressor-cascade flow (Fernandez et al., 2016).

4. Diagnostics, benchmarks, and representative flows

ILES performance is assessed less by formal closure coefficients than by energy budgets, spectra, transition markers, wall statistics, and integral forces. The Taylor–Green vortex is especially prominent because it exposes whether a method reproduces correct energy decay, dissipation peaks, and inertial-range spectra. In the DG-Boltzmann study, the method reproduced the peak dissipation near $k^{-5/3}$1 at $k^{-5/3}$2, and high-order runs on $k^{-5/3}$3 elements with $k^{-5/3}$4 agreed closely with a $k^{-5/3}$5 pseudospectral reference (Ata et al., 2 Jan 2026). NSFR showed similar behavior in Taylor–Green, with oversampling of the velocity field to $k^{-5/3}$6 equidistant points per direction required to avoid a spurious pile-up of turbulent kinetic energy at the smallest resolved scales (Brillon et al., 2024).

External aerodynamic and transitional benchmarks test a different aspect of ILES: its ability to capture laminar separation, natural transition, vortex shedding, and near-wall stresses without explicit SGS support. The high-order IEDG method captured laminar separation, TS-wave amplification, KH instability, and turbulent reattachment over a NACA 65-(18)10 compressor cascade at $k^{-5/3}$7 and $k^{-5/3}$8, while using significantly fewer degrees of freedom than a second-order finite-volume LES with an explicit SGS model (Fernandez et al., 2016). The high-order projection-based upwind method reproduced mean pressure coefficient distributions, lift and drag, and separation/reattachment locations for Eppler 387 airfoil flows at $k^{-5/3}$9 and $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$0, while improving $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$1, $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$2, and $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$3 in turbulent boundary layers (Lederer et al., 2024). The wingtip-vortex SVV-iLES captured jetting, core pressure, and secondary-vortex interactions for a NACA 0012 wingtip at $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$4 (Lombard et al., 2015).

Separated and wall-bounded turbulence provide further tests of the coupling between numerical stabilization and physical predictions. In periodic hills at $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$5, a stabilized continuous finite element ILES produced good mean velocities and Reynolds stresses with fewer degrees of freedom than the reference numerical solution, but its accuracy depended on the interaction between time step, mesh size, and the stabilization parameter $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$6; coarse meshes were especially sensitive (Saavedra et al., 2022). For the Imperial Front Wing section in ground effect, Kennedy–Gruber entropy-stable DG ILES captured transition and laminar separation bubbles more effectively than an explicit Vreman LES baseline while permitting a larger stable time step at similar cost per iteration (Ntoukas et al., 2024). In a flow over a sphere at $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$7, DG-Boltzmann ILES predicted mean drag coefficients $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$8 and $$E_k(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\mathbf{u}\cdot\mathbf{u}\,d\Omega,\qquad \zeta(t)=\frac{1}{\rho_\infty\,\Omega}\int_\Omega \frac{1}{2}\rho\,\boldsymbol{\omega}\cdot\boldsymbol{\omega}\,d\Omega,$$9, Strouhal numbers $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$0 and $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$1, and the expected laminar separation, transition, and helical-like shedding structure on a $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$2-element mesh with $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$3 and $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$4 (Ata et al., 2 Jan 2026).

5. Modern reinterpretations and application domains

Recent work increasingly treats ILES not as a monolithic alternative to SGS modeling but as a numerically encoded filter whose action can be analyzed, guided, or combined with data-driven decisions. In two-dimensional Kraichnan turbulence, a machine-learning classifier was trained on DNS-derived a priori eddy viscosities and used to switch between an energy- and enstrophy-preserving Arakawa scheme and an upwind-biased ILES scheme. The result was a selective ILES that applied dissipation only where the DNS-informed classifier indicated positive eddy viscosity, thereby improving angle-averaged kinetic energy spectra and vorticity structure functions relative to uniform ILES or uniform Arakawa discretization (Maulik et al., 2019). A related spectral-difference study learned an element-local, time-local modal filter that maps DNS to ILES behavior. For short restart windows in Taylor–Green, the learned filter showed that higher polynomial orders are less dissipative, with equivalent filter widths $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$5 for $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$6–5, $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$7 for $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$8, and $$\epsilon_1(t)=-\frac{dE_k}{dt},\qquad \epsilon_2(t)=\frac{2\mu}{\rho_\infty}\,\zeta(t),\qquad \epsilon_{\mathrm{num}}=\epsilon_1-\epsilon_2.$$9 for $\epsilon_{\mathrm{num}}$0, and that the implicit filter supports both forward transfer and appreciable local backscatter (Clinco et al., 2024).

Application domains are correspondingly broad. In turbomachinery, high-order FR ILES was validated on the T106c low-pressure turbine and VKI LS89 high-pressure turbine, with the mean heat-transfer coefficient on LS89 showing that $\epsilon_{\mathrm{num}}$1–60 is suitable for a satisfactory prediction of the mean flow using a relatively coarse spanwise spacing (Wang, 2024). In astrophysics, ILES is often the default because shock-capturing methods are simple and robust for highly compressible flows; however, explicit SGS models remain valuable where unresolved turbulent energy, diffusion, or dynamo action must be represented consistently (Schmidt-Brückner, 8 Sep 2025). Cosmological simulations with adaptive mesh refinement emphasize that, unlike explicit SGS LES, pure ILES has no explicit unresolved kinetic-energy reservoir and therefore cannot enforce global conservation of resolved-plus-unresolved turbulent energy across refinement changes (Schmidt et al., 2013). Global solar-convection simulations with EULAG-MHD show that ILES can achieve converged non-rotating large-scale convection with sufficient radial resolution, yet in rotating cases the effective viscosity changes with resolution strongly enough to alter the differential-rotation regime from solar-like to anti-solar (Guerrero et al., 2022).

6. Limitations, controversies, and practical guidance

One common misconception is that any stable high-order method automatically furnishes a good ILES. The collected studies do not support that view. Numerical damping must be spectrally selective, and diagnostics must verify that energy decay, dissipation, spectra, and coherent structures are not being distorted. Upwind fluxes may be suitable, but more dissipative alternatives such as LLF can be counterproductive in under-resolved turbulence (Ata et al., 2 Jan 2026). Oversampling or dealiasing may be essential for trustworthy spectra and energy budgets, particularly in split-form or collocated high-order schemes (Brillon et al., 2024).

A second misconception is that explicit SGS models are always inferior to ILES. The evidence is mixed by design. In NSFR free-shear turbulence, explicit eddy-viscosity models did not improve results and often overdamped the smallest resolved scales (Brillon et al., 2024). In forced homogeneous isotropic turbulence, by contrast, deconvolutional neural-network SGS models outperformed ILES, dynamic Smagorinsky, and dynamic mixed models in energy spectra, SGS-flux PDFs, structure functions, and coherent structures, while the ILES baseline overestimated small-scale energy (Yuan et al., 2020). This suggests that the adequacy of ILES is problem-dependent and tied to the numerical filter’s match to the actual unresolved transfer.

A third controversy concerns physical admissibility. One line of argument, developed for Euler-equation turbulence, contends that ILES based on minimal subgrid removal violates a maximal entropy production admissibility criterion and therefore underpredicts dissipation and intermittency relative to dynamic SGS LES (Glimm et al., 2018). The astrophysical review literature presents a more moderate conclusion: because most astrophysical solvers already possess significant numerical diffusion, the additional effect of SGS models is often small, but explicit modeling can improve convergence for processes such as magnetic-field amplification in binary neutron star mergers or turbulent metal diffusion in mesh-free codes (Schmidt-Brückner, 8 Sep 2025). These positions are not identical, but they converge on one point: ILES is not a substitute for careful budget analysis.

Practical guidance in the cited work is correspondingly technical. For high-order DG-Boltzmann ILES, choose $\epsilon_{\mathrm{num}}$2, set $\epsilon_{\mathrm{num}}$3, use upwind rather than LLF fluxes, over-integrate the nonlinear collision term, and prefer higher polynomial order before drastic mesh refinement (Ata et al., 2 Jan 2026). For NSFR, a positive correction parameter $\epsilon_{\mathrm{num}}$4 can enlarge the stable CFL, Roe face upwinding suppresses spurious pressure-dilatation oscillations, and spectra should be computed from oversampled velocity fields (Brillon et al., 2024). For projection-based upwind wall-bounded ILES, adaptive damping of only the high-order facet-jump content is beneficial in highly under-resolved channel and hill flows (Lederer et al., 2023). In stabilized finite-element ILES, the time step is itself part of the implicit SGS mechanism because it enters the stabilization parameter $\epsilon_{\mathrm{num}}$5; on coarse meshes, reducing $\epsilon_{\mathrm{num}}$6 can reduce the dissipation needed for acceptable results (Saavedra et al., 2022).

Taken together, these studies define ILES less as a single method than as a design principle: unresolved turbulence is represented through numerically induced transfer and dissipation, but only successful if that transfer is narrow-band, energetically credible, and empirically validated against the flow physics of interest.