Open AccessReview

Bridging Large Eddy Simulation and Reduced-Order Modeling of Convection-Dominated Flows through Spatial Filtering: Review and Perspectives

Department of Mathematics, University of Houston, 3551 Cullen Blvd, Houston, TX 77204, USA

Department of Mechanical, Aerospace and Biomedical Engineering, University of Tennessee, 1512 Middle Drive, Knoxville, TN 37996, USA

Department of Mathematics, Emory University, 400 Dowman Drive, Atlanta, GA 30322, USA

⁴

Department of Computer Science, Emory University, 400 Dowman Drive, Atlanta, GA 30322, USA

⁵

Department of Mathematics, Virginia Tech, 225 Stanger Street, Blacksburg, VA 24061, USA

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Fluids 2024, 9(8), 178; https://doi.org/10.3390/fluids9080178

Submission received: 29 June 2024 / Revised: 29 July 2024 / Accepted: 31 July 2024 / Published: 4 August 2024

(This article belongs to the Special Issue Recent Advances in Fluid Mechanics: Feature Papers, 2024)

Download

Browse Figures

Versions Notes

Abstract

Reduced-order models (ROMs) have achieved a lot of success in reducing the computational cost of traditional numerical methods across many disciplines. In fluid dynamics, ROMs have been successful in providing efficient and relatively accurate solutions for the numerical simulation of laminar flows. For convection-dominated (e.g., turbulent) flows, however, standard ROMs generally yield inaccurate results, usually affected by spurious oscillations. Thus, ROMs are usually equipped with numerical stabilization or closure models in order to account for the effect of the discarded modes. The literature on ROM closures and stabilizations is large and growing fast. In this paper, instead of reviewing all the ROM closures and stabilizations, we took a more modest step and focused on one particular type of ROM closure and stabilization that is inspired by large eddy simulation (LES), a classical strategy in computational fluid dynamics (CFD). These ROMs, which we call LES-ROMs, are extremely easy to implement, very efficient, and accurate. Indeed, LES-ROMs are modular and generally require minimal modifications to standard (“legacy”) ROM formulations. Furthermore, the computational overhead of these modifications is minimal. Finally, carefully tuned LES-ROMs can accurately capture the average physical quantities of interest in challenging convection-dominated flows in science and engineering applications. LES-ROMs are constructed by leveraging spatial filtering, which is the same principle used to build classical LES models. This ensures a modeling consistency between LES-ROMs and the approaches that generated the data used to train them. It also “bridges” two distinct research fields (LES and ROMs) that have been disconnected until now. This paper is a review of LES-ROMs, with a particular focus on the LES concepts and models that enable the construction of LES-inspired ROMs and the bridging of LES and reduced-order modeling. This paper starts with a description of a versatile LES strategy called evolve–filter–relax (EFR) that has been successfully used as a full-order method for both incompressible and compressible convection-dominated flows. We present evidence of this success. We then show how the EFR strategy, and spatial filtering in general, can be leveraged to construct LES-ROMs (e.g., EFR-ROM). Several applications of LES-ROMs to the numerical simulation of incompressible and compressible convection-dominated flows are presented. Finally, we draw conclusions and outline several research directions and open questions in LES-ROM development. While we do not claim this review to be comprehensive, we certainly hope it serves as a brief and friendly introduction to this exciting research area, which we believe has a lot of potential in the practical numerical simulation of convection-dominated flows in science, engineering, and medicine.

Keywords:

large eddy simulation; reduced-order modeling; spatial filtering; machine learning; incompressible fluids; compressible fluids; cardiovascular modeling; atmospheric modeling

1. Introduction

In a couple of seminal papers published in 1941 [1,2], Kolmogorov showed that the vortices and eddies in a fluid flow span an increasingly large range of scales as the inertial forces become dominant over the viscous forces. The aim of a direct numerical simulation (DNS) is to capture all these flow structures with a computational mesh sufficiently refined to resolve even the smallest scales. For a large number of important applications in science, engineering, and medicine, a DNS requires extremely fine meshes, leading to exorbitant computational costs for current computing technology. Examples of such applications include atmospheric flow for weather and climate predictions, flow around wind turbines for efficient energy production, and blood flow in larger arteries to predict the progression of a disease or plan for surgery.

The goal to reduce the computational cost of DNS while maintaining a reasonable level of accuracy has motivated a vast body of literature. Among the many proposed approaches, Reynolds averaged Navier–Stokes (RANS) and large eddy simulation (LES) have become popular and widespread across different sciences and industries. In broad terms, RANS models average the Navier–Stokes equations in different ways (quite often in time), while LES techniques filter them (usually in space). This considerably lowers the computational cost but requires closures to represent certain terms emerging from the averaging/filtering process. For closure, RANS models typically rely on empirical descriptions calibrated with available data. While one important line of research seeks to improve the accuracy of RANS models by systematically informing them with larger datasets [3], this review focuses on LES, which is a mathematically and physically justified approach [4,5,6].

LES techniques directly resolve a portion of the flow scales and require a model to account for the remaining (small) unresolved scales that are not directly captured due to mesh under-refinement. Much of the literature has been dedicated to the derivation and performance analysis of LES models, which inevitably feature parameters. The combination of properly tuned LES models and high mesh resolution has been shown to yield accurate computational results in fields ranging from blood flow studies [7,8,9] to weather and climate forecasts [10,11,12,13,14,15,16,17,18,19]. However, while feasible (unlike DNS), these highly accurate LES simulations still carry a high computational cost. These hefty computational costs could be acceptable in exploratory studies, but they are incompatible with pressing deadlines and/or multiquery contexts.

Multiquery contexts arise from needing to, e.g., assess the uncertainty in the computed solution, optimally control a system, or solve an inverse problem for parameter identification or optimization. They are as ubiquitous as tight deadlines in science, engineering, and medicine. Unable to be patient for the long computational times of accurate simulations, LES practitioners typically trade accuracy (high-resolution meshes and costly parameter tuning) for computational efficiency (coarse meshes and parameters set empirically). This trade-off limits the full potential of computational simulations and has possible detrimental consequences on the reliability of a prediction (e.g., weather and climate), the structural integrity of a physical asset (e.g., a wind turbine), or the life of an individual (e.g., surgical planning for cardiovascular disease). Drastically reducing the computational cost of highly accurate LES is a significant research goal that would enable tremendous advancements across sciences and industries.

A viable way to rectify the imbalance between computational efficiency and accuracy is by bridging two fields that have so far remained disconnected: LES and advanced reduced-order modeling (ROM) techniques. ROMs are methods specifically designed to speed up computations in multiquery settings using state-of-the-art linear algebra techniques and error estimation. The natural connections between these two fields emerge around the mathematical operations of spatial (low-pass) filtering and approximate deconvolution, both in the classical sense [4,20,21,22] and reinterpreted through a more modern, data-driven lens [23,24,25,26,27,28]. Filtering can be redefined through the integration of physics-based models and data-driven methodologies, fostering the evolution of hybrid solvers that respect the underlying physics while harnessing the robust structure and pattern recognition capabilities of machine learning. With a research effort on the complex computational mathematics behind spatial filtering and approximate deconvolution for LES and ROM, fast, accurate, and robust simulations of convection-dominated flows will be possible with minimal user set-up. The modular bidirectional interaction between ROM and LES will yield a converged computational twin usable for detailed flow analysis (LES counterpart) and swift decision making (ROM counterpart), dynamically updating itself with incoming measurements and synthetic data.

This review is structured as follows. Section 2 introduces general concepts in LES modeling and then, for clarity of presentation, focuses on one of the most popular and successful models, the EFR strategy, with details on the application of the EFR algorithm to the incompressible Navier–Stokes equations and the weakly compressible Euler equations. Section 3 presents results obtained with the EFR method as a full-order method for complex applications arising in hemodynamics and physics of the atmosphere and discusses related open problems. In Section 4, we show how key concepts in the construction of LES models, i.e., spatial filters and approximate deconvolution, can be leveraged to construct LES-ROMs. An important feature in Section 4 is the discussion of preliminary numerical analysis results obtained for LES-ROMs. Selected achievements and prospective applications of LES-ROMs for incompressible and compressible fluids are reported in Section 5. Section 6 draws conclusions from this review and summarizes the open perspectives.

2. LES as a Full-Order Model

In this section, we briefly outline some of the main ideas, concepts, and tools used in the development of LES models. We focus on spatial filtering and approximate deconvolution since these mathematical concepts have played an important role in connecting (“bridging”) two research fields, LES and reduced-order modeling, that have been treated separately until recently. We also focus on one particular class of LES models, i.e., the evolve–filter–relax strategy, which has also been central in “bridging” LES and reduced-order modeling. Finally, we present these LES concepts and models for both incompressible and compressible flows. We made this choice to demonstrate that these concepts and models are truly universal and are not restricted to particular classes of flows or computational settings. Of course, our LES presentation is inherently brief and focused on spatial filtering, approximate deconvolution, and the evolve–filter–relax strategy. For a deeper and more comprehensive introduction to LES, the reader is referred to, e.g., the research monographs [4,5,6,29,30].

In order to explain the main ideas of LES, let us briefly summarize some key concepts in Kolmogorov’s theory [1,2].

The kinetic energy (KE) of a fluid, which is the kinetic energy associated with eddies in the flow, defined as

\frac{1}{2}

ρ \int_{Ω} {∥ u ∥}^{2} d ω

, where

ρ

is the fluid density,

u

the fluid velocity,

Ω

the spatial domain where the fluid is flowing, and

| | \cdot {| |}_{L^{2}}

denotes the

L^{2}

norm, is injected in the system at large scales (low wave numbers). Since the large-scale eddies are unstable, they break down, transferring the energy to smaller and smaller eddies. Finally, the KE is dissipated by the viscous forces at small scales (high wave numbers). This process is usually referred to as energy cascade. See Figure 1. The scale at which the viscous forces dissipate energy takes the name of the Kolmogorov scale:

η = {(\frac{ν^{3}}{ε})}^{1 / 4},

(1)

where

ν

is the kinematic viscosity of the fluid and

ε

is the time average of the rate at which the energy is dissipated (see, e.g., [31]). In order to formally define

ε

, let

[0, t_{f}]

be the time interval of interest for the fluid flow. Then,

ε = \underset{t_{f} \to \infty}{lim sup} \frac{1}{t_{f} | Ω |} \int_{0}^{t_{f}} {ν | | \nabla u | |}_{L^{2}}^{2} d t,

where

| Ω |

denotes the measure of the domain.

For a flow in a developed turbulent regime, the dissipation rate has to be of the same magnitude as the production rate, which is the rate at which the KE is supplied to the small scales. A common way to express

ε

is in terms of the macro-scale variables associated with the flow. For this, let U and L be characteristic macroscopic speed and length, respectively. Then, we can write

ε \sim U^{3} / L

[32]. If we plug this into (1), we obtain

η \sim R e^{- 3 / 4} L,

(2)

where

R e

is the Reynolds number

R e = \frac{U L}{ν} .

(3)

The Reynolds number, which weights the inertial forces over the viscous forces, is often used as an indication of whether a flow is turbulent or not. Roughly speaking, at “high” Reynolds numbers, inertial forces dominate and flows tend to be turbulent.

Equation (2) explains the difficulty associated with the numerical solution of flows at high Reynolds numbers: in order to correctly capture the dissipated energy, one needs a grid with spacing

h \sim η

. This approach, called DNS, leads in general to a large number of unknowns and prohibitive computational costs as the Reynolds number increases. If one tries to reduce the computational cost simply by making the mesh coarser (i.e., using a mesh with size h, which fails to resolve the Kolmogorov scale), the energy physically dissipated at the small scales is not properly dissipated at the numerical level. The under-diffusion in the simulation leads to nonphysical computed velocities. In the best-case scenario, nonphysical oscillations in the velocity field eventually make the simulation crash, which points to an obvious problem in the simulation set-up. In the worst-case scenario, the simulation does not crash and one accepts a nonphysical solution as good.

A possible way to reduce the computational cost without sacrificing accuracy is to solve for the flow average using a mesh coarser than that required by DNS and introduce a model for the effects not captured by the flow average alone. Among such alternatives, RANS and LES are the most widely used for practical applications in science and engineering. In the RANS approach, none of the eddies in the flow are directly resolved from the flow equations, meaning that the effect of all eddies, from largest to smallest, is modeled. See Figure 1. In a sense, RANS could be considered the opposite extreme to DNS. RANS models are formally obtained from decomposing instantaneous quantities (velocity and pressure) into a time-averaged component and a fluctuating (over the average) component and applying an averaging procedure to the governing equations. Extra stresses and scalar fluxes emerging from this averaging procedure require closure, i.e., additional equations or assumptions to close the new system of equations. In the LES approach, one resolves only the larger scales by using a “coarse” mesh (again with respect to DNS) and models the effect of the smaller scales that are not directly solved. See Figure 1. LES models, which are a versatile compromise between DNS and RANS, are devised from applying a low-pass filter to the governing equations. Such low-pass filtering can be viewed as a space averaging.

Traditionally, LES models introduce dissipation via momentum fluxes that are linearly dependent upon the rate of strain of the resolved scales. To clarify this, let us consider the incompressible Navier–Stokes equations (NSEs): denoting the velocity by

u

and the pressure by p, the NSEs read

\begin{matrix} ρ \frac{\partial u}{\partial t} + ρ \nabla \cdot (u \otimes u) - \nabla \cdot (- p I + 2 μ ϵ (u)) & = f in Ω \times (0, t_{f}), \end{matrix}

(4)

\begin{matrix} \nabla \cdot u & = 0 in Ω \times (0, t_{f}), \end{matrix}

(5)

where

ϵ (u) = (\nabla u + \nabla u^{T}) / 2

is the strain-rate tensor,

μ

is the dynamic viscosity, and

f

is a forcing term. In real applications, appropriate initial and boundary conditions define the problem to solve, as we will see below. With a traditional LES model, one ultimately solves the following problem: find

\bar{u}

and

\bar{p}

such that

\begin{matrix} ρ \frac{\partial \bar{u}}{\partial t} + ρ \nabla \cdot (\bar{u} \otimes \bar{u}) - \nabla \cdot (- \bar{p} I + 2 (μ + μ_{a}) ϵ (\bar{u})) & = \bar{f} in Ω \times (0, t_{f}), \end{matrix}

(6)

\begin{matrix} \nabla \cdot \bar{u} & = 0 in Ω \times (0, t_{f}), \end{matrix}

(7)

where

\bar{u}

is the filtered velocity,

\bar{p}

is the filtered pressure,

\bar{f}

is the filtered forcing, and

μ_{a}

is an artificial viscosity introduced by the LES model. This is known as eddy viscosity closure [4,5]. Different LES models arise for different definitions of

μ_{a}

One of the most widely used eddy viscosity models is the Smagorinsky model [33]. Introduced in the early 1960s, it defines the artificial viscosity as

\begin{matrix} μ_{a} = ρ {(C_{s} α)}^{2} \sqrt{2 ϵ (\bar{u}) : ϵ (\bar{u})}, C_{s}^{2} = C_{k} \sqrt{\frac{C_{k}}{C_{ϵ}}}, \end{matrix}

(8)

where

α

is a filter width (typically comparable with the mesh size) and

C_{k}

and

C_{ϵ}

are model parameters. The values of such parameters are not universal and need to be tuned for each specific application. The Smagorinsky model has enjoyed much success for several reasons. It is relatively simple and computationally inexpensive compared to other LES models. It is easy to implement and can be used with a wide range of numerical methods, also thanks to the tunable parameters. However, it has one important limitation: the assumption of local balance between the subgrid scale energy production and dissipation. Since such equilibrium conditions do not hold in many practical applications, the Smagorinsky model often results in over-diffusive simulations. A large body of research has been motivated by the aim of improving upon the Smagorinsky model.

Some alternative methods introduce artificial diffusion that can be solution-dependent (see, e.g., [34,35,36,37,38]) or residual-based (see, e.g., [39,40,41,42,43,44]) so that the artificial viscosity vanishes where the solution is smooth and/or decreases as the mesh is refined. Other methods add a set of equations to the discrete governing equations formulated on a coarse mesh. This extra problem can be devised in different ways, for example, by a functional splitting of the solved and unresolved scales as in variational multiscale methods (see, e.g., [45,46,47,48]) or by resorting to the concept of “suitability” of weak solutions [49]. All these methods entail adding one or more extra terms to the governing equations, in the spirit of (6), and have had great accomplishments. However, one practical downside, in this case, is the need to implement ad hoc extensions that are not part of standard solvers, requiring access to legacy source codes (if not rewriting the solver). An attractive option is, therefore, to add the extra problem sequentially to the fluid dynamics model, since in this case the implementation of the LES model would not require any major modification of a legacy solver. This can be achieved by applying a differential nonlinear low-pass filter to the (nonphysical) solution computed with the fluid dynamics model and a “coarse” mesh. Since this particular LES strategy, i.e., nonlinear spatial filtering, has already made an impact in ROMs, below we focus on this approach. This will improve the clarity of our presentation on how to establish a two-way connection between (i.e., “bridging”) LES and ROMs.

2.1. Nonlinear Spatial Filtering for LES

Let us consider a generic fluid model (e.g., the momentum equation of the NSE (4))

\frac{\partial U}{\partial t} + M (U) = 0

(9)

and discretize it in time with, e.g., the backward Euler scheme for simplicity:

\frac{U^{n + 1} - U^{n}}{Δ t} + M (U^{n + 1}) = 0 .

(10)

Here,

U

denotes a generic variable (e.g.,

U = ρ u

in (9) represents the variable in (4)) and

M (\cdot)

denotes nontrivial operators that characterize the fluid model. For simplicity, we do not include any forcing term in (9).

If we further apply space discretization to (10) and use a mesh of size

h ≫ η

, the computed solution typically features nonphysical oscillations. Such oscillations can be attenuated or removed through a filtering operation. This is the idea at the core of a sequential algorithm called evolve–filter–relax (EFR). The EFR algorithm applied to (9) reads as follows:

-: Evolve: Find intermediate variable $V^{n + 1}$ such that

$\begin{matrix} \frac{V^{n + 1} - U^{n}}{Δ t} + M (V^{n + 1}) = 0 . \end{matrix}$

(11)

For this step, one could adopt the same space discretization technique used for (10) and, hence, the same solver.
-: Filter: Find filtered variable ${\bar{V}}^{n + 1}$ such that

$\begin{matrix} {\bar{V}}^{n + 1} = F V^{n + 1}, \end{matrix}$

(12)

where $F$ is a generic spatial filter. Different types of spatial filters will be described throughout this paper (e.g., in (16)). The solver for this step can be added as a module to the solver for (10).
-: Relax: Set

$\begin{matrix} U^{n + 1} = (1 - χ) V^{n + 1} + χ {\bar{V}}^{n + 1}, \end{matrix}$

(13)

with relaxation parameter $χ \in [0, 1]$ . Here, $U^{n + 1}$ is the end-of-step approximation of $U$ at time $t^{n + 1}$ , which “corrects” the intermediate approximation $V^{n + 1}$ .

The appeal of the EFR algorithm lies not only in modularity but also in flexibility, since there exists a large variety of spatial filters that one could use (e.g., Gaussian, box, sharp cutoff) [5,6].

The main idea of using spatial filtering can be illustrated with a schematic like in Figure 2. The challenge is to represent a flow variable, e.g., the velocity,

u

(represented by the maroon curve in Figure 2), on a given coarse mesh of size h. As mentioned in Section 1, having to work with a coarse mesh is generally the case in realistic engineering, geophysical, and, often, medical applications of convection-dominated flows. As illustrated in Figure 2, if

u

has small-scale components of a size below the mesh size, h, then these components will not be accurately represented on the given mesh. More importantly, the inaccuracies in representing

u

on the given coarse mesh are generally propagated and amplified in convection-dominated, nonlinear problems (e.g., the NSEs), yielding spurious numerical oscillations. To eliminate/alleviate these spurious oscillations, the central idea in LES is extremely simple: instead of trying to approximate all the spatial scales in

u

on the given coarse mesh, try to approximate only the large scales in

u

. This can be achieved by using a low pass spatial filter,

F

, to filter

u

. The filtered flow variable,

\bar{u} : = F u

, which is represented by the orange curve in Figure 2, retains only the large spatial structures in

u

. Thus,

\bar{u}

can be represented accurately on the given coarse mesh. Furthermore, the spatial filtering eliminates the small scale, oscillatory component of

u

. This, in a nutshell, is the guiding principle of LES [4,5].

The EFR algorithms, which are also known as filter stabilizations (or regularizations [6]) because they reduce or eliminate unphysical fluctuations in the computed solution, can be categorized as eddy viscosity models for certain discretization choices [50,51]. The connection between the EFR algorithm and LES modeling is easily seen by shifting the index

n + 1

to n in (12) and (13) and plugging them into (11) to obtain

\begin{matrix} \frac{V^{n + 1} - V^{n}}{Δ t} + M (V^{n + 1}) + \frac{χ}{Δ t} (I - F) V^{n} = 0, \end{matrix}

(14)

where

I

is the identity operator. Equation (14) gives us an implicit discretization of problem (9) with the backward Euler scheme and an additional, explicitly treated dissipation term (compare (14) with (10)). Let us assume that

χ = χ_{0} Δ t

, where

χ_{0}

is a time-independent constant. Then, (14) can be seen as a time-stepping scheme for the following problem:

\frac{\partial V}{\partial t} + M (V) + χ_{0} (I - F) V = 0 .

(15)

Thus, the EFR algorithm (11)–(13) can be interpreted as a splitting scheme for problem (15). In turn, this can be read as problem (9) with an additional term designed to be dissipative.

Stabilization based on linear filters [52,53,54] has been widely studied (see, e.g., [55,56,57]). However, the breaking down of eddies into smaller ones until they get damped is a highly nonlinear process. Hence, it seems more appropriate to use nonlinear filters to select the eddies to be damped. A well-established [58,59,60,61,62,63,64] nonlinear filter for (12) is

F = {(I + L)}^{- 1}, L = - \nabla \cdot (δ \nabla), with δ = α^{2} a (V^{n + 1}),

(16)

where

α

can be interpreted as the filtering radius and

a (\cdot) \in (0, 1]

is the so-called indicator function. This is a critical function that determines where and how much eddy viscosity is added, and it is crucial to the success of the EFR method. We discuss possible definitions of

a (\cdot)

in Section 2.1.3. Notice that

F

can be interpreted in view of the theory of maximal monotone operators and their Yosida regularized operators (see [65], Ch. 7). This provides a solid mathematical foundation to the methodology.

The rest of this section focuses on the application of the EFR algorithm to the incompressible NSE (in Section 2.1.1) and the weakly compressible Euler equations (in Section 2.1.2). However, it could be easily applied to other fluid dynamics models, like, e.g, the one-layer [66] and two-layer [67] quasi-geostrophic equations.

2.1.1. The EFR Method for the Incompressible Navier–Stokes Equations

For the sake of completeness, let us consider Equations (4) and (5) supplemented with initial conditions and boundary conditions. Precisely, we are given initial data

u_{0} (x)

such that

u (x, 0) = u_{0} (x), x \in Ω,

and boundary data

u_{D}

and

g

such that

\begin{matrix} u (x, t) & = u_{D} (x, t), x \in Γ_{D}, \end{matrix}

(17)

\begin{matrix} - p n + 2 μ ϵ (u) \cdot n & = g (x, t), x \in Γ_{N}, \end{matrix}

(18)

where

Γ_{D} \cup Γ_{N} = \partial Ω

Γ_{D} \cap Γ_{N} = \emptyset

, and

n

is the outward unit vector normal to

Γ_{N}

. A mix of Dirichlet and Neumann boundary conditions as in (17) and (18) occurs in many practical problems. Other possible boundary conditions (Robin type) are addressed in Remark 1.

Algorithm (11)–(13), with a filter as in (16), applied to problem (4)–(5) reads as follows:

-: Evolve: find intermediate velocity $v^{n + 1}$ and pressure $q^{n + 1}$ such that

$\begin{matrix} \frac{ρ}{Δ t} v^{n + 1} + ρ \nabla \cdot (u^{n} \otimes v^{n + 1}) - 2 μ \nabla \cdot ϵ (v^{n + 1}) + \nabla q^{n + 1} & = f^{n + 1} + \frac{ρ}{Δ t} v^{n}, \end{matrix}$

(19)

$\begin{matrix} \nabla \cdot v^{n + 1} & = 0, \end{matrix}$

(20)

with boundary conditions

$\begin{matrix} v^{n + 1} (x) & = u_{D} (x, t^{n + 1}), x \in Γ_{D}, \end{matrix}$

(21)

$\begin{matrix} - q^{n + 1} n + 2 μ ϵ (v^{n + 1}) \cdot n & = g (x, t^{n + 1}), x \in Γ_{N} . \end{matrix}$

(22)

Again, for simplicity, we used the backward Euler scheme for time discretization. Of course, other time discretization schemes are possible. Notice that, consistent with the time-discretization order, we also linearized the problem with a first-order time extrapolation of the convective velocity. We also set $v^{0} = u_{0}$ .

-: Relax: set

$\begin{matrix} u^{n + 1} & = (1 - χ) v^{n + 1} + χ {\bar{v}}^{n + 1}, \end{matrix}$

(27)

$\begin{matrix} p^{n + 1} & = q^{n + 1} . \end{matrix}$

(28)

This particular version of the EFR method for incompressible flows was first studied in [62].

Note that filter problem (23), (24) can be conveniently rewritten as a generalized Stokes problem. In fact, by multiplying Equation (23) by

ρ / Δ t

and rearranging the terms, we obtain

\begin{matrix} \frac{ρ}{Δ t} {\bar{v}}^{n + 1} - \nabla \cdot (μ_{a} \nabla {\bar{v}}^{n + 1}) + \nabla {\bar{q}}^{n + 1} & = \frac{ρ}{Δ t} v^{n + 1}, μ_{a} = ρ \frac{α^{2}}{Δ t} a (v^{n + 1}), \end{matrix}

(29)

where

{\bar{q}}^{n + 1} = ρ λ^{n + 1} / Δ t

. Problem (29), (24) can be seen as a time-dependent Stokes problem with a nonconstant artificial viscosity

μ_{a}

, discretized by the backward Euler scheme. Unlike

δ

, which has the dimension of a length squared,

μ_{a}

is dimensionally a dynamic viscosity. The big advantage of this rewriting is that a solver for problem (29), (24) can then be obtained by adopting a standard linearized Navier–Stokes solver. Another advantage is that unlike

λ^{n + 1}

{\bar{q}}^{n + 1}

has the same dimensional units as

q^{n + 1}

. Hence, if one replaces (23) and (24) with (29) and (24), it is possible to relax the pressure too. In [59], it was proposed to replace (28) with

p^{n + 1} = (1 - χ) q^{n + 1} + χ {\bar{q}}^{n + 1}

Remark 1.

It can happen that a Robin boundary condition is enforced on part of the domain:

\begin{matrix} - p n + 2 μ ϵ (u) \cdot n + γ u (x, t) = g_{R} (x, t), x \in Γ_{R}, \end{matrix}

where

g_{R}

and γ are given. This kind of boundary condition can arise in certain cases, such as domain decomposition algorithms (see, e.g., [68]) or the use of “transpiration” techniques for fluid–structure interaction problems (see, e.g., [69]). The corresponding boundary condition for the evolve step is

\begin{matrix} - q^{n + 1} n + 2 μ ϵ (v^{n + 1}) \cdot n + γ v^{n + 1} = g_{R} (x, t^{n + 1}), x \in Γ_{R}, \end{matrix}

while in the filter step—still driven by consistency arguments—one could impose

\begin{matrix} - {\bar{q}}^{n + 1} n + 2 μ ϵ ({\bar{v}}^{n + 1}) \cdot n + γ {\bar{v}}^{n + 1} = γ v^{n + 1}, x \in Γ_{R} . \end{matrix}

Remark 2.

While Equation (24) is required to prove that the solution of the EFR algorithm (19)–(28) converges to the solution of the NSEs (4) and (5) [6], it can be disregarded in the implementation of the algorithms. Releasing constraint (24) leads to substantial simplification and computational time savings since there is one less variable (i.e., the Lagrange multiplier to enforce the incompressibility constraint). In fact, one would solve a simplified filter problem: find filtered velocity

{\bar{v}}^{n + 1}

such that

\begin{matrix} - \nabla \cdot (δ \nabla {\bar{v}}^{n + 1}) + {\bar{v}}^{n + 1} & = v^{n + 1}, \end{matrix}

(30)

As noted in [70], the incompressibility is exactly preserved by the simplified differential filter (30) only for periodic conditions. In all the other cases, the end-of-step velocity

u^{n + 1}

does not strictly satisfy mass conservation. However, it can be shown numerically that at the discrete level, the mass conservation error is very low (see, e.g., [71]).

EFR algorithms have been well studied for the incompressible Navier–Stokes equations with a finite element [50,59,60,61,62,63,70], spectral element [52,53,54], or finite volume method [58]. They can be seen as splitting schemes for a nonlinear version of the Leray model. The original Leray model was introduced in 1934 as a theoretical tool to prove the existence of weak solutions of the NSEs [72]. As a computational tool, Leray regularization was first used in [73] as a stabilization strategy for under-resolved simulations of turbulent flows with classical numerical discretizations [6]. It was also shown [39,74] that when a differential filter is used, the Leray model is similar to the NS-

α

model of Foias, Holm, and Titi [75]. The EFR algorithm with nonlinear differential filters was shown to provide numerical results that are more precise in localizing where the eddy viscosity is needed and are overall more accurate than results obtained with Smagorinsky-type models while featuring comparable computational efficiency [60,61,62,63,70]. In [58,59], the results given by the EFR algorithm for certain filter choices were shown to be very accurate when compared to measurements for complex 3D incompressible flow problems.

2.1.2. The EFR Method for the Weakly Compressible Euler Equations

Despite the encouraging results for incompressible flows, EFR algorithms have been used to simulate compressible flows [51] only recently. While ideas similar to filter stabilization have been applied to the Euler equations [76,77,78,79], its full potential for compressible flows remains to be tapped into.

Below, we apply the EFR algorithm (11)–(13) to the weakly compressible Euler equations in the formulation most used in the literature of atmospheric simulations (see, e.g, [80,81,82] and references therein). Let us state such a formulation first.

Let

ρ

be the air density,

u

the wind velocity, p the atmospheric pressure, and

θ

the potential temperature. Temperature

θ

is a quantity of interest in meteorology and is defined as

θ = T / π

, where T is the absolute temperature and

π

is the Exner pressure. Since one is typically interested in fluctuations of

ρ

and p with respect to a constant (in time) and hydrostatically balanced background state

ρ_{0}

and

p_{0}

, we write

p = p_{0} + p^{'}

and

ρ = ρ_{0} + ρ^{'}

, with

\nabla p_{0} + ρ_{0} g \hat{k} = 0

, where g is the gravitational constant and

\hat{k}

is the unit vector aligned with the vertical axis z. See, e.g., [83] for more details on this. Then, the conservation of mass, momentum, and potential temperature in spatial domain

Ω

and over a time interval

(0, t_{f}]

can be written as

\begin{matrix} \frac{\partial ρ}{\partial t} + \nabla \cdot (ρ u) = 0 & in Ω \times (0, t_{f}], \end{matrix}

(31)

\begin{matrix} \frac{\partial (ρ u)}{\partial t} + \nabla \cdot (ρ u \otimes u) + \nabla p^{'} + ρ^{'} g \hat{k} = f_{u} & in Ω \times (0, t_{f}] . \end{matrix}

(32)

\begin{matrix} \frac{\partial (ρ θ)}{\partial t} + \nabla \cdot (ρ u θ) = f_{θ} & in Ω \times (0, t_{f}], \end{matrix}

(33)

where

f_{u}

and

f_{θ}

are possible forcing terms. We close system (31)–(33) with a thermodynamics equation of state for p:

\begin{matrix} p = p_{g} {(ρ R θ / p_{g})}^{c_{p} / c_{v}}, \end{matrix}

(34)

where

p_{g} = 10^{5}

Pa is a standard pressure, R is the specific gas constant of dry air, and

c_{p}

(respectively,

c_{v}

) is the specific heat capacity at constant pressure (respectively, volume). Obviously, system (31)–(34) have to be supplemented with proper initial and boundary conditions. However, we do not specify them, as they would be treated as explained in Section 2.1.1. Other perspectives on this are included in Section Open Problems.

EFR algorithm (11)–(13), with filtering as in (16), applied to problem (31)–(34) after time discretization by the backward Euler scheme reads as follows:

-: Evolve: find density $ρ^{n + 1}$ , density fluctuation $ρ^{', n + 1}$ , and intermediate variables $v^{n + 1}$ , $q^{n + 1}$ , $q^{', n + 1}$ , and $β^{n + 1}$ such that

$\begin{matrix} \frac{ρ^{n + 1} - ρ^{n}}{Δ t} + \nabla \cdot (ρ^{n + 1} v^{n + 1}) = 0, \end{matrix}$

(35)

$\begin{matrix} \frac{ρ^{n + 1} v^{n + 1} - ρ^{n} u^{n}}{Δ t} + \nabla \cdot (ρ^{n + 1} v^{n + 1} \otimes v^{n + 1}) + \nabla q^{', n + 1} + ρ^{', n + 1} g \nabla ρ^{n + 1} = f_{u}^{n + 1}, \end{matrix}$

(36)

$\begin{matrix} \frac{ρ^{n + 1} β^{n + 1} - ρ^{n} θ^{n}}{Δ t} + \nabla \cdot (ρ^{n + 1} v^{n + 1} β^{n + 1}) = f_{θ}^{n + 1}, \end{matrix}$

(37)

$\begin{matrix} q^{n + 1} = p_{g} {(ρ^{n + 1} R β^{n + 1} / p_{g})}^{c_{p} / c_{v}}, \end{matrix}$

(38)

$\begin{matrix} q^{n + 1} = p_{0} + q^{', n + 1}, \end{matrix}$

(39)

$\begin{matrix} ρ^{n + 1} = ρ_{0} + ρ^{', n + 1} . \end{matrix}$

(40)
-: Filter: find filtered variables ${\bar{v}}^{n + 1}$ and ${\bar{β}}^{n + 1}$ such that

$\begin{matrix} - \nabla \cdot (δ_{v} \nabla ({\bar{v}}^{n + 1})) + {\bar{v}}^{n + 1} = v^{n + 1}, δ_{v} = α^{2} a (v^{n + 1}), \end{matrix}$

(41)

$\begin{matrix} - \nabla \cdot (δ_{β} \nabla ({\bar{β}}^{n + 1})) + {\bar{β}}^{n + 1} = β^{n + 1}, δ_{β} = α^{2} a (β^{n + 1}) . \end{matrix}$

(42)
-: Relax: find end-of-step $u^{n + 1}$ , $p^{n + 1}$ , $p^{', n + 1}$ , and $θ^{n + 1}$ such that

$\begin{matrix} u^{n + 1} = (1 - χ) v^{n + 1} + χ {\bar{v}}^{n + 1}, \end{matrix}$

(43)

$\begin{matrix} θ^{n + 1} = (1 - ξ) β^{n + 1} + ξ {\bar{β}}^{n + 1}, \end{matrix}$

(44)

$\begin{matrix} p^{n + 1} = p_{g} {(ρ^{n + 1} R θ^{n + 1} / p_{g})}^{c_{p} / c_{v}}, \end{matrix}$

(45)

$\begin{matrix} p^{n + 1} = p_{0} + p^{', n + 1}, \end{matrix}$

(46)

with relaxation parameters $χ, ξ \in [0, 1]$ .

This particular version of the EFR method for compressible flows was introduced in [51]. Therein, the same filter was applied to the intermediate velocity and intermediate potential temperature. However, different variables have different physical properties and, thus, different types of filtering may be a better choice.

Similar to what was performed for (23), one can multiply (41) by

ρ^{n + 1} / Δ t

to obtain a Stokes-like problem:

\frac{ρ^{n + 1}}{Δ t} ({\bar{v}}^{n + 1} - v^{n + 1}) - \nabla \cdot (μ_{a} \nabla {\bar{v}}^{n + 1}) = 0, μ_{a} = ρ^{n + 1} \frac{α^{2}}{Δ t} a (v^{n + 1}),

(47)

where

μ_{a}

is dimensionally a dynamic viscosity, i.e., a proper eddy viscosity.

2.1.3. Indicator Function

Once the filter is selected, the success of the EFR algorithm depends on the reliability of the indicator function. This function should take values close to zero where its argument (i.e., velocity or potential temperature) does not need regularization and close to one where the argument needs to be regularized. Linear filters, which are equivalent to setting

a (v) = 1,

(48)

have limited efficacy, since they introduce the same amount of regularization everywhere in the domain and likely lead to overdiffusion. Much attention has been devoted to devising accurate indicator functions.

Before presenting a brief survey of the indicator functions available in the literature, we remark that one can recover a Smagorinsky-like model from the EFR algorithm by taking the following indicator function [84]:

a_{S} (v) = \frac{| \nabla v |}{{∥ \nabla v ∥}_{\infty}} .

(49)

This indicator function is convenient because of its strong monotonicity properties. However, it is known to introduce excessive diffusion, just like (48). In fact, notice that

a_{S}

takes large values in the regions of the domain with shear flow, which is not turbulent.

Below, we group indicator functions into two categories (physical phenomenology-based and mathematics-based) and present their strengths and limitations.

Physical Phenomenology-Based Indicator Functions

In [62], the authors propose indicator functions based on physical quantities that are known to vanish for coherent flow structures, i.e., structures that persist in their form for “long” periods of time. One of the most popular methods to identify coherent vortices is the Q criterion [85], which classifies as persistent and coherent vortices those regions where

Q (u, u) = \frac{1}{2} (\nabla^{s} u : \nabla^{s} u - ϵ (u) : ϵ (u)) > 0,

where

\nabla^{s} u

is the spin tensor (skew-symmetric part of

\nabla u

\nabla^{s} u = \nabla u - ϵ (u) = \frac{1}{2} (\nabla u - \nabla u^{T}) .

(50)

It is easy to see that

Q > 0

in the regions of the domain where spin dominates deformation. An indicator is obtained by rescaling

Q (u, u)

so that the condition

Q (u, u) > 0

implies

a (u) ≃ 0

, which means that regularization is not needed. A possibility is the following:

\begin{matrix} a_{Q} (u) = \frac{1}{2} - \frac{1}{π} arctan (ϵ^{- 1} \frac{Q (u, u)}{| Q (u, u) | + ϵ^{2}}), \end{matrix}

where

ϵ

is a parameter introduced to avoid the denominator vanishing.

A second indicator uses an eddy viscosity coefficient formula proposed by Vreman [86] that vanishes for 320 types of flow structures known to be coherent. The Vreman-based indicator function reads

\begin{matrix} a_{V} (u) = \sqrt{\frac{B (u)}{{| \nabla u |}_{F}^{4}}}, \end{matrix}

where the subindex F refers to the Frobenius norm and

B (u)

is defined as

\begin{matrix} B = β_{11} β_{22} - β_{12}^{2} + β_{11} β_{33} - β_{13}^{2} + β_{22} β_{33} - β_{23}^{2}, β_{i j} (u) = \sum_{m = 1, 2, 3} \frac{\partial u_{i}}{\partial x_{m}} \frac{\partial u_{j}}{\partial x_{m}} . \end{matrix}

Since

0 \leq B (u) / {| \nabla u |}_{F}^{4} \leq 1

a_{V} (u) \in [0, 1]

. The Vreman-based indicator function was shown to be successful in [60].

Another physics-based indicator function uses the relative helicity density

R H

, which is a local quantity related to the macroscopic helicity H:

\begin{matrix} R H = \frac{u \cdot w}{| u | | w |}, H = \frac{1}{| Ω |} \int_{Ω} u \cdot w d Ω, \end{matrix}

where

w = \nabla \times u

denotes vorticity. From the NSE in rotational form, it is possible to see that local high helicity suppresses local turbulent dissipation caused by the breakdown of eddies into smaller ones. A helicity-based indicator can be created by ensuring that values of

R H

near one imply

a (u) ≃ 0

, i.e.,

\begin{matrix} a_{H} (u) = 1 - |\frac{u \cdot w}{| u | | w | + ϵ^{2}}|, \end{matrix}

where again

ϵ

is a parameter that prevents the denominator from vanishing.

Notice that other (more selective) indicator functions can also be obtained by taking the geometric average of two (or more) indicator functions.

All these physical phenomenology-based indicator functions have the advantage of requiring only algebraic operations on

u

and its derivatives. Hence, their implementation may be quite straightforward. However, the major drawback is that they do not allow for a rigorous convergence theory to verify the robustness of the associated filtering method. This limitation is overcome by the mathematics-based (rather than physics-based) indicators discussed next.

Mathematics-Based Indicator Functions

In order to define a class of mathematics-based indicator functions, let us recall a basic result about linear filters

F

, which are invertible, self-adjoint, compact operators from a Hilbert space V to itself. From the spectral theorem (see, e.g., [65]), we have

\begin{matrix} F x = \sum_{i = 0}^{\infty} λ_{i} 〈 x, e_{i} 〉 e_{i}, F^{- 1} y = \sum_{i = 0}^{\infty} \frac{1}{λ_{i}} 〈 y, e_{i} 〉 e_{i}, \end{matrix}

where

λ_{i}

are the eigenvalues of

F

and

e_{i}

are the corresponding eigenfunctions, which form an orthonormal basis for V. Since

F

is compact, its eigenvalues accumulate to 0 and the inverse operator

F^{- 1}

is unbounded.

Let

D

be a bounded, regularized approximation of

F^{- 1}

, whose action on

y

is given by

D y = \sum_{i = 0}^{\infty} ϕ (\frac{1}{λ_{i}}) 〈 y, e_{i} 〉 e_{i}

with

ϕ (\frac{1}{λ_{i}}) ≃ \{\begin{matrix} \frac{1}{λ_{i}} & if i is “ small ”, \\ 0 & if i is “ large ” . \end{matrix}

Then, by composing

F

and

D

, we can obtain a low-pass filter, so that

∥ x - D F x ∥ is \{\begin{matrix} small & if x is “ smooth ”, \\ large & if x is not “ smooth ”, \end{matrix}

where “smooth” is intended with respect to the eigenfunction of the operator

F

. This means that

x

is smooth if

〈 x, e_{i} 〉

is significantly different from zero only for small values of i. The composition of the two operators

F

and

D

can be interpreted as a low-pass filter. This motivates the definition of indicator function

\begin{matrix} a_{D} (u) = \frac{|u - D (F (u))|}{{||u - D (F (u))||}_{\infty}} . \end{matrix}

(51)

A well-studied choice for the linear filter

F

in (51) is the linear Helmholtz filter, i.e., (16) with

δ

constant in space and time. Let us denote it with

F_{H}

. A popular choice for

D

is the Van Cittert deconvolution operator

D_{N}

, defined as

D_{N} = \sum_{n = 0}^{N} {(I - F)}^{n} .

Note that the evaluation of

a_{D}

with

D = D_{N}

requires the application of the linear filter

F

a total of

N + 1

times. This is an argument in favor of limiting the value of N, since large N would increase the computational cost. Another argument is provided by interpreting

D_{N}

as the N-th iteration of a Richardson scheme to solve the filter problem. See [59] for more details. In practice, one takes

N = 0, 1

, corresponding to

D_{0} = I

and

D_{1} = 2 I - F

, respectively. For these choices of N, the indicator function (51) becomes

\begin{matrix} a_{D_{0}} (u) = \frac{|u - F (u)|}{max (1, {||u - F (u)||}_{\infty})}, a_{D_{1}} (u) = \frac{|u - 2 F (u) + F (F (u))|}{max (1, {||u - 2 F (u) + F (F (u))||}_{\infty})} . \end{matrix}

(52)

Indicator function (51) with

D = D_{N}

and

F = F_{H}

, which we will refer to as the approximate deconvolution (AD) indicator function, was proposed in [61]. However, the idea of using van Cittert approximate deconvolution in fluid models to increase accuracy is well established and has strong mathematical foundations [87,88,89]. It is possible to prove [89] that

\begin{matrix} ϕ - D_{N} (F_{H} (ϕ)) = {(- 1)}^{N + 1} δ^{2 N + 2} Δ^{N + 1} F_{H}^{N + 1} ϕ . \end{matrix}

(53)

From this, it follows that, indeed,

a_{D_{N}} (u)

takes values close to zero in the regions of the domain where

u

is smooth.

In addition to allowing for rigorous mathematical analysis, the AD indicator function has been shown to be more selective than indicator

a_{S} (\cdot)

(49) in identifying the regions of the domain where diffusion is needed for both incompressible [61] and compressible flows [51]. Hence, the EFR algorithm with indicator function

a_{D} (\cdot)

corrects the overdiffusivity of

a_{S} (\cdot)

. For incompressible flows,

a_{D} (\cdot)

has also been shown to be particularly accurate for realistic 3D problems [8,58,59,64].

2.2. Machine Learning for LES

In this section, we briefly outline several machine learning strategies that have been recently used in the construction of LES models. While not as well developed as classical LES models, the machine learning strategies are a dynamic research area, with a lot of potential for LES, reduced-order modeling, and CFD in general. Thus, in this section, we outline several machine learning approaches used in LES and reduced-order modeling.

Recent years have seen an escalating interest in leveraging machine learning techniques for CFD, as elucidated in [3,90]. This subsection is devoted to highlighting some opportunities offered by machine learning for LES modeling, starting with a literature review. Noteworthy endeavors within various research groups have demonstrated the feasibility of training deep learning-based closure models [23,91,92,93], while others have focused on the acquisition of closure models through symbolic regression techniques [94,95,96]. It is imperative to underscore that a relatively underexplored domain is the integration of blending and classification methodologies, particularly utilizing indicator functions as presented in the previous subsection.

There is an existing body of knowledge detailing the efficacy of various closure models for specific types of flows. Furthermore, current efforts are marked by a growing emphasis on generating extensive turbulent flow datasets [97]. However, a critical gap remains in establishing comprehensive maps between underlying patterns and flow structures, along with the identification of the best-suited, verified, tested, and mathematically analyzed closure models. On a more granular level, once trained, the same maps can also be used to find the most suitable filters or filtering parameters. This approach allows the utilization of indicator functions within certain latent spaces to accurately identify a low-pass filter or its associated shape and parameters, which can then be applied to mainstream closure approaches. In essence, machine learning tools can be leveraged not only to identify a filter but also to build a hybrid closure model atop it. Furthermore, it is important to note that references [98,99] demonstrate how a single hypothesis or map can be used to create a simple selection mechanism. This mechanism proves valuable for choosing between different structural and functional closures or even various numerical methods.

A hypothesis segregation concept proposed in [98] uses functional and structural closure models and in [99] tests upwind and central numerical methods. The key idea is to train a modular neural network for structure identification offline and use that model to define an expert system to adjust underlying closure models and numerical approaches in a weighted manner temporally or spatially online as needed. A key aspect is the modular nature of offline–online split approaches, where the same trained model defines filter strengths, radii, and shapes for use in a structural model. See Figure 3 for a schematic representation.

We also advocate for the judicious use of machine learning tools to facilitate the generation of diverse maps within the latent space [100,101,102]. Such learned maps serve as foundational elements for constructing an expert system capable of discerning flow regimes and their associated closure models from an array of existing models. This approach offers a departure from the conventional approach of relying purely on black-box closure models or resource-intensive symbolic regression practices.

By harnessing the pattern and structure recognition capabilities inherent in machine learning, one can effectively employ existing structural or functional closure models in a dynamic and adaptive manner. For example, an approach with two types of artificial neural networks for ML has been proposed [103]: a classifier to identify the contribution of each building block in the flow and a predictor to estimate wall shear stress by combining these building-block flows. Such approaches with structure recognition capabilities are poised to enhance the efficiency and accuracy of flow regime identification, thereby presenting a compelling alternative to traditional methodologies. The use of flow physics-based indicator functions brings the additional advantage of exploiting existing structural and functional knowledge in a dynamic and adaptive framework. We conclude this section by referring to recent survey articles that provide further insights into machine learning for turbulence modeling [26,104].

3. Applications of LES for FOM

In this section, we present some of the achievements of the LES concepts and models presented in Section 2, specifically spatial filtering, approximate deconvolution, and the evolve–filter–relax strategy. We present these achievements for both incompressible flows (Section 3.1) and compressible flows (Section 3.2). We believe that presenting applications for both incompressible and compressible flows is important in order to demonstrate that these LES concepts and models are truly universal and are not restricted to particular classes of flows or computational settings. Furthermore, by showcasing the achievements of spatial filtering, approximate deconvolution, and the evolve–filter–relax strategy in challenging applications in biomedicine and weather modeling, rather than in academic benchmarks, we also illustrate the potential that these concepts and models have in reduced-order modeling.

3.1. Incompressible Flows

The use of LES models for incompressible flows covers a wide range of applications. Among these, we want to mention here biomedical applications, with particular reference to computational hemodynamics. This is a field of CFD that has progressively received more attention with the availability of patient-specific data and reliable methods of image processing-enabled numerical methods from the benchmark to the bedside/operating room stage [105]. Using CFD in the clinical routine for the assessment of cardiovascular diseases—and sometimes also for the planning of therapeutic options—is nowadays a consolidated practice, with some outstanding examples involving entrepreneurial activities [106,107]. This circumstance clearly raises some specific challenges that go beyond the mathematical modeling stage, involving the efficiency required by clinical timelines and quality certification/uncertainty quantification.

In general, human blood flow does not feature a turbulent regime in most of the vascular districts of interest. However, in some cases, highly disturbed flow is possible either for physiological reasons (like in the ascending aorta during the systolic phase, when the aortic valve is open and the heart pumps blood into the system) or for pathological conditions, like in the presence of significant stenoses. In some cases, turbulence may be induced by devices required by therapeutical options. In this section, we focus on an application involving a pathological condition called type B aortic dissection (AoD) because it highlights the practical usefulness of the EFR method in a complex real-life setting and the associated challenges.

3.1.1. Computational Hemodynamics in Type B Aortic Dissections

AoD is a pathology consisting of a splitting of the natural lumen of the descending aorta into two separate channels: true lumen and false lumen—see Figure 4. This splitting may have different reasons, including traumatic ones, and represents a danger for the patient: the tissue of the two lumens is weaker than the physiological one, so it may break, causing life-threatening hemorrhages. Clinically, there is a strong interest in predicting the evolution of the false lumen in time, in particular whether it remains stable or grows, the latter condition being the most dangerous for the patient. While drugs are an option for a stable false lumen, surgery is the usual choice for cases in which growth is expected. Since surgery is invasive, it should be considered only when necessary.

Currently, decisions are made only on the geometry of the dissection. The addition of hemodynamics information, in particular of quantities that are difficult or impossible to measure in vivo, is anticipated to improve significantly the reliability of the predictions [8]. In particular, the wall shear stress (WSS), i.e., the tangential component of the normal stress at the wall, may play a significant role for prognostic purposes. It is defined as

τ = 2 μ ϵ (u) \cdot n - 2 μ (n \cdot ϵ (u) \cdot n) n on Γ,

(54)

where

Γ

is a portion of

\partial Ω

of interest and

n

is its outward normal unit vector. Hence, CFD of aortic dissections may become a critical component in the clinical routine. However, the complexity of patient-specific geometries and the high convective fields induced by the small entry tears between the true and the false lumen (see Figure 5) pose significant challenges. It is important to stress that in clinical applications, there is the need to run many patient-specific cases, as required by the concept of computer-aided clinical trials (CACTs), to obtain statistically significant data. This means that beyond accuracy, computational efficiency is a critical aspect. In the case of AoD, a one-shot simulation might be possible with DNS, but the computational cost is not compatible with the large number of patients required by CACTs. The trade-off between efficiency and accuracy—where accuracy has to be intended in the clinical sense, i.e., the quantification of indexes to enable correct clinical decisions—is of primary importance, and the EFR scheme, in this context, seems to be the ideal approach.

Geometry

The geometries of AoD may be significantly diverse—see Figure 6. The location of the first tear connecting the true and the false lumen (also known as primal entry tear—PET) has a critical importance. In fact, clinical guidelines suggest that the evolution of the dissection depends on the distance between the PET and the left subclavian artery. A short distance generally indicates the evolution into complicated clinical cases requiring surgery [108]. However, there are cases with a distal PET evolving in cases to be surgically operated [8]. For certain sites of interest in the cardiovascular system (e.g., the coronary arteries), the geometry can be generally templatized, i.e., the patient-specific anatomy can be regarded as a map or a minor modification of a template geometry. In AoDs, it is very hard to find a realistic template. While AI-based methodologies can be used for the segmentation and reconstruction of coronaries, at this time, the classical level-set method is the method of choice for complex geometries like AoD. In fact, the lack of sufficiently large datasets prevents reliable training of data-based methodologies, like neural networks. The geometries presented in this paper were reconstructed with the Vascular Modeling ToolKit [109], a free, open-source tool based on the level-set method, using images from computer tomography. The segmentation was carried out by an expert, trained for the identification of the tears. The mesh was successively obtained by the NetGen mesh generator, with an in-home Python script to guide the refinement of the mesh close to the boundaries (the so-called mesh boundary layers). The introduction of a mesh more refined close to the boundary is a common practice in biomedical applications whenever one of the main quantities of interest is computed at the boundary, like the WSS. The meshes we deemed to be appropriate for simulations of AoD with EFR turbulence modeling have a number of elements in the range

6 \times 10^{5} - 6 \times 10^{6}

Boundary Conditions and Backflows

The patient-specific domain

Ω

generally features physical boundaries

Γ_{w}

represented by the arterial walls and artificial boundaries represented by the regions where the blood enters (

Γ_{i n}

) and leaves (

Γ_{o u t}

) the region of interest

Ω

. The arterial walls include also the flow divider separating the true and the false lumens. On the wall, we postulate the nonslip condition prescribing the continuity of the velocity field of the wall and the blood. Since we assume a rigid structure, this leads to the homogeneous Dirichlet condition

u = 0, on Γ_{w} .

The role of fluid–structure interaction, in particular at the flap separating the two lumens, has recently attracted the attention of the biomedical engineering community [111,112]. We considered rigid structures for simplicity. However, nothing prevents the extension of the framework presented here to the case of compliant structures (see, e.g., [113]).

When considering the artificial inflow/outflow portions of the boundaries, computational hemodynamics applied to clinical problems generally suffers from a long-term condition, sometimes called the “defective boundary conditions problem” [114] or, in the specific case of medical applications, the “image-legacy problem”. Medical doctors used to keep track of images but not of possible hemodynamic conditions to be used as boundary data. This means that we generally have patient-specific geometrical data but with a substantial lack, if not an absence, of patient-specific boundary conditions. To leverage the abundance of geometrical data, many different approaches have been considered in the literature to fill this gap. Standard procedures consider (a) literature-based data for flow rates at the inflow, with an empirical selection of velocity profiles fitting these data; and (b) special boundary conditions at the outflow, based on a lumped parameter modeling of the peripheral circulation. We briefly detail these conditions hereafter.

Inflow conditions. The available data may generally refer to the flow rate or to the velocity (as a function of time) at one point of the inflow. While this is clearly not enough as a boundary condition, a popular and practical approach consists of assuming a velocity profile constructed around the available data. For instance, one can assume a parabolic profile (corresponding to the well-known Poiseuille solution) or a Womersley profile (corresponding to the time-dependent Womersley solution) [115] that fit the available flow rate or pointwise velocity. This clearly introduces a bias in the solution since the choice of the profile—even if educated—is arbitrary. To mitigate this aspect, an artificial extension of the inflow tract is added to the computational domain, called flow extension. An analysis of the error introduced by this approach can be found in [116]. More sophisticated and mathematically sound approaches were introduced with the idea of identifying the best profile through an optimization approach (see [116,117,118,119,120]), yet they require an additional computational cost (and nonstandard solvers)). In the simulation considered in our work on AoD [8,121], we assumed flow rate data provided by the literature, with flow extensions at the inlet (the ascending aorta) and the selection of a constant velocity profile. An analysis of the impact of idealized velocity profiles in AoD simulations can be found in [122].
Outflow Conditions. For outflow conditions on the portion of the domain denoted by $Γ_{o u t}$ , a common approach is to combine the 3D model based on the incompressible Navier–Stokes Equations (4) and (5) with surrogate—dimensionally reduced—models representing the downstream circulations, in what has been called the geometrically multiscale approach [123,124,125,126]. It is important to stress that in this particular case of AoD, we may have some outflow sections referring to collateral vessels like the renal or the mesenteric arteries where not only do we lack patient-specific data but also it is hard to find data in the literature. Yet, the inclusion of the renal flow is critical for a reliable assessment of the hemodynamics in an AoD. A popular approach, in this case, is to resort to the introduction of a lumped parameter model called three-element Windkessel (Windkessel was a device used by firefighters to pump water from a reservoir, converting a periodic action into a quasi-steady flow, which is exactly what happens in the peripheral circulation, where the pulsatile aortic flow is eventually converted into a steady flow in the capillaries [127]), representing the downstream circulation at each outflow boundary (see, e.g., [127]). This approach leads to the prescription of a traction condition of the form

$p n - μ (\nabla u + \nabla^{T} u) \cdot n = P_{w k} n, on Γ_{o u t},$

(55)

where $P_{w k}$ is a function of time defined as

$P_{w k} (t) = R_{1} Q (t) + e^{- \frac{t}{R_{2} C}} (P_{w k} (0) - R_{1} Q (0)) + \int_{0}^{t} \frac{e^{\frac{τ - t}{R_{2} C}}}{C} Q (τ) d τ .$

Here, $Q (t) = \int_{Γ_{o u t}} u \cdot n d γ (t)$ , and $R_{1}$ , $R_{2}$ , and C are parameters representing the proximal and distal flow resistances and the compliance, respectively, of the downstream vascular district. The condition (55) turns out to be a sort of “Robin” condition (see Remark 1) for the NSEs, mixing values of velocity and pressure, as the result of the multiscale coupling between the NSEs and the Windkessel model. Each outflow section features its own set of parameters $R_{1}$ , $R_{2}$ , and C. The calibration of these parameters is usually performed with a combination of patient-specific/literature data with some empirical criteria (see, e.g., [128,129]). Although critical, this step is not central to the topics of this paper and was mentioned only for completeness.

It is worth pointing out that in aortic simulations, a numerical instability may occur on the boundary

Γ_{o u t}

where traction conditions, like (55), are prescribed. This circumstance goes under the name of “backflow instability,” and it is related to a lack of energy control in the solution when a significant amount of fluid enters through

Γ_{o u t}

. In detail, when computing the energy of the fluid, the boundedness of the Navier–Stokes solution can be proved under the condition that

\int_{Γ_{o u t}} ρ {| u |}^{2} u \cdot n d γ \geq 0 .

(56)

Unfortunately, the occurrence of incoming flow (that at outflow sections we can call backflows) such that

u \cdot n < 0

may impair the control on the boundedness of the data, eventually leading to the so-called backflow instabilities [130]. On the other hand, backflows are physiological in the aortic circulation. The numerical treatment of this problem has been discussed in several contributions, see, e.g., [130,131,132,133,134,135]. A natural, simple approach in the specific case of Windkessel conditions (55) consists of artificially increasing the proximal resistance

R_{1}

, as this results in damping the impact of the backflow on the energy balance. As we will see below, when using LES-EFR models, this turned out to be unnecessary.

EFR in Action

The occurrence of highly disturbed flow in the human ascending aorta is the consequence of the high inlet velocity at the aortic valve [136]. However, the alternation of fast and slow transients associated with the different phases of the heartbeat (“systole” when the aortic valve is open, “diastole” when it is closed) generally prevents the transition to turbulence in physiological cases. Nonetheless, in the presence of small tears where the blood is strongly accelerated, we may have the transition to turbulence. Some studies on AoD report peak Reynolds numbers exceeding 10,000 [137] and consider

k - ω

modeling for the turbulence. Other studies [138] find peak Reynolds numbers below 5000 and consider laminar flow appropriate. This is a consequence of the fact that in computational hemodynamics, the patient-specific diversity and the high level of uncertainty may lead to different findings. As we mentioned earlier, CACTs require simulations for a reasonably large number of patients. Hence, computational efficiency is a criterion of paramount relevance, as is the capability to capture transitional and turbulent regimes. The use of LES modeling, and in particular EFR, not only is appropriate for the expected Reynolds numbers but also allows for a significant reduction of the computational costs because it enables the use of coarse meshes that do not necessarily solve all the significant scales.

In [8], we compared a DNS for a patient-specific case with a mesh with 1,063,000 elements vs. an EFR simulation with a mesh with 654,000 elements. Space discretization was performed with classical Taylor–Hood inf–sup compatible finite elements within the C++ library LifeV [139]. The results obtained were similar and equivalent from the clinical point of view. The reduction of the computational time, however, was significant (around

80 %

For the Dirichlet conditions on the wall and at the inflow, we used the conditions for

v

and

\bar{v}

, as specified in Section 2.1.1. For the Windkessel conditions (55), after the time discretization, we have two options:

Implicit, i.e., we actually evaluate $Q^{n + 1} = \int_{Γ_{o u t}} u^{n + 1} \cdot n d γ$ so that we resort to Robin boundary conditions similar to the ones discussed in Remark 1;
Explicit, i.e., we use a time extrapolation of the velocity consistent with the time-discretization accuracy. For a first-order time advancing, for instance, one has $Q^{n + 1} \approx \int_{Γ_{o u t}} u^{n} \cdot n d γ$ , leading to classical traction (Neumann) conditions.

In the former case, for the filtered velocity, we follow the strategy indicated in Remark 1, and in the latter one, we adopt the homogeneous Neumann conditions specified in (26).

Filter Radius and Relaxation Parameter

The selection of the parameters

α

(radius) and

χ

(relaxation) is obviously critical and still an open problem. Ideally, one would like to have an automatic selection of these parameters. This calls for the identification of general criteria to be applied to the specific case of interest. In the case of AoD, we identified many different geometrical and fluid dynamic conditions that make the application of general criteria quite problematic. With extensive experience in different cases, we identified some basic principles and heuristic criteria (see [59,110,121]).

The radius

α

is calibrated to be related to the mesh size h. However, we found that setting

α = h_{m a x}

, i.e., the maximum size of the elements of the finite element mesh, may generate unnecessary dissipation. We decided to set

α = h_{m i n}

(smallest edges of elements in the mesh) to confine the effects of the diffusion introduced by the EFR method.

When selecting

α

proportional to the mesh size, a consistency argument for the numerical EFR solution with the original equations suggests we take

χ \propto Δ t

[59]. In particular, it has been observed that if we set

χ = \frac{μ Δ t}{c_{0} ρ η α^{2}} max (h - η, 0),

(57)

we obtain that the viscous term of the EFR model is comparable to the viscous term in a standard NS solver when the mesh size is the Kolmogorov scale

η

While this approach to set

χ

is theoretically intriguing, in the case of AoD we need to consider the additional problem of backflow instabilities. The following result was proven [121]:

Theorem 1.

If the data of the problem are regular enough, assuming homogeneous boundary data, for

χ \to 1^{-}

and a suitable choice of the radius α, then

∥ u^{n} ∥_{L^{2}} \leq {∥ u_{0} ∥}_{L^{2}} .

This theorem has an important consequence: the energy of the system, even in the presence of backflows, decreases since the diffusive term introduced by the EFR scheme controls the term (56) responsible for the instability. In order to prevent backflow instabilities, the calibration of

χ

based on (57) needs to be adjusted. This is illustrated in Table 1, referring to a geometry designed to trigger backflows (a domain represented by a portion of a torus with a decreasing radius).

In general, a larger

α

allows for a smaller

χ

and vice versa. However, the stabilization of backflows is obtained for values of the relaxation parameter independent of

α

and with a sublinear dependence on the time step.

A specific sensitivity analysis to elucidate the impact of the value of

α

on the simulation is reported below.

Results

With the method described above, we simulated a number of patients in search of possible hemodynamics predictors that could improve the clinical guidelines to diagnose the evolution of the dissection. This work, still in progress [140], is in collaboration with the Department of Surgery at Emory University (Dr. B. Leshnower). While the biomedical implications of our results are not the focus of this paper, we summarize here some major findings as evidence of the practical success of the EFR methodology. In [8], we considered a longitudinal study of a patient presented at Emory University in 2006 and followed for several years. Imaging in 2010 evidenced a significant growth of the false lumen, even though the PET was distal. We co-registered the images in 2006 and 2010 to quantify the false lumen growth and correlate this growth with the hemodynamics computed with our method in the 2006 geometry. The results pointed out a potential role of the time-averaged WSS (TAWSS) defined as

TAWSS \equiv \frac{1}{T} ∥\int_{0}^{T} τ (t) d t∥,

where T is the heartbeat duration. The growth is proportional to the TAWSS, but only below a certain threshold. Other hemodynamic indexes were found to be correlated with the possible thrombotic activity leading to a local shrinking of the false lumen (see [110] for more details). After analyzing several more patients with both growth and no growth, we are now focusing on the difference between the TAWSS in the false and true lumens as a potential predictor of pathology evolution. See Figure 7.

Sensitivity Analysis for the Filtering Radius

Sensitivity analysis can be carried out in different ways. A possible approach is an experimental, yet educated, sampling of the solution for different values of the parameters. In general, this requires several samplings, which, in our case, means running the simulation for several values of

α

Local sensitivity analysis consists of quantifying the dependence of the solution on a parameter by computing the corresponding (Gateaux) derivatives. These are the so-called sensitivities. The numerical solution of these equations for a particular value of the parameter allows us to determine in what parts of the computational domain the solution is more sensitive to such a parameter. This is a local analysis since it applies to a particular value of the parameter. A complete derivation of the sensitivity equations to the radius

α

for Leray linear and nonlinear spatial filtering problems with AD indicator functions (48) and

a_{D_{0}}

can be found in [141].

The results in [141] show that the velocity sensitivity can be a reliable predictor of the regions in the domain with high sensitivity to a specific parameter value. In simulating the flow past an obstacle, for instance, the overdiffusive effects immediately past the obstacle due to an excessively large value of

α

are correctly identified by the velocity sensitivity. This leads to a possible “adaptive” algorithm for the calibration of the radius. Denoting by

s_{u} \equiv \frac{\partial_{u}}{\partial_{α}}

and by

s_{p} \equiv \frac{\partial_{p}}{\partial_{α}}

the sensitivities of velocity and pressure, one could solve the LES model and the sensitivity equations with a “large” value

α_{0}

of the radius to attain stability and then refine the solution with a Taylor expansion

u (α) \approx u (α_{0}) + s_{u} (α_{0}) (α - α_{0}), p (α) \approx p (α_{0}) + s_{p} (α_{0}) (α - α_{0}),

(58)

with

α < α_{0}

. In the same paper, the sensitivity of the solution to a scheme with (48) is demonstrated to be significantly higher than with the nonlinear AD filter, pointing out the superiority of the latter.

For specific cardiovascular applications, a global sensitivity analysis was performed too [142]. The purpose of this analysis is to understand what inputs affect the solution the most. The inputs are regarded as stochastic processes. Quantities of interest obtained from the numerical solution of the hemodynamics through our EFR method are consequently regarded as stochastic processes too. The sensitivity to the different inputs is quantified in terms of the variance of the outputs as a function of the input variances. In [8], three inputs were considered: the geometry of the domain, the inflow rate, and the radius

α

. For each of these variables, a probabilistic distribution was postulated. For the radius, a Gaussian distribution with average

h_{m i n}

and variance

0.25 h_{m i n}^{2}

was deemed to be a correct description of the range of interest for

α

. The uncertainty in the geometry comes from the noise in the images and the consequent errors in the patient-specific reconstruction. Two cases were considered: (i) an idealized aortic arch, where the uncertainty relies on the geometrical parameters of the ideal geometry; (ii) a patient-specific case, where the geometry was reconstructed in a time range of 6 years. In case (ii), after an image registration procedure, the geometry was parametrized by one parameter specifying the morphing instant from the initial to the final geometry. For this parameter, we postulated a uniform probabilistic distribution. To quantify the uncertainty, we used polynomial chaos expansion (PCE). In PCE, a variable-of-interest function of random variables is regarded as a polynomial expansion on a basis incorporating the stochastic features of the inputs. This basis is orthogonal with respect to a suitable probability measure. The coefficients of this (truncated) expansion can be computed with a pseudo-spectral approach through appropriate quadrature rules. This allows us to compute the so-called Sobol’ indexes, which are variance-based indicators to characterize the dependence of the output variance on the stochastic inputs [143,144] and rank the influence of the inputs, especially for nonlinear models [145,146]. In short, the Sobol’ index for a specific input i on a quantity of interest f is defined as

S_{i}^{f} \equiv \frac{Var (f_{i}^{P C})}{Var (f^{P C})},

where

Var (f^{P C})

denotes the variance of the stochastic process f represented by the polynomial chaos expansion and due to all the stochastic inputs, while

Var (f_{i}^{P C})

refers only to the variance induced by the input i.

In [142], the outputs of interest are the turbulent kinetic energy (TKE), the TAWSS, and the oscillatory shear index (OSI). The OSI is a hemodynamic index ranging between 0 and 1/2, which quantifies how long during the heartbeat the WSS is opposite to the direction of the flow. It is defined as

OSI \equiv \frac{1}{2} (1 - \frac{∥ \int_{0}^{T} τ (t) d t ∥}{\int_{0}^{T} ∥ τ (t) ∥ d t}) .

Through the application of PCE, we learned a couple of important points:

The geometry is by far the most important factor in the results: an accurate geometrical reconstruction of the region of interest is critical for any biomedical analysis. See Figure 8, rightmost panels.
The impact of the radius on the TAWSS and the OSI is minimal. See Figure 8, leftmost and center panels. This means that the selection of the radius with the empirical rules used in our simulations is not expected to have a major impact on the clinical conclusions of the computational analysis.

This corroborates the confidence in our EFR modeling for AoD.

3.1.2. Open Problems

The EFR method gave excellent results in the simulation of AoD, providing the trade-off between the computational efficiency and the accuracy needed to process a significant number of patients. Yet, there are some aspects of the method that still need to be elucidated.

A first important point is the optimal selection of the parameters

α

and

χ

. As mentioned earlier, the selection is largely based on empirical and physical criteria. An automatic, somehow optimal, selection would be certainly desirable. As noted, a possible approach for the radius comes from the local sensitivity analysis and the approximate expansion (58). An initial conservative choice of the radius can be refined at each time step so long as some indicators do not suggest the presence of instabilities. This approach has the clear drawback of requiring the solution of the sensitivity equations, which may be quite expensive in some applications, like AoD hemodynamics. Yet, it can provide the backbone for the definition of an adaptive scheme once the sensitivity computations can be properly surrogated. Likewise, the initial setting of

χ

can be based on linear dependence on

Δ t

and physical arguments, but the occurrence of backflows requires successive refinement. All these aspects need to be further investigated.

Another important problem in hemodynamics applications is related to “data assimilation” procedures that may improve the overall reliability of numerical simulations. With more accurate measurement tools, like 4D MRI, and more storage space available, more hemodynamics data are progressively more available, in particular velocity data in some points of the region of interest. In general, these data are not on the boundary, so they do not provide boundary conditions, yet they give important information about the blood dynamics. However, these data can barely replace the wealth of information provided by a numerical simulation, particularly for hemodynamic indexes like the TAWSS or the OSI. In a cooperative view, where data and models are not alternative but complementary, we need to find a solution both fitting the data (possibly loosely in the presence of high noise) and fulfilling the equations. The problem has been well-known for years (see, e.g., [147,148,149]), yet a “standard” solution may require even higher—unfeasible—computational costs. In this respect, the combination of physically informed neural networks (PINNs) [150] and our EFR scheme may provide a potential breakthrough that will be investigated in the years to come.

3.2. Compressible Flows

As mentioned in Section 2.1.2, despite many encouraging results for incompressible flows, EFR algorithms have been applied to simulate compressible flows only very recently in [51]. Therein, algorithm (35)–(46) was introduced and tested on classical benchmarks for atmospheric flow. Here, we report the main results to support the following findings: (i) for a given mesh, the EFR algorithm with the AD indicator function is less dissipative than a Smagorinsky-like method (EFR with indicator function

a_{S}

in (49)), (ii) the reason for point (i) is the fact that

a_{D}

is more selective in identifying the regions of the domain where artificial diffusion is needed, and (iii) the additional cost for using

a_{D}

instead of

a_{S}

is marginal. For all the results reported here, we set

χ = ξ = 1

so that the artificial diffusion introduced by the EFR algorithm can easily be calculated (see [51]). See Section Open Problems below for more on this.

To illustrate the above point (i), we first use results obtained for the well-known rising thermal bubble benchmark. This is a 2D benchmark that entails perturbing with a bubble of warm air a neutrally stratified atmosphere with uniform background potential temperature. In the time interval

(0, 1020]

s, the air warmer than the ambient rises due to buoyancy and deforms due to shearing motion into a mushroom shape. For the set-up, we refer to, e.g., [51,151,152]. The computational domain, which is

5 \times 10

Km², is discretized with two meshes with uniform resolution: a coarser one with

h = 62.5

m and a finer one with

h = 31.25

m. The time step is set to

Δ t = 0.1

s for all the simulations. Figure 9 shows the spatial distribution of

θ^{'}

t = 1020

s computed with the EFR algorithm, indicator functions

a_{S}

a_{D}

, and the two meshes. On a given mesh, the Rayleigh–Taylor instability at the edge of the bubble is more developed when using

a_{D}

instead of

a_{S}

, which indicates that

a_{D}

introduces less artificial viscosity than

a_{S}

. This is corroborated by the results for a more complex benchmark called density current [153,154] (see also [155,156,157] for LES of density currents in incompressible flows), in which a neutrally stratified atmosphere with uniform background

θ

is perturbed with a bubble of cold air. In the time interval

(0, 900]

s, the cold air descends due to negative buoyancy and when the cold air reaches the ground, it rolls up and forms a front. By time

t = 900

s, a three-rotor structure emerges as a result of the shear generated at the top boundary of the cold front. The results obtained with the EFR method and indicator functions

a_{S}

and

a_{D}

are shown in Figure 10. Note that

θ^{'}

computed with

a_{D}

is sharper and less smeared than

θ^{'}

computed with

a_{S}

. Additionally, the plots of the indicator functions show that

a_{S}

t = 900

s takes larger values over wider regions than

a_{D}

, illustrating the above point (ii).

The computational time to run the simulations in Figure 10 on a common laptop is reported in Table 2. Obviously, the total computational cost increases when switching from

a_{S}

a_{D}

because while

a_{S}

requires a simple post-processing of the velocity field,

a_{D_{0}}

in (52) requires one application of the linear Helmholtz filter. However, note that that increase amounts to only 0.006 s per time step and the filter step with

a_{D}

takes about half of the time needed for the evolve step. In total, the simulation takes only about 13 min.

We remark that a DNS for this benchmark would require over 5000 billion DOFs. The results presented here show that the EFR method with a deconvolution-based indicator function is much less dissipative than Smagorinsky-like LES methods while featuring comparable computational efficiency. These findings suggest that EFR is a promising approach for compressible flows, as the studies on incompressible flows suggested.

Open Problems

To obtain the results shown in Figure 9 and Figure 10, we set

χ = ξ = 1

(i.e., no relaxation),

N = 0

, and performed a costly tuning of

α

to obtain results that compare well with the literature [43,81,151,158,159]. A clear understanding of the complex interplay of these parameters and their effect on the solution accuracy is missing for compressible flows. In Section 3.1.1, we saw that even for incompressible flows, for which much more work has been undertaken, this understanding is still partial. A full understanding is obviously central to the success of the EFR algorithm with the AD indicator function and the definition of automatic tuning procedures. Another open question is the identification of the most physically accurate and computationally efficient filter

F

for each filtered variable. For this purpose, one could use various filters available in the literature (see, e.g, [5]).

Just like in the case of complex hemodyanmics simulations, care must be taken in setting the boundary conditions for compressible flow problems. Indeed, for compressible flows, the presence of artificial boundaries may trigger unphysical reflections that need to be suppressed by “nonreflecting” conditions. These conditions mimic the evolution of the solution outside the domain along the characteristics in order to suppress artificial inflow components. When nonreflecting boundary conditions are imposed, we need to devise conditions that minimize any artifact in the physical propagation. As the filter adds diffusivity, an initial idea could be to explore noninvasive conditions to minimize the impact of the filter. If this is not enough, one could explore different filter formulations.

Finally, we would like to stress the importance of data assimilation for weather/climate prediction too. In the atmospheric sciences, there is an abundance of data, which are collected either in situ or through remote sensing. Such data (or observations) are routinely incorporated in numerical weather prediction through data assimilation techniques, which seek to optimally combine the weather model with the observations. As mentioned in Section 3.1.2, it would be very interesting to study the computational efficiency and accuracy allowed by a combination of PINN and EFR.

4. LES for Reduced-Order Models

In this section, we build on the LES concepts, ideas, and models (i.e., spatial filtering, approximate deconvolution, and evolve–filter–relax strategy) introduced in Section 2 to construct LES-inspired ROMs, which we call LES-ROMs.

Over the past couple of decades, ROMs (see, e.g., [160,161,162,163,164,165] for reviews) have emerged as the methodology of choice to reduce the computational burden when traditional flow simulations (also called full-order models—FOMs) have to be carried out for several parameter values, as in the case of uncertainty quantification, optimal control, and inverse problems. ROMs replace the FOM with a lower-dimensional approximation that captures the essential behavior of the flow. The reduction in computational time is achieved through a two-step procedure. In the first step, called the offline phase, one constructs a database of several FOM solutions associated with given times and/or physical parameter values. The FOM database is used to generate a reduced basis, which is (hopefully much) smaller than the high-dimensional FOM basis but still preserves the essential features of the flow. One of the most popular reduced bases is proper orthogonal decomposition (POD), which extracts the dominant modes from a FOM database. We remark, however, that the strategies we discuss are not limited to POD and work equally well with other methods (e.g., the reduced basis method [163,166]). In the second step, called the online phase, one uses this reduced basis to quickly compute the solution for newly specified times and/or parameter values. Note that while the offline phase is performed only once, the online phase is performed as many times as needed.

ROMs are very efficient surrogate models when the number of reduced basis functions is small, as is typically the case for diffusion-dominated flows. Unfortunately, though, a large number of basis functions is needed to capture the essential features of convection-dominated flows. If one retains such a large number of modes in order for the ROM to be accurate, then the computational efficiency suffers. If the number of modes is otherwise kept low, a severe loss of information hinders the accurate reconstruction of the solution. In fact, Galerkin projection-based ROMs of turbulent flows are affected by energy stability problems related to the fact that, e.g., POD retains the modes biased toward large, high-energy scales, while the TKE is dissipated at the level of the small scales. A possible way to tackle this challenging problem is to introduce dissipation via a closure model [167,168,169] or stabilization [170,171,172,173,174].

We note, however, that the current ROM closures and stabilizations are typically heuristic and generally lack mathematical support. Traditionally, researchers propose ROM closures and stabilizations that make sense intuitively, test them in limited computational settings, and declare success once they yield more accurate solutions than the corresponding nonstabilized ROM. The main drawback of these ad hoc ROM closures and stabilizations is that while they work well for the specific flow configuration on which they were trained, they can fail in different flow configurations and need to be retrained (which is a costly, tedious process) in order to perform adequately.

Recently, there have been several new research developments aimed at constructing novel ROMs for realistic turbulent flows. The following are three significant research avenues in this direction: LES-ROMs (Section 4.1) aim at “bridging” two distinct research fields, LES and ROMs, to construct efficient and accurate ROMs for turbulent flows. FOM-ROM coupling (Section 4.2) aims at “bridging” the numerical methods used in the FOM and ROM in order to improve the stability and accuracy of the latter. Finally, the numerical analysis of the FOM-ROM framework (Section 4.3) leverages rigorous numerical analysis results to improve the ROM performance.

4.1. LES-ROMs

For the past several years, our groups and collaborators have proposed, analyzed, and investigated a new strategy for constructing ROMs for turbulent flows. This new strategy is based on spatial filtering, which is a mathematically and physically sound principle that is central in LES (see Section 2). This new strategy aims squarely at bridging two research fields, LES and ROMs, that have been treated separately until now. We also emphasize that this new strategy is fundamentally different from current ROM closures and stabilizations, which are often disconnected from the numerical methods used at the FOM level.

These new LES-inspired ROMs, which we call LES-ROMs, are capable of operating over a range of flow problems with minimal user intervention. Again, this is in stark contrast with current ROM closures and stabilizations, which need to be retrained for each new computational setting.

Just as in the FOM setting, the new LES-ROMs are constructed by leveraging a spatial filter

F

to smooth out different terms in the governing equations and eliminate spurious numerical oscillations that generally occur in the convection-dominated, under-resolved regime. In Section 4.1.1, we outline the main types of ROM filters and approximate deconvolution used to construct LES-ROMs. Then, we outline several LES-ROMs: the EFR-ROM (Section 4.1.2), which is arguably the most popular LES-ROM, Leray ROM (Section 4.1.3), approximate deconvolution Leray ROM (Section 4.1.3), and time-relaxation ROM (Section 4.1.3). Other types of LES-ROMs are mentioned in Section 4.1.4.

4.1.1. ROM Filters and Approximate Deconvolution

This section defines the ROM spatial filters and ROM approximate deconvolution used to construct the LES-ROMs presented in Section 4.1.2–Section 4.1.4. As discussed in Section 2, the spatial filter

F

is crucial in the development of LES for FOMs [4,5]. The spatial filter is also increasingly used to construct LES-ROMs for turbulent flows [175,176] because the idea is straightforward:

F

is used to eliminate the spurious numerical oscillations in stabilization strategies. Next, we outline the three spatial filters and the approximate deconvolution operators that are currently used for ROMs. In what follows, let

u

be a generic solution and

u_{r} = \sum_{j = 1}^{r} u_{r, j} φ_{j} (x)

its ROM approximation, with

φ_{j}

for

j = 1, 2, \dots, r

the ROM basis functions. We denote by

{\underset{̲}{u}}_{r} = {(u_{r, 1}, \dots, u_{r, r})}^{⊤}

the vector of ROM coefficients of

u_{r}

ROM Differential Filter

The ROM differential filter (DF) [175,177,178] is probably the most popular ROM spatial filter used to build LES-ROMs. The ROM differential filter is a modular, simple, linear filter, which is a natural extension of the differential filter used to construct the EFR-FOM in Section 2.1. It is defined as follows: find

{\bar{u}}_{r} (x) = \sum_{j = 1}^{r} {\bar{u}}_{r, j} φ_{j} (x)

such that

\begin{matrix} ({\bar{u}}_{r} - α^{2} Δ {\bar{u}}_{r}, φ_{i}) = (u_{r}, φ_{i}) \forall i = 1, \dots r, \end{matrix}

(59)

where

α

is the filter radius. The DF weak form (59) yields the following linear system:

\begin{matrix} (I + α^{2} A) {\underset{̲}{\bar{u}}}_{r} = {\underset{̲}{u}}_{r}, \end{matrix}

(60)

where

{\underset{̲}{\bar{u}}}_{r}

is the vector of ROM coefficients of

{\bar{u}}_{r}

and

I

and A are the identity and ROM stiffness matrices, respectively. We note that

I

in (60) is a consequence of the POD basis functions being orthonormal in the

L^{2}

norm. Thus, the ROM mass matrix is the identity matrix.

We emphasize that (60) is a low-dimensional,

r \times r

linear system, whose computational overhead is negligible (since r is relatively small). Thus, the ROM differential filter will be used to construct LES-ROMs that increase the ROM accuracy without significantly increasing the computational cost.

The differential filter was used in LES of turbulent flows with classical numerical discretizations [5,6,179]. In reduced-order modeling, the ROM DF was used to develop LES-ROMs for the Kuramoto–Sivashinsky equation [177], the NSEs [175,180], and the quasi-geostrophic equations [181].

ROM Higher-Order Algebraic Filter

The second filter we review is the ROM higher-order algebraic filter (HOAF): find

{\bar{u}}_{r} (x) = \sum_{j = 1}^{r} {\bar{u}}_{r, j} φ_{j} (x)

such that

\begin{matrix} (I + α^{2 m} A^{m}) {\underset{̲}{\bar{u}}}_{r} = {\underset{̲}{u}}_{r}, \end{matrix}

(61)

where m is a positive integer. As explained in [54] in the Fourier setting, the exponent m controls the percentage of filtering at different wavenumbers: as m increases, the amount of filtering increases for the high wavenumber components and decreases for the low wavenumber components. Just like the DF (60), the HOAF (61) is a relatively low-dimensional,

r \times r

linear system that can be efficiently solved. Thus, the HOAF also leads to the development of accurate and efficient LES-ROMs. We also note that for

m = 1

, the linear systems (61) and (60) are identical. Thus, DF can be considered a particular case of HOAF with

m = 1

The HOAF (61) was proposed in [176] and was based on the HOAF introduced by Fischer and Mullen [54] in a spectral element method (SEM) setting. In [176], the HOAF (61) was called the higher-order differential filter, to be consistent with the SEM nomenclature. In [182], we showed that the HOAF is related to, but slightly different from, the spatial discretization of a higher-order differential operator (they are the same in the periodic case). For clarity, in this paper, we call the operator in (61) the high-order algebraic filter. In [182], a theoretical and numerical investigation of HOAF is presented.

ROM Projection

The third ROM spatial filter used to construct LES-ROMs is the ROM projection filter, which is defined as follows: for a fixed

r_{1} < r

and a given

u_{r} \in X^{r} = span {φ_{1}, \dots, φ_{r}}

, the ROM projection (Proj) seeks

{\bar{u}}_{r} \in X^{r_{1}} = span {φ_{1}, \dots, φ_{r_{1}}}

such that

\begin{matrix} ({\bar{u}}_{r}, φ_{j}) = (u_{r}, φ_{j}) \forall j = 1, \dots, r_{1} . \end{matrix}

(62)

To our knowledge, the first time the Proj was used as an explicit ROM spatial filter was in [167], where it was leveraged to construct the ROM dynamic subgrid scale (SGS) model. We note that the dynamic SGS is generally considered to be the state-of-the-art closure model in LES [4,5]. Next, the ROM projection was used in [175,178] to construct the Leray ROM (see Section 4.1.3) and EFR-ROM (see Section 4.1.2).

ROM Filter Radius

The filter radius

α

is fundamental in the definitions of both the ROM differential filter (59) and the ROM higher-order algebraic filter (61): the radius of the ROM spatial filters represents the size of the scales that are captured by the filtered variables. Furthermore, the filter radius is also essential in the approximate deconvolution process, which depends on the particular filter used. Thus, a natural question is “How do we choose the ROM filter radius?” For example, the following practical question could be asked:

\begin{matrix} \begin{matrix} Given 10 ROM basis functions, what is the resolved lengthscale ? \end{matrix} \end{matrix}

(63)

The most straightforward answer to question (63) is to choose the ROM filter radius the same way we choose the FOM filter radius. As mentioned in Section 2, the first choice is a computational scale, i.e., a scale based on the underlying numerical discretization. We also mention that other

α (h)

scalings are possible [5,6]. Another choice for the FOM differential filter radius is a physical scale, i.e., based on the physics of the underlying flow. For example, a popular choice in this class is

α

∼

η

, where we recall

η

is the Kolmogorov scale.

One limitation of using the computational scale or the physical scale to define the ROM filter radius is that these scales lack adaptivity, i.e., do not change when we vary the ROM dimension, r. To address this limitation, in [183] we put forth a new energy-based ROM lengthscale, which was defined by using a judicious partition of the energy. We showed that the new energy-based ROM lengthscale displays the correct asymptotic behavior: it converges to the smallest lengthscale in the system when the ROM dimension increases, and it converges to the largest lengthscale in the system when the ROM dimension decreases. Furthermore, in the numerical simulation of the turbulent channel flow at Reynolds number

R e_{τ} = 395

, we compared the EFR-ROM equipped with the new energy-based ROM lengthscale and a classical ROM lengthscale based on dimensional arguments [184]. The numerical results showed that the former is more robust with respect to parameter changes than the latter. Finally, we note that the new energy-based ROM lengthscale could also be useful in the preprocessing strategy advocated in [185,186] as a means to filter out the noise in the input data.

Approximate Deconvolution

Although the AD strategy is central in image processing and inverse problems [187,188] and it has been used in LES development (see Section 2.1.3 and the research monographs [4,5,6,29,30]), it is a relatively new concept in reduced-order modeling. The AD strategy can be formulated as follows: given an approximation of the filtered input variable,

\bar{u} ≐ F u

, where

F

is an invertible spatial filter, find an approximation of the unfiltered input signal,

u

. Of course, the exact deconvolution,

F^{- 1} \bar{u}

, would seem a natural choice. Computationally, however, the exact deconvolution is a very bad idea. Indeed, as shown in [187,188], the noise in the high wavenumber components of

\bar{u}

is amplified by the inverse filter,

F^{- 1}

, in the exact deconvolution. Thus, in practice, AD strategies are used instead of the exact deconvolution [187,188]. In LES of turbulent flows, the AD models were pioneered by Adams and Stolz for classical numerical discretizations [87]. To our knowledge, in reduced-order modeling, the first AD model was proposed in [189], where AD was used to develop a ROM closure model. Recently, the AD strategy was also used to develop the AD Leray ROM [190], where three different AD strategies were investigated: the van Cittert AD approach (Section 3.2 in [190]), the Tikhonov AD approach (Section 3.3 in [190]), and the Lavrentiev AD approach (Section 3.4 in [190]). We also note that the Lavrentiev AD strategy was used in [189] to develop a ROM closure model, and the Tikhonov approach was used in [191,192,193] for ROM regularizations. These three AD methods were investigated numerically in Section 5.2 in [190].

4.1.2. EFR-ROM

The most popular LES-ROM exports the EFR strategy (11)–(13) to the reduced-order level, in which the FOM variables

U

and

V

are replaced by the ROM variables

u_{r}

and

v_{r}

, respectively, and the filter and deconvolution operators are replaced by their ROM counterparts. This yields the evolve–filter–relax ROM (EFR-ROM):

-: Evolve: find $v^{n + 1}$ such that

$\begin{matrix} \frac{v_{r}^{n + 1} - u_{r}^{n}}{Δ t} + M (v_{r}^{n + 1}) = 0 . \end{matrix}$

(64)
-: Filter: find filtered variable ${\bar{v}}_{r}^{n + 1}$ such that

$\begin{matrix} {\bar{v_{r}}}^{n + 1} = F v_{r}^{n + 1} . \end{matrix}$

(65)
-: Relax: set

$\begin{matrix} u_{r}^{n + 1} = (1 - χ) {v_{r}}^{n + 1} + χ {\bar{v_{r}}}^{n + 1}, \end{matrix}$

(66)

with relaxation parameter $χ \in [0, 1]$ .

In the evolve step of the EFR-ROM algorithm, one time step of the standard Galerkin projection ROM (G-ROM) time discretization is leveraged to advance the EFR-ROM approximation at the current time step,

u_{r}^{n}

, to an intermediate EFR-ROM approximation,

v_{r}^{n + 1}

. In the filter step of the EFR-ROM algorithm, one of the ROM spatial filters presented in Section 4.1.1 is used to filter the intermediate EFR-ROM approximation from the evolve step and obtain a smoother ROM approximation without spurious numerical oscillations. Finally, in the relax step of the EFR-ROM algorithm, the EFR-ROM approximation at the next time step is defined as a convex combination of the unfiltered intermediate EFR-ROM approximation obtained in the evolve step,

v_{r}^{n + 1}

, and its filtered counterpart,

{\bar{v}}_{r}^{n + 1}

. Just like in the case of EFR used as FOM (see Section 2.1), the goal of the relax step is to adjust the amount of dissipation introduced in the filter step by using a relaxation parameter,

0 \leq χ \leq 1

. By varying the relaxation parameter

χ

, one can produce a range of filter strengths, from no filtering at all (

χ = 0

) to maximum filtering (

χ = 1

). As pointed out in Section 3.1.1 for the EFR used at the FOM level, the choice of

χ

may be critical in the presence of backflows. At the ROM level, the numerical investigation in [180] has shown that EFR-ROM is sensitive with respect to

χ

The evolve–filter ROM was introduced in [175] and EFR-ROM in [176]. Since then, EFR-ROM has been developed in several directions, for example, numerical simulation of turbulent flows (e.g., the turbulent channel flow) [183] and ocean dynamics [181], the FOM-ROM consistency [180], and feedback control [194].

4.1.3. Leray ROM, Approximate Deconvolution Leray ROM, and Time-Relaxation ROM

In this section, we present several LES-ROMs constructed by using ROM spatial filtering and approximate deconvolution.

Leray ROM

The Leray ROM (L-ROM) [175,195,196] is an LES-ROM in which the spatial filter,

F

, is used to smooth out only the convective velocity of the nonlinear term in the NSE, i.e.,

\begin{matrix} (\frac{\partial u_{r}}{\partial t}, φ_{i}) + R e^{- 1} (\nabla u_{r}, \nabla φ_{i}) + (({\bar{u}}_{r} \cdot \nabla) u_{r}, φ_{i}) = 0, \forall i = 1, \dots, r, \end{matrix}

(67)

where

{\bar{u}}_{r}

is the ROM velocity filtered with one of the ROM spatial filters, i.e., DF (60), HOAF (61), or Proj (62). By introducing the spatial filtering of the convective term, the Leray ROM smooths out the spurious numerical oscillations that affect the standard G-ROM in the convection-dominated, under-resolved regime.

Leray regularization (see Section 2.1) was first used in the context of reduced-order models in [177] for the Kuramoto–Sivashinsky equations. For fluid flows, L-ROM was first used in [175] for the 3D flow past a circular cylinder at

R e = 1000

. Since then, L-ROM has been successfully used as a stabilization technique for various under-resolved flows: the NSEs [71,197], the stochastic NSEs [176,196], the Boussinesq equations [198,199], the quasi-geostrophic equations [66,181], and the turbulent channel flow [182]. We also note that the ROM projection was leveraged in [198,199,200] to build the Leray ROM for challenging applications.

Approximate Deconvolution Leray ROM

The approximate deconvolution Leray ROM (ADL-ROM) [190] is an extension of the Leray ROM that leverages the AD operator

D

. The idea of the ADL-ROM is to use AD to limit the amount of dissipation introduced in the Leray ROM when the filter radius is too large: Given

u_{r}^{n}

, find

u_{r}^{n + 1}

such that

(\frac{u_{r}^{n + 1} - u_{r}^{n}}{Δ t}, φ_{i}) + R e^{- 1} (\nabla u_{r}^{n + 1}, \nabla φ_{i}) + ((D ({\bar{u}}_{r}^{n + 1}) \cdot \nabla) u_{r}^{n + 1}, φ_{i}) = 0, \forall i = 1, \dots, r .

(68)

The numerical investigation in [190] for convection-dominated systems shows that when the filter radius is relatively large, the ADL-ROM is more accurate than the standard Leray ROM and less sensitive to model parameters. In addition, it shows that the increased accuracy allowed by the ADL-ROM (68) is not at the expense of numerical stability.

We emphasize that this is yet another instance in which ideas from completely different fields are synthesized. Indeed, AD is central in the image processing and inverse problems communities and has been successfully used to develop LES models [87,88], as mentioned also in Section 2.1.3.

Time-Relaxation ROM

The time-relaxation ROM [182] is an LES-ROM recently proposed with the aim to increase the ROM numerical stability by adding an extra term to the incompressible NSE. Specifically, the time-relaxation ROM reads as follows: Find

u_{r}

such that

\begin{matrix} (\frac{\partial u_{r}}{\partial t}, φ_{i}) + R e^{- 1} (\nabla u_{r}, \nabla φ_{i}) + ((u_{r} \cdot \nabla) u_{r}, φ_{i}) + (χ_{t} (u_{r} - {\bar{u}}_{r}), φ_{i}) = 0, \forall i = 1, \dots r, \end{matrix}

(69)

where

χ_{t}

is the time-relaxation parameter and

{\bar{u}}_{r}

is the ROM velocity filtered with one of the ROM spatial filters introduced in Section 4.1.1. The role of this extra term is to add numerical stabilization only at the marginally resolved scales without unnecessarily filtering the resolved scales.

In [182], it was shown that the numerical simulation of the turbulent channel flow at Reynolds numbers

R e_{τ} = 180

and

R e_{τ} = 395

performed with the new time-relaxation ROM equipped with the ROM differential filter and the ROM higher-order algebraic filter yields accurate results.

4.1.4. Other LES-ROMs

There are numerous LES-inspired ROMs, and the list is continuously growing. For comprehensive reviews on ROM closures, the reader is referred to [169,201]. We note, however, that many of the ROM closures presented in [169,201] do not use ROM filtering or approximate deconvolution. Examples include the eddy viscosity ROMs and machine learning ROM closures [169,201]. Since the main goal of this paper is to discuss a specific class of LES-inspired ROMs, which are constructed by using ROM spatial filtering and approximate deconvolution, we do not discuss these alternative ROM closures and instead refer the reader to [169,201]. In this section, however, we outline two types of LES-ROMs that are different from those we presented in Section 4.1.2–Section 4.1.3.

The approximate deconvolution ROM (AD-ROM) was introduced in [189]. To our knowledge, this is the first AD model used in reduced-order modeling. The AD-ROM is constructed by using AD (i.e., a Lavrentiev regularization) to construct an accurate ROM closure model. In [189], the AD-ROM was successfully tested in the numerical simulation of the three-dimensional flow past a circular cylinder at a Reynolds number

R e = 1000

The variational multiscale (VMS) methods [202,203] are constructed by using the principle of the locality of energy transfer, which states that energy is transferred mainly between neighboring scales or modes. Since ROMs use hierarchical bases in which the large and small structures are clearly displayed, the VMS framework naturally lends itself to be extended to the ROM setting. Although VMS-ROMs do not explicitly use a ROM spatial filter with an explicit filter radius, they display a connection with the LES-ROMs presented in this paper: In the VMS-ROMs, the decomposition of the flow variables into large scales and small scales can be regarded as the result of using the ROM projection operator [173,204,205]. For more discussion on VMS-ROMs, we refer the reader to [169] (Section IV.A.5).

4.2. LES-ROM Consistency

As explained in [180,206], ROMs are of two types:

FOM-ROM consistent, i.e., the ROM uses the same computational model and the same numerical discretization as the FOM (see [206] Definition 1.1).
FOM-ROM inconsistent, i.e., the ROM uses a computational and/or numerical discretization that are different from those used by the FOM.

We note that most ROMs (and data-driven models in general) are FOM-ROM inconsistent. The reason is simple: the FOM-ROM inconsistent ROMs (which are often called “nonintrusive”) are much easier to use in practice, e.g., because they allow for the use of legacy codes at the FOM level. We emphasize, however, that recent numerical investigations [180,206] have shown that the FOM-ROM consistent ROMs yield more accurate results than the FOM-ROM inconsistent ROMs. For example, in [180] we generated the FOM data by using the EFR-FOM and we considered two types of ROMs: a FOM-ROM consistent ROM that uses the EFR strategy at the ROM level and a FOM-ROM inconsistent ROM that does not use the EFR strategy at the ROM level. Our numerical investigation showed that the FOM-ROM consistent ROM was more accurate than the FOM-ROM inconsistent ROM. Other examples of FOM-ROM consistent ROMs based on the EFR strategy are presented in [71,197] (see [207,208,209] for FOM-ROM consistent ROMs based on different stabilizations).

The LES-ROMs that we present in this section are FOM-ROM consistent since they use an LES strategy both at the FOM and the ROM levels. Thus, based on the preliminary studies in [180] for the EFR method, we believe that LES-ROMs will yield more accurate results than their inconsistent counterparts. We also note that by ensuring the modeling consistency between the FOM and the ROM, we increase the ROM robustness with respect to changes in the FOM parameters. For example, when the FOM resolution needs to be adapted by changing the FOM filter radius, the ROM filter radius is adapted automatically and the ROM does not have to be retrained. We emphasize that this is a departure from most of the current ROM closures and stabilizations, which require retraining (a costly process) every time an LES model parameter is changed at the FOM level.

Of course, as pointed out in Section 5 of [180], there are still many open questions. For example, one could investigate the FOM-ROM consistency when different LES models are used at the FOM and ROM levels, e.g., the Leray model is used at the FOM level and the EFR model is used at the ROM level. One could also investigate the parameter FOM-ROM consistency, which is complementary to the model FOM-ROM consistency. To investigate the parameter FOM-ROM consistency, one could consider the same LES model at the FOM and ROM levels but use different parameters (e.g., different

α

values) in these LES models. For example, one could investigate whether using different

α

χ

values at the FOM and ROM levels (i.e.,

α^{F O M} \neq α^{R O M}

χ^{F O M} \neq χ^{R O M}

) could yield more accurate ROM solutions.

We also emphasize that the FOM-ROM consistency has been investigated mainly numerically, without providing theoretical support. For notable exceptions, see the numerical analysis performed in [207,210] for FOM-ROM consistency of the streamline upwind Petrov–Galerkin (SUPG) stabilization and [206] for FOM-ROM consistency with respect to the discretization of the nonlinearity of the Navier–Stokes equations. We note, however, that the numerical analysis of the FOM-ROM consistency for LES-ROMs is still an open question.

4.3. LES-ROM Numerical Analysis

At the FOM level, extensive mathematical support exists for LES closures and stabilizations with classical numerical discretizations (e.g., the finite element method). For example, the monographs [5,29,30] present the mathematical analysis for many LES models, as well as the numerical analysis of their discretization. Similarly, the monograph [211] presents the mathematical and numerical analysis of classical stabilization strategies. This extensive literature provides answers not only to fundamental numerical analysis questions (e.g., stability and convergence) but also to critical practical challenges (e.g., parameter scalings for critical parameters like the filtering radius). In stark contrast, fundamental questions in the numerical analysis are still wide open for most of the ROM closure models: Is the proposed ROM stable? Does the ROM converge? If so, what does it converge to? What is the ROM rate of convergence?

Only the first steps in the numerical analysis of ROM closures and stabilizations have been taken: numerical analysis (e.g., stability and convergence) of classic ROM closures was performed in [212,213,214,215,216]. Numerical analysis was carried out in [213,214] for the eddy viscosity variational multiscale ROMs and in [215,216] for the Smagorinsky model in a reduced basis method setting. The first numerical analysis of LES-ROMs was performed in [195,196] for the Leray ROM (67) (see also [217] for related work) and in [194] for the EFR-ROM. We have also taken the first steps in proving FOM-ROM parameter scalings for ROM stabilizations and closures [207,213,214,218]. For example, in [207], for the SUPG-ROM, we proved scalings between the FOM mesh size and the SUPG-ROM stabilization parameter.

Finally, we note that the most active research area in ROM closure modeling is in the development of data-driven ROM closures in which available data are utilized to build the ROM closure model. An example of data-driven ROM closure is the data-driven variational multiscale ROM (d2-VMS-ROM) proposed in [205,219,220]. The d2-VMS-ROM and its developments have been investigated numerically in [205,219,220,221,222,223,224,225,226]. However, providing mathematical support for the d2-VMS-ROM, and data-driven ROM closures in general, is an open problem. Recently, in [227], we laid the mathematical foundations of data-driven ROM closures. Specifically, in Theorem 2 of [227], we proved that the d2-VMS-ROM is verifiable, i.e., the d2-VMS-ROM solution is accurate since the d2-VMS-ROM closure model is accurate.

Despite these first steps that our groups and our collaborators have taken, the numerical analysis of ROM closures and stabilizations is nowhere close to the numerical analysis for FOM closures and stabilizations. Laying the mathematical foundations for LES-ROMs and, more importantly, providing numerical analysis guidance to practical LES-ROM choices (e.g., parameter scalings) are still open questions.

5. Applications of LES for ROM

In this section, we present some results of LES with ROMs for test cases with an academic flavor and some real applications. We also draw a road map of other potential applications for more real-world problems that motivated us in the development of the methodologies described above. We anticipate these applications will provide an excellent benchmark for the assessment and improvement of such methodologies. To illustrate the versatility of the LES-ROMs, we show results for both incompressible flows (Section 5.1) and compressible flows (Section 5.2), as we did in Section 3.

5.1. Incompressible Flows

5.1.1. Flow Past a Cylinder

Although simple (or maybe precisely because of that), the 2D [228,229] and 3D [228,230,231] flow past a cylinder benchmarks have been widely used to assess numerical strategies, even those designed for higher Reynolds number flows. This is due to the fact that to obtain accurate results in the case of the time-dependent Reynolds number

0 \leq R e (t) \leq 100

, one needs a mesh with about 200 K elements for the 2D benchmark and roughly 3 million elements for the 3D benchmark. Thus, one can use these benchmarks to demonstrate that accurate results can be obtained also with much coarser meshes if suitable LES models are used. Indeed, for example, it has been shown that with the EFR algorithm, one can obtain accurate results with meshes as coarse as roughly 15 K elements in both 2D and 3D [58,59,61,62,71,180,232]. Below, we present selected 2D results from [197], which touch upon several topics discussed in Section 4.1.

The work in [71,197] stemmed from the need for FOM-ROM consistency described in Section 4.2. We proposed to use the EFR algorithm at both the full- and reduced-order level. In particular, we generated the snapshots with under-refined meshes. This is unlike other works, e.g., [175,176,233], where the snapshots are obtained by DNS. By using the EFR method as FOM and ROM, we adopted the same mathematical framework during both the offline and online stages and thus had a ROM that is fully consistent with the FOM. In [197], we proposed to use the POD basis related to the evolve velocity to approximate the filtered velocity and compute the reduced pressure field with a Poisson pressure equation method [234,235]. The main difference between [197] and [71,180,232] lies in the indicator function: it is linear in [71,180,232] and AD-based (52) in [197]. The ROM in [197] is called hybrid projection/data-driven because it exploits a traditional projection method (G-ROM) for the computation of the intermediate reduced velocity and pressure fields and uses a data-driven technique to compute the reduced coefficients of the indicator function field. The data-driven approach leverages interpolation with radial basis functions [236]. Alternatives to the data-driven approach could be using the same set of reduced coefficients for velocity, pressure, and indicator function or an EIM/DEIM technique [237,238]. The former was considered for RANS in [239] and was shown in [240] to provide less accurate results than a hybrid procedure.

The 2D flow past a cylinder for

0 \leq R e (t) \leq 100

is challenging because the flow is laminar for the first 4 s and then the well-known vortex shedding appears. In [197], the focus was on the time interval

[4, 8]

s, for which 200 high-fidelity snapshots were collected. The first 50 most energetic POD modes were considered, and a convergence test was performed based on three different energy thresholds: 99% (11 modes for the velocity, 5 modes for the pressure, and 20 modes for the indicator function), 99.9% (26 modes for the velocity, 12 modes for the pressure, and 44 modes for the indicator function), and 99.99% (42 modes for the velocity, 24 modes for the pressure, and 50 modes for the indicator function). Figure 11 compares the lift coefficient given by the FOM with the one given by the ROM. We observe that the evolution of lift coefficient computed by the ROM is very accurate when 99.9% or 99.99% of the eigenvalue energy is retained. Similar accuracy was found for the 3D cylinder test. See [197] for more details.

To complete the picture, we need to connect the accuracy with the required computational time. Figure 12 (left) shows the time-averaged relative

L^{2}

error for velocity and pressure versus the corresponding wall time when the number of basis functions for the velocity is varied. The number of basis functions for the pressure is fixed at 24, while the number of basis functions for the indicator function is 50. The relative wall time is given by the ratio between the CPU time needed by the ROM and the CPU time required by a FOM simulation (i.e., 398 s). As one would expect, the errors initially decrease and the relative wall time increases as the number of basis functions increases. The errors reach their minimum value for 46 basis functions for the velocity. When the number of basis functions is increased to 50, two undesirable things happen: the relative wall time is larger than 1, meaning the ROM is more costly than the FOM (an absurdity), and the errors increase. Since 42 velocity modes are enough to retain 99.99% of the eigenvalue energy, we suspect the errors increase because of the noise introduced by the highest modes. We note that the ROM is more computationally efficient in the 3D test. See Figure 12 (right) and [197] for details.

5.1.2. T-Junction

Here, we summarize the results of the numerical investigation performed in [200]. We investigate three of the LES-ROMs discussed in Section 4.1 (i.e., the EFR-ROM, Leray ROM, and time-relaxation ROM) in the numerical simulation of a T-junction problem. When streams of rapidly moving flow merge in a T-junction, large oscillations can occur at the scale of the diameter, D, with a period of

O (D / U)

, where U is the characteristic flow speed. If the streams are at different temperatures, the oscillations result in fluctuations (see Figure 13). This phenomenon, known as thermal striping, can accelerate thermal–mechanical fatigue at the pipe wall in the outlet branch and ultimately cause pipe failure. Since thermal striping is of critical importance in nuclear engineering, the nuclear energy modeling and simulation community established a T-junction benchmark [241] to test the ability of CFD codes to predict thermal striping.

We consider the T-junction problem at the Reynolds number

R e

= 10,000. The FOM is based on a spectral element discretization with about 21 million mesh points. Further details about the computational setting and FOM are given in [200]. We test the standard G-ROM and the three LES-ROMs in the predictive regime, that is, we test these models on a time interval that is larger than the interval where snapshots are collected. We consider

r = 250

for the G-ROM and

r = 100

for the LES-ROMs. To assess the ROMs’ performance, we compare the near-wall temperature at several locations (

(x, \pm 0.5, 0), (x, \pm 0.45, 0), (x, 0, \pm 0.5)

, and

(x, 0, \pm 0.45)

for

x = 2, 3, . . ., 9

) with the FOM data (see Figure 14).

Figure 15 and Figure 16 show the near-wall (i.e.,

z = 0.45

) temperature of the downstream pipe at several x locations for the G-ROM, the Leray ROM, the EFR-ROM, and the time-relaxation ROM. In these plots, the two vertical red lines denote the time interval in which the snapshots were collected. The results indicate that the near-wall temperature of the G-ROM initially agrees with the FOM but becomes unstable in a relatively short time (≈6 convective time units after the initial condition) even with

r = 250

. This indicates that the G-ROM is not accurate even in the reconstructive regime. Furthermore, larger fluctuations are observed for small x values. In contrast, the LES-ROM temperature stays on the overall FOM trajectory, and in the predictive regime (i.e., outside the region between the two vertical red lines) all three LES-ROMs are significantly more accurate than the G-ROM.

5.1.3. Hemodynamics Applications

As mentioned in Section 3.1, blood flow is generally in a laminar regime, and it does not require specific modeling for turbulence. AoDs are an exception where the combination of the specific flow features induced by the tiny entry tears and the need for an efficient solver to retrieve information from many patients of a (computer-assisted) clinical trial call for specific modeling solutions like EFR.

There are other problems in the field of computational hemodynamics where the need for extremely efficient solvers is critical. Typically, these are problems where the use of specific medical devices requires an optimization procedure. While “optimization” can be performed by trial-and-error approaches, translational mathematicians (i.e., mathematicians engaged in bringing state-of-the-art mathematical and computational procedures to the forefront of medical practice) feel the need for introducing automatic or, at least, semi-automatic rigorous, mathematically sound procedures. The cost of several attempts to find the optimal solution, either as the result of empirically educated guesses or iterations of a mathematical optimization procedure, may be prohibitive. In general, in this kind of application, efficiency is privileged with detrimental consequences for accuracy. For this reason, RANS models are often preferred.

An example of great interest in pediatric surgery is related to the so-called total cavopulmonary connection (TCPC). This is a palliative surgery for newborn babies with a severe heart malformation called left ventricle hypoplasia syndrome. This condition (called univentricular circulation) is not compatible with life and requires specific treatments. Years ago, the TCPC was introduced to buy time before a transplant [242]. It consists of an artificial connection between the superior and inferior venae cavae (SVC and IVC, respectively) with the pulmonary artery (PA) to route systemic circulation to the pulmonary circulation. An oversimplified representation of this shape is reported in Figure 17. The functioning ventricle pumps blood into the large circulation. The connection between the venae cavae and the PA has a peculiar cruciform shape, shown in Figure 17, with many interesting implications in terms of fluid dynamics. The biomedical engineering community has investigated these implications in a huge number of excellent contributions, see, e.g., [243,244,245,246,247,248,249,250,251,252]. From the computational point of view, the presence of the colliding fronts injected at the SVC and IVC raises some challenges since even moderate Reynolds numbers may induce flow disturbances. In general, the shape of the cruciform connection was speculated to have a major role in the long-term healthy conditions of the patient. For instance, an even splitting of the flow coming from the IVC (called hepatic flow) is critical and is determined by the offset between the SVC and IVC left after the surgery (see Figure 17). However, a rigorous shape optimization based on mathematical procedures in the operating room is still a dream for many reasons. One is the identification of all the relevant criteria that define a reliable and robust “optimal” solution. Another one is the huge computational cost of numerical optimization. On top of this, a more recent variant of the surgery requires the introduction of a pump in the cross-shaped junction to establish hemodynamic conditions closer to the physiology of biventricular circulation [253]. The regime induced by the pump requires specific turbulence modeling [7]. High computational costs currently lead to the use of RANS models, even though the superiority of LES modeling for the flow regimes characteristic of TCPC is commonly recognized [254]. Our intention is to bring LES-ROM modeling into this field as a potential breakthrough that enables not only accurate LES modeling as the routine approach but also rigorous shape optimization procedures as part of regular surgical planning.

5.1.4. Wind Energy Applications

LES-based wind modeling presents significant challenges, with key questions posed on turbulence, mostly focused on mean-flow kinetic energy (MKE) entrainment and dissipation in large wind farms [255]. Understanding the upper limit of power production is linked to MKE entrainment. Recent studies, such as [256], identified spatial constraints and the influence of the Coriolis parameter on power output through idealized atmospheric simulations run with the weather research and forecasting model. An entrainment-based model for wind farm flow and power prediction was recently proposed in [257]. Wind energy, projected to meet 35% of global electricity needs by 2050, faces significant challenges, including optimizing wind farm layouts and minimizing wake interactions that reduce power output and affect turbine longevity. With this in mind, it is worth noting that enhancing energy production by just 1% could yield 30 TWh annually, equivalent to adding 3600 turbines and USD 1 billion in revenue [258].

The filtering-based approaches reviewed in this paper can provide accurate modeling tools to predict wind farm wakes, addressing spatio-temporal variability, unsteady wake interactions, and atmospheric turbulence. Moreover, improved turbine wake prediction and cost-effective surrogate models can optimize turbine layout and support wind farm deployment [259]. For further insights on LES of wind applications, we refer the reader to references [256,257,260,261,262,263,264]. Moreover, a recent study [265] offers a comprehensive overview of digital twin technology, with a focus on wind energy applications. The study consolidates digital twin definitions, identifies the current state of the art in modeling and simulation techniques, and outlines research needs in the wind energy sector from an industrial perspective. It proposes solutions to these challenges from research institutes’ viewpoints and provides recommendations for stakeholders to facilitate technology adoption, with the methodologies highlighted in our article seen as key enablers for this digital technology adaptation.

5.2. Compressible Flows

ROMs combined with LES modeling has the potential to be a major player in reducing the computational cost associated with weather forecasting. This section is meant to convey an ideal of such potential and show some limitations of current methodologies as an incentive to undertake more work in this research area. Mainly, we report results from [266] obtained with selected data-driven ROMs for the same benchmarks for atmospheric flow described in Section 3.2. However, let us start with a brief literature review.

While POD, which is also referred to as empirical orthogonal function (EOF) analysis in the geophysical fluid dynamics community, has been used for a long time, it is only recently that ROMs have been applied for the simulation of atmospheric flows. EOF has been applied to identify spatio-temporal coherent meteorological patterns, e.g., the Madden–Julian oscillation, the quasi-biennial oscillation, and the El Niño–southern oscillation. See, e.g., [267,268,269,270,271]. The EOF analysis for these phenomena uses data on the global scale and considers time (ranging from several months to many years) as the only parameter. In addition, it is mostly limited to data analysis as a mean to understand the weather system, i.e., EOF is not used for forecasts. Only very recently, a data-driven ROM based on EOF analysis has been used for pattern prediction, specifically to forecast the weekly average sea surface temperature [268]. Other data-driven methods borrowed from machine learning have been applied to global weather forecasting. See, e.g., [272,273,274,275,276,277,278].

The work in [266] is a first attempt to apply ROMs that generate a reduced-order basis (unlike machine learning methods) to both reconstruct and forecast regional atmospheric flows (unlike EOF, which is not used to forecast). Additionally, since the spatial scale is a few kilometers and the time scale is a few hours, there is an obvious difference in resolution from EOF analysis. The considered data-driven ROMs are dynamic mode decomposition (DMD), Hankel dynamic mode decomposition (HDMD), and proper orthogonal decomposition with interpolation (PODI). DMD was specifically designed to predict the future behavior of a system [279,280,281,282] and has been successfully applied to several incompressible fluid dynamics problems [283,284,285]. HDMD enhances the DMD algorithm with time-delay embedding [286,287,288,289,290,291] and is able to predict more accurately and for longer periods of time systems exhibiting strong nonlinear dynamics [290,292]. Successful applications include periodic cavity flow [286,293], electromechanical systems [291], and biological systems [288]. Unlike DMD and HDMD, PODI was not designed to forecast system evolution but rather to interpolate solutions in a parameter space, where time is one of possibly many parameters of interest. To date, PODI has been applied to perform parametric studies for problems in hemodynamics [294], chemical [295] and naval [296,297] engineering, and aeronautics [298].

For the rising thermal bubble, we collect a database of potential temperature perturbations consisting of 204 snapshots, i.e., the computed

θ^{'}

every 5 s. For simplicity, we used a constant eddy viscosity

ν_{a} = 15

m²/s, which is frequently used in the literature (see, e.g., [43,151,159]). Then, 90% of the database (i.e., 184 snapshots) is used for training. In the case of DMD and HDMD, these 184 solutions are the first 184 in the database (associated with the time interval

(0, 920]

s). In the case of PODI, these 184 solutions are selected randomly over the entire time interval

[0, 1020]

s. The remaining 20 solutions form the validation set. This difference in the training and validation sets reflects the different nature of PODI and DMD/HDMD algorithms. To retain 99% of the eigenvalue energy in the database, one needs 17 POD modes for DMD and PODI and 46 modes for HDMD. Figure 18 shows a qualitative comparison of the ROM solutions with the FOM solution. Among the four time instants chosen for the visualization, two of them (

t = 255

and

t = 505

s) correspond to solutions belonging to the training set and allow us to asses the ability of each ROM technique to identify the system dynamics. We see that all three ROMs perform well in system identification. The remaining two times (

t = 980

and

t = 1020

s) are not associated with the training set and thus are used to check the accuracy of the ROM in predicting (for DMD and HDMD) or interpolating (for PODI) the system dynamics. For these times, DMD provides a poor approximation of

θ^{'}

, while HDMD and PODI reconstruct the solution well. The problem with the reconstruction provided by DMD appears to be related to assigning large weights to basis functions associated with previous times. Indeed, in the DMD solution for

t = 1020

s, we can observe the time history of the system dynamics, i.e., the rising of the bubble.

All the simulations whose results are reported in Figure 18 were run on a common laptop. Table 3 reports the computational time needed to construct the reduced basis offline and to perform a simulation online for each of the ROMs we considered. The computational times for DMD and PODI are comparable, which for the online phase is explained by the fact that the DMD and PODI reduced basis have the same size, while HDMD is computationally more expensive.

The density current benchmark proved to be more challenging for the fact that it features mainly vertical dynamics for part of the time interval (bubble descending) and mainly horizontal dynamics for the rest of the time (after the bubble hits the ground and a cold front propagates). To separate the two dynamics [180,197], we collect snapshots associated only with the cold front propagation, i.e., the computed

θ^{'}

every 4 s for

t \geq 280

. Furthermore, in this case, we set the eddy viscosity to a constant value (

μ_{a} = 65

m²/s) taken from the literature [151,154]. Out of the snapshot database,

85 %

is used for training and

t = [808, 900]

s is used as the prediction window. We set again the eigenvalue energy threshold to 99%, which leads to retaining 49 POD modes for DMD, 96 modes for HDMD, and 52 modes for PODI. Figure 19 compares the evolution of the potential temperature perturbation given by DMD, HDMD, and PODI with the evolution computed by the FOM. The top two rows in Figure 19 correspond to solutions belonging to the training set. Once again, we see that all three ROMs perform well in system identification. The bottom two rows in Figure 19 correspond to solutions not within the training set. The DMD solutions are affected by the same problem observed in Figure 18: we can see the evolution of

θ^{'}

in each panel because large weights are assigned to basis functions associated with previous times. The HDMD solutions are a clear improvement, but the PODI solutions compare more favorably with the FOM solutions.

In conclusion, although DMD and HDMD are intended for forecasts, the accuracy in the prediction of the system dynamics is low even when 99% of the eigenvalue energy is retained and the snapshots in the database are tailored to the problem at hand. Thanks to the interpolatory approach, PODI maintains a good level of accuracy during the entire time interval of interest. In [266], it is shown that this is true also when a physical parameter (not just time) is varied within a parametric study and for 3D results. While still preliminary, these results show that a lot more work is needed if we want to use ROMs to cut the computational time required to forecast the weather for long prediction windows and in high-dimensional parameter space. The results for the T-junction benchmark test summarized in Section 5.1 suggest that L-ROM, EFR-ROM, and TR-ROM could likely improve the results presented in this section in terms of accuracy. Additional improvements in terms of efficiency could come from machine learning-based techniques that can better detect and reproduce the nonlinear behavior exhibited by the FOM, e.g., convolutional autoencoders and long short-term memory [299,300,301,302].

6. Concluding Remarks and Outlook

In this paper, we connect two important research fields that have been treated separately until now, LES and ROMs, by reviewing ROMs that are inspired from LES strategies, which we call LES-ROMs. In our review, we tried to emphasize two essential features of LES-ROMs: (i) The LES-ROM achievements in reduced-order modeling of realistic, convection-dominated, under-resolved flows; (ii) Their rigorous, mathematical foundations, which enable a clear physical interpretation of the results.

Convection-dominated flows in the under-resolved regime are central in important engineering, scientific, and medical applications. Standard ROMs generally fail to produce acceptable results in these challenging settings. Thus, alternatives are urgently needed. In our review, we emphasized the significant achievements of LES-ROMs in the efficient and accurate numerical simulation of convection-dominated, under-resolved flows in important applications in engineering (e.g., turbulent flows in aerospace industry), science (e.g., weather/climate modeling), and medicine (e.g., cardiovascular flows). In these critical applications, LES-ROMs have yielded accurate results, while their cost was orders of magnitude lower than the cost of traditional numerical methods (e.g., the finite element method), also called FOMs. We also emphasized that an important advantage of LES-ROMs is that they often are extremely easy to implement. Indeed, one can start with a “legacy” ROM code (e.g., one of the current ROM software libraries [303,304]) and, in a matter of minutes, implement an efficient and accurate LES-ROM. We illustrated the significant achievements of LES-ROMs for both incompressible flows (e.g., the incompressible Navier–Stokes equations used in computational hemodynamics) and compressible flows (e.g., the weakly compressible Euler equations used in weather prediction).

Our review revolves around the mathematical tools that are used to construct LES-ROMs, which not only provide sound support for LES-ROMs but also connect (“bridge”) two important research fields, LES and ROMs, that have been investigated separately. Our presentation focused on spatial filters, both at the FOM level and at the ROM level, and approximate deconvolution (AD) operators, which provide efficient and stable approximations to the inverse of the spatial filter. The FOM spatial filters we discussed are the differential filter and the nonlinear filter. We used these filters to build LES models. The ROM spatial filters we presented are the ROM differential filter, the ROM higher-order algebraic filter, and the ROM projection. We then illustrated how these ROM spatial filters are leveraged to construct LES-ROMs. As for the AD, we presented several strategies (e.g., Lavrentiev, Tikhonov, and Van Cittert) and showed how they can be used to increase the accuracy of both LES models and LES-ROMs.

Equipped with these mathematical tools (i.e., spatial filtering and AD), we then showed how they can be used to construct LES models and LES-ROMs. For clarity of presentation, we focused on one of the most popular and successful models, the evolve–filter–relax (EFR) strategy, and discussed it at both the FOM and the ROM levels. At the FOM level, we showed how the EFR model is constructed for both the incompressible Navier–Stokes equations and the weakly compressible Euler equations. Then, we illustrated some of the EFR’s achievements in incompressible flows (computational hemodynamics) and compressible flows (the rising thermal bubble and the density current benchmarks). At the ROM level, we constructed the EFR-ROM and several other LES-ROMs (the Leray ROM, AD Leray ROM, and time-relaxation ROM). Then, we illustrated some LES-ROMs’ achievements for both incompressible flows (the classical flow past a cylinder and the T-junction test case with thermal striping) and compressible flows (again, the rising thermal bubble benchmark and density currents).

We are aware that this is not an exhaustive review of LES-ROMs. However, in writing this paper, our hope was to provide a brief and friendly introduction to this exciting research area, which we believe has a lot of potential in practical numerical simulation of convection-dominated flows in engineering, science, and medicine. We hope we conveyed to the reader the reason for our excitement: LES-ROMs are extremely easy to implement and use, highly efficient, and accurate in reproducing average flow quantities of interest, even in the case of challenging, critical applications of convection-dominated flows.

Despite the LES-ROMs’ success, there are still many research directions that need to be further explored. This is natural since the LES-ROMs are much more recent than classical LES models. Among the promising research directions in LES-ROM development, we mention the construction, analysis, and investigation of other LES-ROMs, spatial filters, and approximate deconvolution operators. Indeed, at the FOM level, there is a plethora of LES models that are constructed with various types of spatial filters and closure modeling strategies. A few of those have been already used to build LES-ROMs, but there are many more that have not been investigated. Another interesting research direction is the construction of new, data-driven spatial scales (lengthscales) at the ROM level. Indeed, at the FOM level, the spatial filter radius can be defined based on the mesh used in the numerical discretization. At the ROM level, however, how should we define the filter radius? In this paper, we outlined a few possibilities, but much more remains to be undertaken. Probably one of the most important research directions is the investigation of the new LES-ROMs in classical test problems for numerical methods for turbulent flows, e.g., the turbulent channel flow [305]. These test problems will yield new insight into the practical application of the proposed LES-ROMs in the numerical simulation of turbulent flows. We note, however, that the investigations of ROMs for turbulent flows is relatively scarce (see, e.g., [182] for steps in this direction). Finally, an important research direction is continuing to lay the mathematical foundations of LES-ROMs. Only the first steps in this direction have been made, as we briefly discussed in this paper. There is, however, much more that needs to be done to bring the mathematical support of LES-ROMs to the same level as the mathematical support of LES (see, e.g., the research monographs [5,6,29,30] on the mathematical theory of LES).

Author Contributions

Conceptualization, A.Q., O.S., A.V. and T.I.; methodology, A.Q., O.S., A.V. and T.I.; software, A.Q., O.S., A.V. and T.I.; validation, A.Q., O.S., A.V. and T.I.; formal analysis, A.Q., O.S., A.V. and T.I.; investigation, A.Q., O.S., A.V. and T.I.; resources, A.Q., O.S., A.V. and T.I.; data curation, A.Q., O.S., A.V. and T.I.; writing—original draft preparation, A.Q., O.S., A.V. and T.I.; writing—review and editing, A.Q., O.S., A.V. and T.I.; visualization, A.Q., O.S., A.V. and T.I.; supervision, A.Q., O.S., A.V. and T.I.; project administration, A.Q., O.S., A.V. and T.I.; funding acquisition, A.Q., O.S., A.V. and T.I. All authors have read and agreed to the published version of the manuscript.

Funding

A.Q. was supported by the U.S. National Science Foundation through grant number DMS-1953535. O.S. was supported by the U.S. National Science Foundation through grant number DMS-2345048. A.V. was supported by the U.S. National Science Foundation through grant number DMS-2012286 and grant number DMS-2038118, the U.S. National Institutes of Health through grant number R01EB031101, and Emory University Research Committee. T.I. was supported by the U.S. National Science Foundation through grant number DMS-2012253 and grant number CDS&E-MSS-1953113.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

We thank Jiaqi Yang for the simulations of aortic dissections, and Marina Piccinelli and Bradley Leshnower for providing the data and elaborating the images. We also thank Gianluigi Rozza, Michele Girfoglio, Francesco Ballarin, and Maria Strazzullo for carefully reading the manuscript and providing feedback that has improved the quality and clarity of the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Kolmogorov, A.N. The local structure of turbulence in incompressible viscous fluids at very large Reynolds numbers. Dokl. Akad. Nauk. SSSR 1941, 30, 299–303. [Google Scholar]
Kolmogorov, A.N. Dissipation of energy in isotropic turbulence. Dokl. Akad. Nauk. SSSR 1941, 32, 19–21. [Google Scholar]
Duraisamy, K.; Iaccarino, G.; Xiao, H. Turbulence modeling in the age of data. Annu. Rev. Fluid Mech. 2019, 51, 357–377. [Google Scholar] [CrossRef]
Sagaut, P. Large Eddy Simulation for Incompressible Flows, 3rd ed.; Scientific Computation; Springer: Berlin/Heidelberg, Germany, 2006; p. xxx+556. [Google Scholar]
Berselli, L.C.; Iliescu, T.; Layton, W.J. Mathematics of Large Eddy Simulation of Turbulent Flows; Scientific Computation; Springer: Berlin/Heidelberg, Germany, 2006; p. xviii+348. [Google Scholar]
Layton, W.J.; Rebholz, L.G. Approximate Deconvolution Models of Turbulence: Analysis, Phenomenology and Numerical Analysis; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012; Volume 2042. [Google Scholar]
Delorme, Y.; Anupindi, K.; Kerlo, A.; Shetty, D.; Rodefeld, M.; Chen, J.; Frankel, S. Large eddy simulation of powered Fontan hemodynamics. J. Biomech. 2013, 46, 408–422. [Google Scholar] [CrossRef] [PubMed]
Xu, H.; Piccinelli, M.; Leshnower, B.G.; Lefieux, A.; Taylor, W.R.; Veneziani, A. Coupled morphological–hemodynamic computational analysis of type B aortic dissection: A longitudinal study. Ann. Biomed. Eng. 2018, 46, 927–939. [Google Scholar] [CrossRef]
Manchester, E.L.; Pirola, S.; Salmasi, M.Y.; O’Regan, D.P.; Athanasiou, T.; Xu, X.Y. Analysis of turbulence effects in a patient-specific aorta with aortic valve stenosis. Cardiovasc. Eng. Technol. 2021, 12, 438–453. [Google Scholar] [CrossRef] [PubMed]
Caldwell, P.M.; Mametjanov, A.; Tang, Q.; Van Roekel, L.P.; Golaz, J.C.; Lin, W.; Bader, D.C.; Keen, N.D.; Feng, Y.; Jacob, R.; et al. The DOE E3SM coupled model version 1: Description and results at high resolution. J. Adv. Model. Earth Syst. 2019, 11, 4095–4146. [Google Scholar] [CrossRef]
Terai, C.; Caldwell, P.; Klein, S.; Tang, Q.; Branstetter, M. The atmospheric hydrologic cycle in the ACME v0.3 model. Clim. Dyn. 2018, 50, 3251–3279. [Google Scholar] [CrossRef]
Wehner, M.; Reed, K.; Li, F.; Prabhat; Bacmeister, J.; Chen, C.T.; Paciorek, C.; Gleckler, P.J.; Sperber, K.R.; Collins, W.D.; et al. The effect of horizontal resolution on simulation quality in the Community Atmospheric Model, CAM5.1. J. Adv. Model. Earth Syst. 2014, 6, 980–997. [Google Scholar] [CrossRef]
Bacmeister, J.; Wehner, M.; Neale, R.; Gettelman, A.; Hannay, C.; Lauritzen, P.; Caron, J.M.; Truesdale, J.E. Exploratory high-resolution climate simulations using the Community Atmosphere Model (CAM). J. Clim. 2014, 27, 3073–3099. [Google Scholar] [CrossRef]
Delworth, T.; Rosati, A.; Anderson, W.; Adcroft, A.J.; Balaji, V.; Benson, R.; Dixon, K.; Griffies, S.M.; Lee, H.C.; Pacanowski, R.C.; et al. Simulated climate and climate change in the GFDL CM2.5 high-resolution coupled climate model. J. Clim. 2012, 25, 2755–2781. [Google Scholar] [CrossRef]
Love, B.; Matthews, A.; Lister, G. The diurnal cycle of precipitation over the maritime continent in a high-resolution atmospheric model. Quaterly J. R. Meteorol. Soc. 2011, 137, 934–947. [Google Scholar] [CrossRef]
Atlas, R.; Reale, O.; Shen, B.W.; Lin, S.J.; Chern, J.D.; Putman, W.; Lee, T.; Yeh, K.S.; Bosilovich, M.; Radakovich, J. Hurricane forecasting with the high-resolution NASA finite volume general circulation model. Geophys. Res. Lett. 2005, 32. [Google Scholar] [CrossRef]
Iorio, J.; Duffy, P.; Govindasamy, B.; Thompson, S.; Khairoutdinov, M.; Randall, D. Effects of model resolution and subgrid-scale physics on the simulation of precipitation in the continental United States. Clim. Dyn. 2004, 23, 243–258. [Google Scholar] [CrossRef]
Duffy, P.; Govindasamy, B.; Iorio, J.; Milanovich, J.; Sperber, K.; Taylor, K.; Wehner, M.F.; Thompson, S.L. High-resolution simulations of global climate, Part 1: Present climate. Clim. Dyn. 2003, 21, 371–390. [Google Scholar] [CrossRef]
Pope, V.; Stratton, R. The processes governing horizontal resolution sensitivity in a climate model. Clim. Dyn. 2002, 19, 211–236. [Google Scholar]
Leonard, A. Energy Cascade in Large-Eddy Simulations of Turbulent Fluid Flows. In Turbulent Diffusion in Environmental Pollution; Frenkiel, F., Munn, R., Eds.; Advances in Geophysics; Elsevier: Amsterdam, The Netherlands, 1975; Volume 18, pp. 237–248. [Google Scholar] [CrossRef]
Germano, M. Turbulence: The filtering approach. J. Fluid Mech. 1992, 238, 325–336. [Google Scholar] [CrossRef]
Moser, R.D.; Haering, S.W.; Yalla, G.R. Statistical Properties of Subgrid-Scale Turbulence Models. Annu. Rev. Fluid Mech. 2021, 53, 255–286. [Google Scholar] [CrossRef]
Beck, A.; Flad, D.; Munz, C.D. Deep neural networks for data-driven LES closure models. J. Comput. Phys. 2019, 398, 108910. [Google Scholar] [CrossRef]
Sirignano, J.; MacArt, J.F.; Freund, J.B. DPM: A deep learning PDE augmentation method with application to large-eddy simulation. J. Comput. Phys. 2020, 423, 109811. [Google Scholar] [CrossRef]
Xie, C.; Wang, J.; E, W. Modeling subgrid-scale forces by spatial artificial neural networks in large eddy simulation of turbulence. Phys. Rev. Fluids 2020, 5, 054606. [Google Scholar] [CrossRef]
Duraisamy, K. Perspectives on machine learning-augmented Reynolds-averaged and large eddy simulation models of turbulence. Phys. Rev. Fluids 2021, 6, 050504. [Google Scholar] [CrossRef]
Raissi, M.; Yazdani, A.; Karniadakis, G.E. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science 2020, 367, 1026–1030. [Google Scholar] [CrossRef] [PubMed]
Di Leoni, P.C.; Zaki, T.A.; Karniadakis, G.; Meneveau, C. Two-point stress–strain-rate correlation structure and non-local eddy viscosity in turbulent flows. J. Fluid Mech. 2021, 914, A6. [Google Scholar] [CrossRef]
John, V. Large Eddy Simulation of Turbulent Incompressible Flows. In Lecture Notes in Computational Science and Engineering; Springer: Berlin/Heidelberg, Germany, 2004; Volume 34, p. xii+261. [Google Scholar]
Rebollo, T.C.; Lewandowski, R. Mathematical and Numerical Foundations of Turbulence Models and Applications; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
Frisch, U. Turbulence: The Legacy of A.N. Kolmogorov; Cambridge University Press: Cambridge, UK, 1995. [Google Scholar]
Tennekes, H.; Lumley, J. A First Course in Turbulence; MIT Press: Cambridge, MA, USA, 1972. [Google Scholar]
Smagorinsky, J. General Circulation Experiments with the Primitive Equations: I. The basic experiement. Mon. Wea. Rev. 1963, 91, 99–164. [Google Scholar] [CrossRef]
Abgrall, R. Toward the Ultimate Conservative Scheme: Following the Quest. J. Comput. Phys. 2001, 167, 277–315. [Google Scholar] [CrossRef]
Prandtl, L. Turbulent Flow. In Lecture Delivered before the International Congress for Applied Mechanics; US Department of Commerce: Springfield, VA, USA, 1926. Available online: https://ntrs.nasa.gov/api/citations/19930090799/downloads/19930090799.pdf (accessed on 28 June 2024).
Kloeckner, A.; Warburton, T.; Hesthaven, J.S. Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method. Math. Model. Nat. Phenom. 2011, 6, 57–83. [Google Scholar] [CrossRef]
Rispoli, F.; Saavedra, R. A stabilized finite element method based on SGS models for compressible flows. Comp. Meth. Appl. Mech. Eng. 2006, 196, 652–664. [Google Scholar] [CrossRef]
Persson, P.O.; Peraire, J. Sub-cell shock capturing for discontinuous Galerkin methods. In Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 9–12 January 2006; p. 112. [Google Scholar]
Guermond, J.L.; Pasquetti, R.; Popov, B. Entropy viscosity method for nonlinear conservation laws. J. Comput. Phys. 2011, 230, 4248–4267. [Google Scholar] [CrossRef]
Guermond, J.L.; Pasquetti, R. Entropy-based nonlinear viscosity for Fourier approximations of conservation laws. C. R. Acad. Sci. Ser. I 2008, 346, 801–806. [Google Scholar] [CrossRef]
Guermond, J.L.; Popov, B. Viscous regularization of the Euler equations and entropy principles. SIAM J. Appl. Math. 2014, 74, 284–305. [Google Scholar] [CrossRef]
Kurganov, A.; Liu, Y. New adaptive artificial viscosity method for hyperbolic systems of conservation laws. J. Comput. Phys. 2012, 231, 8114–8132. [Google Scholar] [CrossRef]
Marras, S.; Nazarov, M.; Giraldo, F.X. Stabilized high-order Galerkin methods based on a parameter-free dynamic SGS model for LES. J. Comput. Phys. 2015, 301, 77–101. [Google Scholar] [CrossRef]
Wang, Z.; Triantafyllou, M.S.; Constantinides, Y.; Karniadakis, G.E. An entropy-viscosity large eddy simulation study of turbulent flow in a flexible pipe. J. Fluid Mech. 2019, 859, 691–730. [Google Scholar] [CrossRef]
Bazilevs, Y.; Calo, V.; Cottrell, J.A.; Hughes, T.J.R.; Reali, A.; Scovazzi, G. Variational Multiscale Residual-based Turbulence Modeling for Large Eddy Simulation of Incompressible Flows. Comput. Methods Appl. Mech. Eng. 2007, 197, 173–201. [Google Scholar] [CrossRef]
Codina, R. Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. Comput. Methods Appl. Mech. Eng. 2002, 191, 4295–4321. [Google Scholar] [CrossRef]
Codina, R.; Badia, S.; Baiges, J.; Principe, J. Variational Multiscale Methods in Computational Fluid Dynamics. In Encyclopedia of Computational Mechanics, 2nd ed.; John Wiley & Sons, Ltd: Hoboken, NJ, USA, 2017; pp. 1–28. [Google Scholar]
Hughes, T.J.R.; Feijóo, G.; Mazzei, L.; Quincy, J. The variational multiscale method—A paradigm for computational mechanics. Comput. Methods Appl. Mech. Eng. 1998, 166, 3–24. [Google Scholar] [CrossRef]
Guermond, J.L.; Pasquetti, R.; Popov, B. From suitable weak solutions to entropy viscosity. J. Sci. Comput. 2011, 49, 35–50. [Google Scholar] [CrossRef]
Olshanskii, M.; Xiong, X. A connection between filter stabilization and eddy viscosity models. Numer. Methods Partial. Differ. Equations 2013, 29, 2061–2080. [Google Scholar] [CrossRef]
Clinco, N.; Girfoglio, M.; Quaini, A.; Rozza, G. Filter stabilization for the mildly compressible Euler equations with application to atmosphere dynamics simulations. Comput. Fluids 2023, 266, 106057. [Google Scholar] [CrossRef]
Boyd, J.P. Two Comments on Filtering (Artificial Viscosity) for Chebyshev and Legendre Spectral and Spectral Element Methods: Preserving Boundary Conditions and Interpretation of the Filter as a Diffusion. J. Comput. Phys. 1998, 143, 283–288. [Google Scholar] [CrossRef]
Fischer, P.; Mullen, J. Filter-based stabilization of spectral element methods. Comptes Rendus De L’académie Des Sci.-Ser. I-Math. 2001, 332, 265–270. [Google Scholar] [CrossRef]
Mullen, J.S.; Fischer, P.F. Filtering techniques for complex geometry fluid flows. Commun. Numer. Meth. Eng. 1999, 15, 9–18. [Google Scholar] [CrossRef]
Mathew, J.; Lechner, R.; Foysi, H.; Sesterhenn, J.; Friedrich, R. An explicit filtering method for large eddy simulation of compressible flows. Phys. Fluids 2003, 15, 2279–2289. [Google Scholar] [CrossRef]
Visbal, M.; Rizzetta, D. Large eddy simulation on curvilinear grids using compact differencing and filtering schemes. J. Fluids Eng. 2002, 124, 836–847. [Google Scholar] [CrossRef]
Garnier, E.; Adams, N.; Sagaut, P. Large Eddy Simulation for Compressible Flows; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
Girfoglio, M.; Quaini, A.; Rozza, G. A Finite Volume approximation of the Navier-Stokes equations with nonlinear filtering stabilization. Comput. Fluids 2019, 187, 27–45. [Google Scholar] [CrossRef]
Bertagna, L.; Quaini, A.; Veneziani, A. Deconvolution-based nonlinear filtering for incompressible flows at moderately large Reynolds numbers. Int. J. Num. Meth. Fluids 2016, 81, 463–488. [Google Scholar] [CrossRef]
Bowers, A.L.; Rebholz, L.G.; Takhirov, A.; Trenchea, C. Improved accuracy in regularization models of incompressible flow via adaptive nonlinear filtering. Int. J. Numer. Methods Fluids 2012, 70, 805–828. [Google Scholar] [CrossRef]
Bowers, A.; Rebholz, L. Numerical study of a regularization model for incompressible flow with deconvolution-based adaptive nonlinear filtering. Comput. Methods Appl. Mech. Eng. 2013, 258, 1–12. [Google Scholar] [CrossRef]
Layton, W.; Rebholz, L.G.; Trenchea, C. Modular nonlinear filter stabilization of methods for higher Reynolds numbers flow. J. Math. Fluid Mech. 2012, 14, 325–354. [Google Scholar] [CrossRef]
Layton, W.; Röhe, L.; Tran, H. Explicitly uncoupled VMS stabilization of fluid flow. Comput. Methods Appl. Mech. Eng. 2011, 200, 3183–3199. [Google Scholar] [CrossRef]
Viguerie, A.; Veneziani, A. Deconvolution-based stabilization of the incompressible Navier–Stokes equations. J. Comput. Phys. 2019, 391, 226–242. [Google Scholar] [CrossRef]
Brezis, H. Functional Analysis, Sobolev Spaces and Partial Differential Equations; Springer: New York, NY, USA, 2011. [Google Scholar]
Girfoglio, M.; Quaini, A.; Rozza, G. A novel Large Eddy Simulation model for the Quasi-Geostrophic equations in a Finite Volume setting. J. Comput. Appl. Math. 2023, 418, 114656. [Google Scholar] [CrossRef]
Besabe, L.; Girfoglio, M.; Quaini, A.; Rozza, G. Linear and nonlinear filtering for a two-layer quasi-geostrophic ocean model. arXiv 2024, arXiv:2404.11718. [Google Scholar]
Quarteroni, A.; Valli, A. Domain Decomposition Methods for Partial Differential Equations; Oxford Science Publications: Oxford, UK, 1999. [Google Scholar]
Deparis, S.; Fernandez, M.A.; Formaggia, L. Acceleration of a fixed point algorithm for fluid-structure interaction using transpiration conditions. ESAIM Math. Model. Numer. Anal. 2003, 37, 601–616. [Google Scholar] [CrossRef]
Ervin, V.; Layton, W.; Neda, M. Numerical Analysis of Filter-Based Stabilization for Evolution Equations. SIAM J. Numer. Anal. 2012, 50, 2307–2335. [Google Scholar] [CrossRef]
Girfoglio, M.; Quaini, A.; Rozza, G. A POD-Galerkin reduced order model for a LES filtering approach. J. Comput. Phys. 2021, 436, 110260. [Google Scholar] [CrossRef]
Leray, J. Sur le mouvement d‘un fluide visqueux emplissant l’espace. Acta Math. 1934, 63, 193–248. [Google Scholar] [CrossRef]
Geurts, B.J.; Holm, D.D. Regularization modeling for large-eddy simulation. Phys. Fluids 2003, 15, L13–L16. [Google Scholar] [CrossRef]
Guermond, J.L.; Oden, J.T.; Prudhomme, S. Mathematical perspectives on large eddy simulation models for turbulent flows. J. Math. Fluid Mech. 2004, 6, 194–248. [Google Scholar] [CrossRef]
Foiaş, C.; Holm, D.; Titi, E. The Navier-Stokes-alpha model of fluid turbulence. Phys. D 2001, 152/153, 505–519. [Google Scholar] [CrossRef]
Chehab, J.P. Damping, Stabilization and Numerical Filtering for the Modeling and the Simulation of time dependent PDEs. Discret. Contin. Dyn. Syst. Ser. S 2021, 14, 2693–2728. [Google Scholar] [CrossRef]
Hesthaven, J.S.; Warburton, T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications; Springer Publishing Company, Incorporated: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
Holm, D.D. Averaged Lagrangians and the mean effects of fluctuations in ideal fluid dynamics. Phys. D Nonlinear Phenom. 2002, 170, 253–286. [Google Scholar] [CrossRef]
Secchi, P. An alpha model for compressible fluids. Discrete Contin. Dyn. Syst. S. 2008, 3, 351–359. [Google Scholar] [CrossRef]
Giraldo, F.X.; Restelli, M. A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases. J. Comput. Phys. 2008, 227, 3849–3877. [Google Scholar] [CrossRef]
Marras, S.; Moragues, M.; Vázquez, M.; Jorba, O.; Houzeaux, G. A Variational Multiscale Stabilized Finite Element Method for the Solution of the Euler Equations of Nonhydrostatic Stratified Flows. J. Comput. Phys. 2013, 236, 380–407. [Google Scholar] [CrossRef]
Marras, S.; Kelly, J.; Moragues, M.; Müller, A.; Kopera, M.; Vázquez, M.; Giraldo, F.; Houzeaux, G.; Jorba, O. A Review of Element-Based Galerkin Methods for Numerical Weather Prediction: Finite Elements, Spectral Elements, and Discontinuous Galerkin. Arch. Comput. Methods Eng. 2016, 23, 673–722. [Google Scholar] [CrossRef]
Marchuk, G.I. Numerical Methods in Weather Prediction; Academic Press: Cambridge, MA, USA, 1974; pp. 1–277. [Google Scholar]
Borggaard, J.; Iliescu, T.; Roop, J.P. A bounded artificial viscosity large eddy simulation model. SIAM J. Num. Anal. 2009, 47, 622–645. [Google Scholar] [CrossRef]
Hunt, J.; Wray, A.; Moin, P. Eddies Stream and Convergence Zones in Turbulent Flows; Technical Report CTR-S88, CTR Report; NASA: Washington, DC, USA, 1988.
Vreman, A. An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications. Phys. Fluids 2004, 16, 3670–3681. [Google Scholar] [CrossRef]
Stolz, S.; Adams, N. An approximate deconvolution procedure for large-eddy simulation. Phys. Fluids 1999, 11, 1699–1701. [Google Scholar] [CrossRef]
Stolz, S.; Adams, N.; Kleiser, L. An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded flows. Phys. Fluids 2001, 13, 997–1015. [Google Scholar] [CrossRef]
Dunca, A.; Epshteyn, Y. On the Stolz-Adams deconvolution model for the large-eddy simulation of turbulent flows. SIAM J. Math. Anal. 2005, 37, 1890–1902. [Google Scholar] [CrossRef]
Brunton, S.L.; Noack, B.R.; Koumoutsakos, P. Machine learning for fluid mechanics. Annu. Rev. Fluid Mech. 2020, 52, 477–508. [Google Scholar] [CrossRef]
Ling, J.; Kurzawski, A.; Templeton, J. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. J. Fluid Mech. 2016, 807, 155–166. [Google Scholar] [CrossRef]
Beck, A.; Kurz, M. A perspective on machine learning methods in turbulence modeling. GAMM-Mitteilungen 2021, 44, e202100002. [Google Scholar] [CrossRef]
Sirignano, J.; MacArt, J.F. Deep learning closure models for large-eddy simulation of flows around bluff bodies. J. Fluid Mech. 2023, 966, A26. [Google Scholar] [CrossRef]
Schmelzer, M.; Dwight, R.P.; Cinnella, P. Discovery of algebraic Reynolds-stress models using sparse symbolic regression. Flow Turbul. Combust. 2020, 104, 579–603. [Google Scholar] [CrossRef]
Reissmann, M.; Hasslberger, J.; Sandberg, R.D.; Klein, M. Application of gene expression programming to a-posteriori LES modeling of a Taylor Green vortex. J. Comput. Phys. 2021, 424, 109859. [Google Scholar] [CrossRef]
Wu, C.; Zhang, Y. Enhancing the shear-stress-transport turbulence model with symbolic regression: A generalizable and interpretable data-driven approach. Phys. Rev. Fluids 2023, 8, 084604. [Google Scholar] [CrossRef]
Ren, P.; Erichson, N.B.; Subramanian, S.; San, O.; Lukic, Z.; Mahoney, M.W. SuperBench: A Super-Resolution Benchmark Dataset for Scientific Machine Learning. arXiv 2023, arXiv:2306.14070. [Google Scholar]
Maulik, R.; San, O.; Jacob, J.D.; Crick, C. Sub-grid scale model classification and blending through deep learning. J. Fluid Mech. 2019, 870, 784–812. [Google Scholar] [CrossRef]
Maulik, R.; San, O.; Jacob, J.D. Spatiotemporally dynamic implicit large eddy simulation using machine learning classifiers. Phys. Nonlinear Phenom. 2020, 406, 132409. [Google Scholar] [CrossRef]
Kim, B.; Azevedo, V.C.; Thuerey, N.; Kim, T.; Gross, M.; Solenthaler, B. Deep fluids: A generative network for parameterized fluid simulations. In Computer Graphics Forum; Wiley Online Library: Hoboken, NJ, USA, 2019; Volume 38, pp. 59–70. [Google Scholar]
Kontolati, K.; Goswami, S.; Karniadakis, G.E.; Shields, M.D. Learning in latent spaces improves the predictive accuracy of deep neural operators. arXiv 2023, arXiv:2304.07599. [Google Scholar]
Du, P.; Parikh, M.H.; Fan, X.; Liu, X.Y.; Wang, J.X. CoNFiLD: Conditional Neural Field Latent Diffusion Model Generating Spatiotemporal Turbulence. arXiv 2024, arXiv:2403.05940. [Google Scholar]
Lozano-Durán, A.; Bae, H.J. Machine learning building-block-flow wall model for large-eddy simulation. J. Fluid Mech. 2023, 963, A35. [Google Scholar] [CrossRef]
Vadrot, A.; Yang, X.I.; Abkar, M. Survey of machine-learning wall models for large-eddy simulation. Phys. Rev. Fluids 2023, 8, 064603. [Google Scholar] [CrossRef]
Veneziani, A.; Vergara, C. Inverse problems in cardiovascular mathematics: Toward patient-specific data assimilation and optimization. Int. J. Numer. Methods Biomed. Eng. 2013, 29, 723–725. [Google Scholar] [CrossRef] [PubMed]
Taylor, C.A.; Fonte, T.A.; Min, J.K. Computational fluid dynamics applied to cardiac computed tomography for noninvasive quantification of fractional flow reserve: Scientific basis. J. Am. Coll. Cardiol. 2013, 61, 2233–2241. [Google Scholar] [CrossRef]
Taylor, C.A.; Petersen, K.; Xiao, N.; Sinclair, M.; Bai, Y.; Lynch, S.R.; UpdePac, A.; Schaap, M. Patient-specific modeling of blood flow in the coronary arteries. Comput. Methods Appl. Mech. Eng. 2023, 417, 116414. [Google Scholar] [CrossRef]
Weiss, G.; Wolner, I.; Folkmann, S.; Sodeck, G.; Schmidli, J.; Grabenwöger, M.; Carrel, T.; Czerny, M. The location of the primary entry tear in acute type B aortic dissection affects early outcome. Eur. J.-Cardio-Thorac. Surg. 2012, 42, 571–576. [Google Scholar] [CrossRef]
Izzo, R.; Steinman, D.; Manini, S.; Antiga, L. The vascular modeling toolkit: A python library for the analysis of tubular structures in medical images. J. Open Source Softw. 2018, 3, 745. [Google Scholar] [CrossRef]
Xu, H. Efficient Modeling of the Incompressible Flow with Moderate Large Reynolds Numbers Using a Deconvolution-Based Leray Model: Analysis, Uncertainty Quantification and Application in Aortic Dissections. Ph.D. Thesis, Georgia Institute of Technology, Atlanta, GA, USA, 2020. [Google Scholar]
Bäumler, K.; Vedula, V.; Sailer, A.M.; Seo, J.; Chiu, P.; Mistelbauer, G.; Chan, F.P.; Fischbein, M.P.; Marsden, A.L.; Fleischmann, D. Fluid–structure interaction simulations of patient-specific aortic dissection. Biomech. Model. Mechanobiol. 2020, 19, 1607–1628. [Google Scholar] [CrossRef] [PubMed]
Zimmermann, J.; Bäumler, K.; Loecher, M.; Cork, T.E.; Marsden, A.L.; Ennis, D.B.; Fleischmann, D. Hemodynamic effects of entry and exit tear size in aortic dissection evaluated with in vitro magnetic resonance imaging and fluid–structure interaction simulation. Sci. Rep. 2023, 13, 22557. [Google Scholar] [CrossRef]
Girfoglio, M.; Quaini, A.; Rozza, G. Fluid-structure interaction simulations with a LES filtering approach in solids4Foam. Commun. Appl. Ind. Math. 2022, 13, 13–28. [Google Scholar] [CrossRef]
Heywood, J.G.; Rannacher, R.; Turek, S. Artificial boundaries and flux and pressure conditions for the incompressible Navier–Stokes equations. Int. J. Numer. Methods Fluids 1996, 22, 325–352. [Google Scholar] [CrossRef]
Womersley, J.R. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J. Physiol. 1955, 127, 553. [Google Scholar] [CrossRef] [PubMed]
Veneziani, A.; Vergara, C. An approximate method for solving incompressible Navier–Stokes problems with flow rate conditions. Comput. Methods Appl. Mech. Eng. 2007, 196, 1685–1700. [Google Scholar] [CrossRef]
Formaggia, L.; Gerbeau, J.F.; Nobile, F.; Quarteroni, A. Numerical treatment of defective boundary conditions for the Navier–Stokes equations. SIAM J. Numer. Anal. 2002, 40, 376–401. [Google Scholar] [CrossRef]
Veneziani, A.; Vergara, C. Flow rate defective boundary conditions in haemodynamics simulations. Int. J. Numer. Methods Fluids 2005, 47, 803–816. [Google Scholar] [CrossRef]
Formaggia, L.; Veneziani, A.; Vergara, C. A new approach to numerical solution of defective boundary value problems in incompressible fluid dynamics. SIAM J. Numer. Anal. 2008, 46, 2769–2794. [Google Scholar] [CrossRef]
Formaggia, L.; Veneziani, A.; Vergara, C. Flow rate boundary problems for an incompressible fluid in deformable domains: Formulations and solution methods. Comput. Methods Appl. Mech. Eng. 2010, 199, 677–688. [Google Scholar] [CrossRef]
Xu, H.; Baroli, D.; Di Massimo, F.; Quaini, A.; Veneziani, A. Backflow stabilization by deconvolution-based large eddy simulation modeling. J. Comput. Phys. 2020, 404, 109103. [Google Scholar] [CrossRef]
Armour, C.H.; Guo, B.; Pirola, S.; Saitta, S.; Liu, Y.; Dong, Z.; Xu, X.Y. The influence of inlet velocity profile on predicted flow in type B aortic dissection. Biomech. Model. Mechanobiol. 2021, 20, 481–490. [Google Scholar] [CrossRef] [PubMed]
Formaggia, L.; Nobile, F.; Quarteroni, A.; Veneziani, A. Multiscale modelling of the circulatory system: A preliminary analysis. Comput. Vis. Sci. 1999, 2, 75–83. [Google Scholar] [CrossRef]
Vignon-Clementel, I.E.; Figueroa, C.A.; Jansen, K.E.; Taylor, C.A. Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Comput. Methods Appl. Mech. Eng. 2006, 195, 3776–3796. [Google Scholar] [CrossRef]
Formaggia, L.; Quarteroni, A.; Veneziani, A. Multiscale models of the vascular system. Cardiovasc. Math. Model. Simul. Circ. Syst. 2009, 395–446. [Google Scholar]
Quarteroni, A.; Veneziani, A.; Vergara, C. Geometric multiscale modeling of the cardiovascular system, between theory and practice. Comput. Methods Appl. Mech. Eng. 2016, 302, 193–252. [Google Scholar] [CrossRef]
Peiró, J.; Veneziani, A. Reduced models of the cardiovascular system. In Cardiovascular Mathematics: Modeling and sImulation of the Circulatory System; Springer: Berlin/Heidelberg, Germany, 2009; pp. 347–394. [Google Scholar]
Romarowski, R.M.; Lefieux, A.; Morganti, S.; Veneziani, A.; Auricchio, F. Patient-specific CFD modelling in the thoracic aorta with PC-MRI–based boundary conditions: A least-square three-element Windkessel approach. Int. J. Numer. Methods Biomed. Eng. 2018, 34, e3134. [Google Scholar] [CrossRef] [PubMed]
Pirola, S.; Guo, B.; Menichini, C.; Saitta, S.; Fu, W.; Dong, Z.; Xu, X.Y. 4-D flow MRI-based computational analysis of blood flow in patient-specific aortic dissection. IEEE Trans. Biomed. Eng. 2019, 66, 3411–3419. [Google Scholar] [CrossRef]
Bertoglio, C.; Caiazzo, A.; Bazilevs, Y.; Braack, M.; Esmaily, M.; Gravemeier, V.; L. Marsden, A.; Pironneau, O.; Vignon-Clementel, I.; Wall, W. Benchmark problems for numerical treatment of backflow at open boundaries. Int. J. Numer. Methods Biomed. Eng. 2018, 34, e2918. [Google Scholar] [CrossRef]
Bruneau, C.H.; Fabrie, P. Effective downstream boundary conditions for incompressible Navier–Stokes equations. Int. J. Numer. Methods Fluids 1994, 19, 693–705. [Google Scholar] [CrossRef]
Bruneau, C.H.; Fabrie, P. New efficient boundary conditions for incompressible Navier-Stokes equations: A well-posedness result. ESAIM: Math. Model. Numer. Anal. 1996, 30, 815–840. [Google Scholar] [CrossRef]
Arbia, G.; Vignon-Clementel, I.; Hsia, T.Y.; Gerbeau, J.F.; Modeling of Congenital Hearts Alliance (MOCHA) Investigators. Modified Navier–Stokes equations for the outflow boundary conditions in hemodynamics. Eur. J.-Mech.-B/Fluids 2016, 60, 175–188. [Google Scholar] [CrossRef]
Bertoglio, C.; Caiazzo, A. A tangential regularization method for backflow stabilization in hemodynamics. J. Comput. Phys. 2014, 261, 162–171. [Google Scholar] [CrossRef]
Bertoglio, C.; Caiazzo, A. A Stokes-residual backflow stabilization method applied to physiological flows. J. Comput. Phys. 2016, 313, 260–278. [Google Scholar] [CrossRef]
Nichols, W.W.; O’Rourke, M.; Edelman, E.R.; Vlachopoulos, C. McDonald’s Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles; CRC Press: Boca Raton, FL, USA, 2022. [Google Scholar]
Chen, D.; Müller-Eschner, M.; von Tengg-Kobligk, H.; Barber, D.; Böckler, D.; Hose, R.; Ventikos, Y. A patient-specific study of type-B aortic dissection: Evaluation of true-false lumen blood exchange. Biomed. Eng. Online 2013, 12, 65. [Google Scholar] [CrossRef] [PubMed]
Fatma, K.; Carine, G.C.; Marine, G.; Philippe, P.; Valérie, D. Numerical modeling of residual type B aortic dissection: Longitudinal analysis of favorable and unfavorable evolution. Med. Biol. Eng. Comput. 2022, 60, 769–783. [Google Scholar] [CrossRef] [PubMed]
Bertagna, L.; Deparis, S.; Formaggia, L.; Forti, D.; Veneziani, A. The LifeV library: Engineering mathematics beyond the proof of concept. arXiv 2017, arXiv:1710.06596. [Google Scholar]
Yang, J.; Piccinelli, M.; Leshnower, B.; Veneziani, A. Predicting The Evolution Of Type-b Aortic Dissection: Combining Computational Hemodynamics And Data Analysis. In Proceedings of the 8th International Conference on Computational and Mathematical Biomedical Engineering—CMBE2024, Arlington, VA, USA, 24–26 June 2024. [Google Scholar]
Bertagna, L.; Quaini, A.; Rebholz, L.G.; Veneziani, A. On the sensitivity to the filtering radius in Leray models of incompressible flow. In Contributions to Partial Differential Equations and Applications; Springer: Berlin/Heidelberg, Germany, 2018; pp. 111–130. [Google Scholar]
Xu, H.; Baroli, D.; Veneziani, A. Global sensitivity analysis for patient-specific aortic simulations: The role of geometry, boundary condition and large Eddy simulation modeling parameters. J. Biomech. Eng. 2021, 143, 021012. [Google Scholar] [CrossRef]
Sobol, I.M. Sensitivity estimates for nonlinear mathematical models. Math. Model. Comput. Exp. 1993, 1, 407–414. [Google Scholar]
Sobol, I.M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simul. 2001, 55, 271–280. [Google Scholar] [CrossRef]
Sudret, B. Global sensitivity analysis using polynomial chaos expansions. Reliab. Eng. Syst. Saf. 2008, 93, 964–979. [Google Scholar] [CrossRef]
Crestaux, T.; Le Maıtre, O.; Martinez, J.M. Polynomial chaos expansion for sensitivity analysis. Reliab. Eng. Syst. Saf. 2009, 94, 1161–1172. [Google Scholar] [CrossRef]
D’Elia, M.; Mirabella, L.; Passerini, T.; Perego, M.; Piccinelli, M.; Vergara, C.; Veneziani, A. Applications of variational data assimilation in computational hemodynamics. Model. Physiol. Flows 2012, 363–394. [Google Scholar]
D’Elia, M.; Perego, M.; Veneziani, A. A variational data assimilation procedure for the incompressible Navier-Stokes equations in hemodynamics. J. Sci. Comput. 2012, 52, 340–359. [Google Scholar] [CrossRef]
D’Elia, M.; Veneziani, A. Uncertainty quantification for data assimilation in a steady incompressible Navier-Stokes problem. ESAIM Math. Model. Numer. Anal. 2013, 47, 1037–1057. [Google Scholar] [CrossRef]
Cai, S.; Mao, Z.; Wang, Z.; Yin, M.; Karniadakis, G.E. Physics-informed neural networks (PINNs) for fluid mechanics: A review. Acta Mech. Sin. 2021, 37, 1727–1738. [Google Scholar] [CrossRef]
Ahmad, N.N.; Lindeman, J. Euler solutions using flux-based wave decomposition. Int. J. Numer. Meth. Fluids 2007, 54, 47–72. [Google Scholar] [CrossRef]
Feng, Y.; Miranda-Fuentes, J.; Jacob, J.; Sagaut, P. Hybrid lattice Boltzmann model for atmospheric flows under anelastic approximation. Phys. Fluids 2021, 33, 036607. [Google Scholar] [CrossRef]
Carpenter, R.; Droegemeier, K.; Woodward, P.; Hane, C. Application of the piecewise parabolic method (PPM) to meteorological modeling. Mon. Wea. Rev. 1990, 118, 586–612. [Google Scholar] [CrossRef]
Straka, J.; Wilhelmson, R.; Wicker, L.; Anderson, J.; Droegemeier, K. Numerical solution of a nonlinear density current: A benchmark solution and comparisons. Int. J. Num. Meth. Fluids 1993, 17, 1–22. [Google Scholar] [CrossRef]
Özgökmen, T.; Iliescu, T.; Fischer, P.; Srinivasan, A.; Duan, J. Large Eddy Simulation of Stratified Mixing in Two-Dimensional Dam-Break Problem in a Rectangular Enclosed Domain. Ocean. Model. 2007, 16, 106–140. [Google Scholar] [CrossRef]
Özgökmen, T.; Iliescu, T.; Fischer, P. Large Eddy Simulation of Stratified Mixing in a Three-Dimensional Lock-Exchange System. Ocean. Model. 2009, 26, 134–155. [Google Scholar] [CrossRef]
Özgökmen, T.; Iliescu, T.; Fischer, P. Reynolds number dependence of mixing in a lock-exchange system from direct numerical and large eddy simulations. Ocean. Model. 2009, 30, 190–206. [Google Scholar] [CrossRef]
Ahmad, N.N. High-Resolution Wave Propagation Method for Stratified Flows. In Proceedings of the AIAA Aviation Forum, Atlanta, GA, USA, 25–29 June 2018. [Google Scholar] [CrossRef]
Girfoglio, M.; Quaini, A.; Rozza, G. Validation of an OpenFOAM^®-based solver for the Euler equations with benchmarks for mesoscale atmospheric modeling. AIP Adv. 2023, 13, 055024. [Google Scholar] [CrossRef]
Benner, P.; Schilders, W.; Grivet-Talocia, S.; Quarteroni, A.; Rozza, G.; Miguel Silveira, L. Model Order Reduction: Volume 1: System- and Data-Driven Methods and Algorithms; De Gruyter: Berlin, Germany, 2021. [Google Scholar]
Benner, P.; Schilders, W.; Grivet-Talocia, S.; Quarteroni, A.; Rozza, G.; Miguel Silveira, L. Model Order Reduction: Volume 2: Snapshot-Based Methods and Algorithms; De Gruyter: Berlin, Germany, 2021. [Google Scholar]
Benner, P.; Schilders, W.; Grivet-Talocia, S.; Quarteroni, A.; Rozza, G.; Miguel Silveira, L. Model Order Reduction: Volume 3: Applications; De Gruyter: Berlin, Germany, 2021. [Google Scholar]
Hesthaven, J.S.; Rozza, G.; Stamm, B. Certified Reduced Basis Methods for Parametrized Partial Differential Equations; Springer: Berlin/Heidelberg, Germany, 2016; Volume 590. [Google Scholar]
Malik, M.H. Reduced Order Modeling for Smart Grids’ Simulation and Optimization. Ph.D. Thesis, École Centrale de Nantes. Universitat Politécnica de Catalunya, Barcelona, Spain, 2017. [Google Scholar]
Rozza, G.; Huynh, D.B.P.; Patera, A.T. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics. Arch. Comput. Methods Eng. 2008, 15, 229–275. [Google Scholar] [CrossRef]
Quarteroni, A.; Manzoni, A.; Negri, F. Reduced Basis Methods for Partial Differential Equations: An Introduction; Springer: Berlin/Heidelberg, Germany, 2015; Volume 92. [Google Scholar]
Wang, Z.; Akhtar, I.; Borggaard, J.; Iliescu, T. Proper orthogonal decomposition closure models for turbulent flows: A numerical comparison. Comput. Meth. Appl. Mech. Eng. 2012, 237–240, 10–26. [Google Scholar] [CrossRef]
Aubry, N.; Holmes, P.; Lumley, J.L.; Stone, E. The dynamics of coherent structures in the wall region of a turbulent boundary layer. J. Fluid Mech. 1988, 192, 115–173. [Google Scholar] [CrossRef]
Ahmed, S.E.; Pawar, S.; San, O.; Rasheed, A.; Iliescu, T.; Noack, B.R. On closures for reduced order models—A spectrum of first-principle to machine-learned avenues. Phys. Fluids 2021, 33, 091301. [Google Scholar] [CrossRef]
Carlberg, K.; Bou-Mosleh, C.; Farhat, C. Efficient non-linear model reduction via a least-squares Petrov–Galerkin projection and compressive tensor approximations. Int. J. Num. Meth. Eng. 2011, 86, 155–181. [Google Scholar] [CrossRef]
Grimberg, S.; Farhat, C.; Youkilis, N. On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows. J. Comput. Phys. 2020, 419, 109681. [Google Scholar] [CrossRef]
Baiges, J.; Codina, R.; Idelsohn, S. Explicit Reduced Order Models for the stabilized finite element approximation of the incompressible Navier-Stokes equations. Int. J. Num. Meth. Fluids 2013, 72, 1219–1243. [Google Scholar] [CrossRef]
Reyes, R.; Codina, R. Projection-based reduced order models for flow problems: A variational multiscale approach. Comput. Methods Appl. Mech. Eng. 2020, 363, 112844. [Google Scholar] [CrossRef]
Parish, E.J.; Wentland, C.; Duraisamy, K. The adjoint Petrov–Galerkin method for non-linear model reduction. Comput. Meth. Appl. Mech. Eng. 2020, 365, 112991. [Google Scholar] [CrossRef]
Wells, D.; Wang, Z.; Xie, X.; Iliescu, T. An evolve-then-filter regularized reduced order model for convection-dominated flows. Int. J. Num. Meth. Fluids 2017, 84, 598–615. [Google Scholar] [CrossRef]
Gunzburger, M.; Iliescu, T.; Mohebujjaman, M.; Schneier, M. An evolve-filter-relax stabilized reduced order stochastic collocation method for the time-dependent Navier-Stokes equations. SIAM-ASA J. Uncertain. 2019, 7, 1162–1184. [Google Scholar] [CrossRef]
Sabetghadam, F.; Jafarpour, A. α regularization of the POD-Galerkin dynamical systems of the Kuramoto–Sivashinsky equation. Appl. Math. Comput. 2012, 218, 6012–6026. [Google Scholar] [CrossRef]
Wells, D. Stabilization of POD-ROMs. Ph.D. Thesis, Virginia Tech, Blacksburg, VA, USA, 2015. Available online: http://vtechworks.lib.vt.edu/bitstream/handle/10919/52960/Wells_DR_D_2015.pdf?sequence=1&isAllowed=y (accessed on 28 June 2024).
Germano, M. Differential filters of elliptic type. Phys. Fluids 1986, 29, 1757–1758. [Google Scholar] [CrossRef]
Strazzullo, M.; Girfoglio, M.; Ballarin, F.; Iliescu, T.; Rozza, G. Consistency of the full and reduced order models for evolve-filter-relax regularization of convection-dominated, marginally-resolved flows. Int. J. Numer. Methods Eng. 2022, 123, 3148–3178. [Google Scholar] [CrossRef]
Girfoglio, M.; Quaini, A.; Rozza, G. A linear filter regularization for POD-based reduced-order models of the quasi-geostrophic equations. C. R. Mech. 2023, 351, 1–21. [Google Scholar] [CrossRef]
Tsai, P.H.; Fischer, P.; Iliescu, T. A Time-Relaxation Reduced Order Model for the Turbulent Channel Flow. arXiv 2023, arXiv:2312.13272. [Google Scholar]
Mou, C.; Merzari, E.; San, O.; Iliescu, T. An energy-based lengthscale for reduced order models of turbulent flows. Nucl. Eng. Des. 2023, 412, 112454. [Google Scholar] [CrossRef]
Holmes, P.; Lumley, J.L.; Berkooz, G. Turbulence, Coherent Structures, Dynamical Systems and Symmetry; Cambridge University Press: Cambridge, UK, 1996. [Google Scholar]
Aradag, S.; Siegel, S.; Seidel, J.; Cohen, K.; McLaughlin, T. Filtered POD-based low-dimensional modeling of the 3D turbulent flow behind a circular cylinder. Int. J. Num. Meth. Fluids 2011, 66, 1–16. [Google Scholar] [CrossRef]
Farcas, I.; Munipalli, R.; Willcox, K.E. On filtering in non-intrusive data-driven reduced-order modeling. In Proceedings of the AIAA AVIATION 2022 Forum, Chicago, IL, USA, 27 June–1 July 2022; p. 3487. [Google Scholar]
Bertero, M.; Boccacci, P. Introduction to Inverse Problems in Imaging; Institute of Physics Publishing: Bristol, UK, 1998; p. xii+351. [Google Scholar]
Hansen, P.C. Discrete Inverse Problems: Insight and Algorithms. Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2010; Volume 7. [Google Scholar]
Xie, X.; Wells, D.; Wang, Z.; Iliescu, T. Approximate Deconvolution Reduced Order Modeling. Comput. Methods Appl. Mech. Eng. 2017, 313, 512–534. [Google Scholar] [CrossRef]
Sanfilippo, A.; Moore, I.R.; Ballarin, F.; Iliescu, T. Approximate deconvolution Leray reduced order model. Finite Elem. Anal. Des. 2023, 226, 104021. [Google Scholar] [CrossRef]
Cordier, L.; Abou El Majd, B.; Favier, J. Calibration of POD reduced-order models using Tikhonov regularization. Int. J. Num. Meth. Fluids 2010, 63, 269–296. [Google Scholar] [CrossRef]
Wang, Y.; Navon, I.M.; Wang, X.; Cheng, Y. 2D Burgers equation with large Reynolds number using POD/DEIM and calibration. Int. J. Num. Meth. Fluids 2016, 82, 909–931. [Google Scholar] [CrossRef]
Weller, J.; Lombardi, E.; Iollo, A. Robust model identification of actuated vortex wakes. Phys. D 2009, 238, 416–427. [Google Scholar] [CrossRef]
Strazzullo, M.; Ballarin, F.; Iliescu, T.; Canuto, C. New Feedback Control and Adaptive Evolve-Filter-Relax Regularization for the Navier-Stokes Equations in the Convection-Dominated Regime. arXiv 2023, arXiv:2307.00675. [Google Scholar]
Xie, X.; Wells, D.; Wang, Z.; Iliescu, T. Numerical Analysis of the Leray Reduced Order Model. J. Comput. Appl. Math. 2018, 328, 12–29. [Google Scholar] [CrossRef]
Gunzburger, M.; Iliescu, T.; Schneier, M. A Leray regularized ensemble-proper orthogonal decomposition method for parameterized convection-dominated flows. IMA J. Numer. Anal. 2020, 40, 886–913. [Google Scholar] [CrossRef]
Girfoglio, M.; Quaini, A.; Rozza, G. A hybrid projection/data-driven reduced order model for the Navier-Stokes equations with nonlinear filtering stabilization. J. Comput. Phys. 2023, 486, 112127. [Google Scholar] [CrossRef]
Kaneko, K.; Tsai, P.H.; Fischer, P. Towards model order reduction for fluid-thermal analysis. Nucl. Eng. Des. 2020, 370, 110866. [Google Scholar] [CrossRef]
Tsai, P.H.; Fischer, P. Parametric model-order-reduction development for unsteady convection. Front. Phys. 2022, 10, 903169. [Google Scholar] [CrossRef]
Tsai, P.H. Parametric Model Order Reduction Development for Navier-Stokes Equations from 2D Chaotic to 3D Turbulent Flow Problems. Ph.D. Thesis, University of Illinois at Urbana-Champaign, Champaign, IL, USA, 2023. [Google Scholar]
Sanderse, B.; Stinis, P.; Maulik, R.; Ahmed, S.E. Scientific machine learning for closure models in multiscale problems: A review. arXiv 2024, arXiv:2403.02913. [Google Scholar]
Hughes, T.J.R. Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput. Methods Appl. Mech. Engrg. 1995, 127, 387–401. [Google Scholar] [CrossRef]
Hughes, T.J.R.; Mazzei, L.; Oberai, A.; Wray, A. The multiscale formulation of large eddy simulation: Decay of homogeneous isotropic turbulence. Phys. Fluids 2001, 13, 505–512. [Google Scholar] [CrossRef]
Codina, R.; Reyes, R.; Baiges, J. A posteriori error estimates in a finite element VMS-based reduced order model for the incompressible Navier-Stokes equations. Mech. Res. Commun. 2021, 112, 103599. [Google Scholar] [CrossRef]
Mou, C.; Koc, B.; San, O.; Rebholz, L.G.; Iliescu, T. Data-Driven Variational Multiscale Reduced Order Models. Comput. Methods Appl. Mech. Eng. 2021, 373, 113470. [Google Scholar] [CrossRef]
Ingimarson, S.; Rebholz, L.G.; Iliescu, T. Full and Reduced Order Model Consistency of the Nonlinearity Discretization in Incompressible Flows. Comput. Meth. Appl. Mech. Eng. 2022, 401, 115620. [Google Scholar] [CrossRef]
Giere, S.; Iliescu, T.; John, V.; Wells, D. SUPG Reduced Order Models for Convection-Dominated Convection-Diffusion-Reaction Equations. Comput. Methods Appl. Mech. Eng. 2015, 289, 454–474. [Google Scholar] [CrossRef]
Zoccolan, F.; Strazzullo, M.; Rozza, G. Stabilized weighted reduced order methods for parametrized advection-dominated optimal control problems governed by partial differential equations with random inputs. arXiv 2023, arXiv:2301.01975. [Google Scholar]
Zoccolan, F.; Strazzullo, M.; Rozza, G. A streamline upwind Petrov-Galerkin reduced order method for advection-dominated partial differential equations under optimal control. Comput. Methods Appl. Math. 2024. [Google Scholar] [CrossRef]
Pacciarini, P.; Rozza, G. Stabilized reduced basis method for parametrized advection–diffusion PDEs. Comput. Meth. Appl. Mech. Eng. 2014, 274, 1–18. [Google Scholar] [CrossRef]
Roos, H.G.; Stynes, M.; Tobiska, L. Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems, 2nd ed.; Springer Series in Computational Mathematics; Springer: Berlin/Heidelberg, Germany, 2008; Volume 24. [Google Scholar]
Borggaard, J.; Iliescu, T.; Wang, Z. Artificial Viscosity Proper Orthogonal Decomposition. Math. Comput. Modelling 2011, 53, 269–279. [Google Scholar] [CrossRef]
Iliescu, T.; Wang, Z. Variational multiscale proper orthogonal decomposition: Convection-dominated convection-diffusion-reaction equations. Math. Comput. 2013, 82, 1357–1378. [Google Scholar] [CrossRef]
Iliescu, T.; Wang, Z. Variational Multiscale Proper Orthogonal Decomposition: Navier-Stokes Equations. Num. Meth. P.D.E.s 2014, 30, 641–663. [Google Scholar] [CrossRef]
Ballarin, F.; Rebollo, T.C.; Ávila, E.D.; Mármol, M.G.; Rozza, G. Certified reduced basis VMS-Smagorinsky model for natural convection flow in a cavity with variable height. Comput. Math. Appl. 2020, 80, 973–989. [Google Scholar] [CrossRef]
Rebollo, T.C.; Ávila, E.D.; Mármol, M.G.; Ballarin, F.; Rozza, G. On a certified Smagorinsky reduced basis turbulence model. SIAM J. Numer. Anal. 2017, 55, 3047–3067. [Google Scholar] [CrossRef]
Azaïez, M.; Rebollo, T.C.; Rubino, S. A cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations. J. Comput. Phys. 2021, 425, 109916. [Google Scholar] [CrossRef]
John, V.; Moreau, B.; Novo, J. Error analysis of a SUPG-stabilized POD-ROM method for convection-diffusion-reaction equations. Comput. Math. Appl. 2022, 122, 48–60. [Google Scholar] [CrossRef]
Xie, X.; Mohebujjaman, M.; Rebholz, L.G.; Iliescu, T. Data-Driven Filtered Reduced Order Modeling of Fluid Flows. SIAM J. Sci. Comput. 2018, 40, B834–B857. [Google Scholar] [CrossRef]
Ivagnes, A.; Stabile, G.; Mola, A.; Iliescu, T.; Rozza, G. Pressure data-driven variational multiscale reduced order models. J. Comput. Phys. 2023, 476, 111904. [Google Scholar] [CrossRef]
Ivagnes, A.; Stabile, G.; Mola, A.; Iliescu, T.; Rozza, G. Hybrid data-driven closure strategies for reduced order modeling. Appl. Math. Comput. 2023, 448, 127920. [Google Scholar] [CrossRef]
Koc, B.; Mohebujjaman, M.; Mou, C.; Iliescu, T. Commutation error in reduced order modeling of fluid flows. Adv. Comput. Math. 2019, 45, 2587–2621. [Google Scholar] [CrossRef]
Mohebujjaman, M.; Rebholz, L.G.; Iliescu, T. Physically-constrained data-driven correction for reduced order modeling of fluid flows. Int. J. Num. Meth. Fluids 2019, 89, 103–122. [Google Scholar] [CrossRef]
Mou, C.; Liu, H.; Wells, D.R.; Iliescu, T. Data-Driven Correction Reduced Order Models for the Quasi-Geostrophic Equations: A Numerical Investigation. Int. J. Comput. Fluid Dyn. 2020, 34, 147–159. [Google Scholar] [CrossRef]
Xie, X.; Webster, C.; Iliescu, T. Closure Learning for Nonlinear Model Reduction Using Deep Residual Neural Network. Fluids 2020, 5, 39. [Google Scholar] [CrossRef]
Ahmed, S.E.; San, O.; Rasheed, A.; Iliescu, T.; Veneziani, A. Physics guided machine learning for variational multiscale reduced order modeling. SIAM J. Sci. Comput. 2023, 45, B283–B313. [Google Scholar] [CrossRef]
Koc, B.; Mou, C.; Liu, H.; Wang, Z.; Rozza, G.; Iliescu, T. Verifiability of the Data-Driven Variational Multiscale Reduced Order Model. J. Sci. Comput. 2022, 93, 1–26. [Google Scholar] [CrossRef]
Turek, S.; Schäfer, M. Benchmark computations of laminar flow around cylinder. In Flow Simulation with High-Performance Computers II; Hirschel, E., Ed.; Notes on Numerical Fluid Mechanics; Vieweg: Decatur, IL, USA, 1996; Volume 52. [Google Scholar]
John, V. Reference values for drag and lift of a two dimensional time-dependent flow around a cylinder. Int. J. Numer. Methods Fluids 2004, 44, 777–788. [Google Scholar] [CrossRef]
John, V. On the efficiency of linearization schemes and coupled multigrid methods in the simulation of a 3D flow around a cylinder. Int. J. Numer. Methods Fluids 2006, 50, 845–862. [Google Scholar] [CrossRef]
Bayraktar, E.; Mierka, O.; Turek, S. Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow. Int. J. Comput. Sci. Eng. 2012, 7, 253–266. [Google Scholar] [CrossRef]
Girfoglio, M.; Quaini, A.; Rozza, G. Pressure stabilization strategies for a LES filtering Reduced Order Model. Fluids 2021, 6, 302. [Google Scholar] [CrossRef]
Xie, X.; Bao, F.; Webster, C. Evolve Filter Stabilization Reduced-Order Model for Stochastic Burgers Equation. Fluids 2018, 3, 84. [Google Scholar] [CrossRef]
Stabile, G.; Rozza, G. Finite volume POD-Galerkin stabilised reduced order methods for the parametrised incompressible Navier‚ÄìStokes equations. Comput. Fluids 2018, 173, 273–284. [Google Scholar] [CrossRef]
Akhtar, I.; Nayfeh, A.H.; Ribbens, C.J. On the stability and extension of reduced-order Galerkin models in incompressible flows. Theor. Comp. Fluid Dyn. 2009, 23, 213–237. [Google Scholar] [CrossRef]
Lazzaro, D.; Montefusco, L. Radial Basis Functions for the Multivariate Interpolation of Large Scattered Data Sets. J. Comput. Appl. Math. 2002, 140, 521–536. [Google Scholar] [CrossRef]
Barrault, M.; Nguyen, N.C.; Maday, Y.; Patera, A.T. An “Empirical Interpolation” Method: Application to Efficient Reduced-Basis Discretization of Partial Differential Equations. Comptes Rendus Math. 2004, 339, 667–672. [Google Scholar] [CrossRef]
Chaturantabut, S.; Sorensen, D. Nonlinear Model Reduction via Discrete Empirical Interpolation. SIAM J. Sci. Comput. 2010, 32, 2737–2764. [Google Scholar] [CrossRef]
Lorenzi, S.; Cammi, A.; Luzzi, L.; Rozza, G. POD-Galerkin method for finite volume approximation of Navier-Stokes and RANS equations. Comput. Methods Appl. Mech. Eng. 2016, 311, 151–179. [Google Scholar] [CrossRef]
Hijazi, S.; Stabile, G.; Mola, A.; Rozza, G. Data-Driven POD-Galerkin reduced order model for turbulent flows. J. Comput. Phys. 2020, 416, 109513. [Google Scholar] [CrossRef]
Obabko, A.V.; Fischer, P.F.; Tautges, T.J.; Karabasov, S.; Goloviznin, V.M.; Zaytsev, M.A.; Chudanov, V.V.; Pervichko, V.A.; Aksenova, A.E. CFD Validation in OECD/NEA T-Junction Benchmark; Technical Report; Argonne National Lab. (ANL): Argonne, IL, USA, 2011. [Google Scholar]
de Leval, M.R.; Kilner, P.; Gewillig, M.; Bull, C.; McGoon, D.C. Total cavopulmonary connection: A logical alternative to atriopulmonary connection for complex Fontan operations: Experimental studies and early clinical experience. J. Thorac. Cardiovasc. Surg. 1988, 96, 682–695. [Google Scholar] [CrossRef]
Sharma, S.; Goudy, S.; Walker, P.; Panchal, S.; Ensley, A.; Kanter, K.; Tam, V.; Fyfe, D.; Yoganathan, A. In vitro flow experiments for determination of optimal geometry of total cavopulmonary connection for surgical repair of children with functional single ventricle. J. Am. Coll. Cardiol. 1996, 27, 1264–1269. [Google Scholar] [CrossRef]
Migliavacca, F.; Kilner, P.J.; Pennati, G.; Dubini, G.; Pietrabissa, R.; Fumero, R.; de Leval, M.R. Computational fluid dynamic and magnetic resonance analyses of flow distribution between the lungs after total cavopulmonary connection. IEEE Trans. Biomed. Eng. 1999, 46, 393–399. [Google Scholar] [CrossRef] [PubMed]
Ensley, A.E.; Lynch, P.; Chatzimavroudis, G.P.; Lucas, C.; Sharma, S.; Yoganathan, A.P. Toward designing the optimal total cavopulmonary connection: An in vitro study. Ann. Thorac. Surg. 1999, 68, 1384–1390. [Google Scholar] [CrossRef] [PubMed]
Dubini, G.; Migliavacca, F.; Pennati, G.; De Leval, M.R.; Bove, E.L. Ten years of modelling to achieve haemodynamic optimisation of the total cavopulmonary connection. Cardiol Young 2004, 14, 48–52. [Google Scholar] [CrossRef] [PubMed]
Hsia, T.Y.; Figliola, R.; Modeling of Congenital Hearts Alliance (MOCHA) Investigators; Bove, E.; Dorfman, A.; Taylor, A.; Giardini, A.; Khambadkone, S.; Schievano, S.; Hsia, T.Y.; et al. Multiscale modelling of single-ventricle hearts for clinical decision support: A Leducq Transatlantic Network of Excellence. Eur. J.-Cardio-Thorac. Surg. 2016, 49, 365–368. [Google Scholar] [CrossRef]
Schiavazzi, D.E.; Kung, E.O.; Marsden, A.L.; Baker, C.; Pennati, G.; Hsia, T.Y.; Hlavacek, A.; Dorfman, A.L.; Modeling of Congenital Hearts Alliance (MOCHA) Investigators. Hemodynamic effects of left pulmonary artery stenosis after superior cavopulmonary connection: A patient-specific multiscale modeling study. J. Thorac. Cardiovasc. Surg. 2015, 149, 689–696. [Google Scholar] [CrossRef]
Pekkan, K.; Kitajima, H.D.; De Zelicourt, D.; Forbess, J.M.; Parks, W.J.; Fogel, M.A.; Sharma, S.; Kanter, K.R.; Frakes, D.; Yoganathan, A.P. Total cavopulmonary connection flow with functional left pulmonary artery stenosis: Angioplasty and fenestration in vitro. Circulation 2005, 112, 3264–3271. [Google Scholar] [CrossRef]
Tang, E.; Restrepo, M.; Haggerty, C.M.; Mirabella, L.; Bethel, J.; Whitehead, K.K.; Fogel, M.A.; Yoganathan, A.P. Geometric characterization of patient-specific total cavopulmonary connections and its relationship to hemodynamics. JACC Cardiovasc. Imaging 2014, 7, 215–224. [Google Scholar] [CrossRef] [PubMed]
Khiabani, R.H.; Restrepo, M.; Tang, E.; De Zélicourt, D.; Sotiropoulos, F.; Fogel, M.; Yoganathan, A.P. Effect of flow pulsatility on modeling the hemodynamics in the total cavopulmonary connection. J. Biomech. 2012, 45, 2376–2381. [Google Scholar] [CrossRef] [PubMed]
Dasi, L.P.; KrishnankuttyRema, R.; Kitajima, H.D.; Pekkan, K.; Sundareswaran, K.S.; Fogel, M.; Sharma, S.; Whitehead, K.; Kanter, K.; Yoganathan, A.P. Fontan hemodynamics: Importance of pulmonary artery diameter. J. Thorac. Cardiovasc. Surg. 2009, 137, 560–564. [Google Scholar] [CrossRef] [PubMed]
Rodefeld, M.D.; Boyd, J.H.; Myers, C.D.; LaLone, B.J.; Bezruczko, A.J.; Potter, A.W.; Brown, J.W. Cavopulmonary assist: Circulatory support for the univentricular Fontan circulation. Ann. Thorac. Surg. 2003, 76, 1911–1916. [Google Scholar] [CrossRef] [PubMed]
Sarfare, S.; Ali, M.S.; Palazzolo, A.; Rodefeld, M.; Conover, T.; Figliola, R.; Giridharan, G.; Wampler, R.; Bennett, E.; Ivashchenko, A. Computational Fluid Dynamics Turbulence Model and Experimental Study for a Fontan Cavopulmonary Assist Device. J. Biomech. Eng. 2023, 145, 111008. [Google Scholar] [CrossRef] [PubMed]
Meneveau, C. Big wind power: Seven questions for turbulence research. J. Turbul. 2019, 20, 2–20. [Google Scholar] [CrossRef]
Antonini, E.G.; Caldeira, K. Spatial constraints in large-scale expansion of wind power plants. Proc. Natl. Acad. Sci. USA 2021, 118, e2103875118. [Google Scholar] [CrossRef] [PubMed]
Bempedelis, N.; Laizet, S.; Deskos, G. Turbulent entrainment in finite-length wind farms. J. Fluid Mech. 2023, 955, A12. [Google Scholar] [CrossRef]
Howland, M.F.; Quesada, J.B.; Martinez, J.J.P.; Larrañaga, F.P.; Yadav, N.; Chawla, J.S.; Sivaram, V.; Dabiri, J.O. Collective wind farm operation based on a predictive model increases utility-scale energy production. Nature Energy 2022, 7, 818–827. [Google Scholar] [CrossRef]
Pawar, S.; Sharma, A.; Vijayakumar, G.; Bay, C.J.; Yellapantula, S.; San, O. Towards multi-fidelity deep learning of wind turbine wakes. Renew. Energy 2022, 200, 867–879. [Google Scholar] [CrossRef]
Luzzatto-Fegiz, P.; Caulfield, C.c.P. Entrainment model for fully-developed wind farms: Effects of atmospheric stability and an ideal limit for wind farm performance. Phys. Rev. Fluids 2018, 3, 093802. [Google Scholar] [CrossRef]
Sedaghatizadeh, N.; Arjomandi, M.; Kelso, R.; Cazzolato, B.; Ghayesh, M.H. Modelling of wind turbine wake using large eddy simulation. Renew. Energy 2018, 115, 1166–1176. [Google Scholar] [CrossRef]
Mehta, D.; Van Zuijlen, A.; Koren, B.; Holierhoek, J.; Bijl, H. Large Eddy Simulation of wind farm aerodynamics: A review. J. Wind. Eng. Ind. Aerodyn. 2014, 133, 1–17. [Google Scholar] [CrossRef]
Meneveau, C. The top-down model of wind farm boundary layers and its applications. J. Turbul. 2012, N7. [Google Scholar] [CrossRef]
Porté-Agel, F.; Wu, Y.T.; Lu, H.; Conzemius, R.J. Large-eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms. J. Wind Eng. Ind. Aerod. 2011, 99, 154–168. [Google Scholar] [CrossRef]
Stadtmann, F.; Rasheed, A.; Kvamsdal, T.; Johannessen, K.A.; San, O.; Kölle, K.; Tande, J.O.; Barstad, I.; Benhamou, A.; Brathaug, T.; et al. Digital twins in wind energy: Emerging technologies and industry-informed future directions. IEEE Access 2023, 11, 110762–110795. [Google Scholar] [CrossRef]
Hajisharifi, A.; Girfoglio, M.; Quaini, A.; Rozza, G. A comparison of data-driven reduced order models for the simulation of mesoscale atmospheric flow. Finite Elem. Anal. Des. 2024, 228, 104050. [Google Scholar] [CrossRef]
Lario, A.; Maulik, R.; Schmidt, O.T.; Rozza, G.; Mengaldo, G. Neural-network learning of SPOD latent dynamics. J. Comput. Phys. 2022, 468, 111475. [Google Scholar] [CrossRef]
Pawar, S.; San, O. Equation-Free Surrogate Modeling of Geophysical Flows at the Intersection of Machine Learning and Data Assimilation. J. Adv. Model. Earth Syst. 2022, 14, e2022MS003170. [Google Scholar] [CrossRef]
Schmidt, O.T.; Mengaldo, G.; Balsamo, G.; Wedi, N.P. Spectral empirical orthogonal function analysis of weather and climate data. Mon. Weather. Rev. 2019, 147, 2979–2995. [Google Scholar] [CrossRef]
Chen, N.; Majda, A.J.; Giannakis, D. Predicting the cloud patterns of the Madden-Julian Oscillation through a low-order nonlinear stochastic model. Geophys. Res. Lett. 2014, 41, 5612–5619. [Google Scholar] [CrossRef]
Chen, N.; Majda, A.J.; Sabeerali, C.T.; Ajayamohan, R.S. Predicting monsoon intraseasonal precipitation using a low-order nonlinear stochastic model. J. Clim. 2018, 31, 4403–4427. [Google Scholar] [CrossRef]
Pathak, J.; Hunt, B.; Girvan, M.; Lu, Z.; Ott, E. Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. Phys. Rev. Lett. 2018, 120, 024102. [Google Scholar] [CrossRef] [PubMed]
Rasp, S.; Thuerey, N. Data-driven medium-range weather prediction with a resnet pretrained on climate simulations: A new model for weatherbench. J. Adv. Model. Earth Syst. 2021, 13, e2020MS002405. [Google Scholar] [CrossRef]
Schultz, M.G.; Betancourt, C.; Gong, B.; Kleinert, F.; Langguth, M.; Leufen, L.H.; Mozaffari, A.; Stadtler, S. Can deep learning beat numerical weather prediction? Philos. Trans. R. Soc. A 2021, 379, 20200097. [Google Scholar] [CrossRef] [PubMed]
Weyn, J.A.; Durran, D.R.; Caruana, R. Can machines learn to predict weather? Using deep learning to predict gridded 500-hPa geopotential height from historical weather data. J. Adv. Model. Earth Syst. 2019, 11, 2680–2693. [Google Scholar] [CrossRef]
Lam, R.; Sanchez-Gonzalez, A.; Willson, M.; Wirnsberger, P.; Fortunato, M.; Pritzel, A.; Ravuri, S.; Ewalds, T.; Alet, F.; Eaton-Rosen, Z.; et al. GraphCast: Learning skillful medium-range global weather forecasting. arXiv 2022, arXiv:2212.12794. [Google Scholar]
Pathak, J.; Subramanian, S.; Harrington, P.; Raja, S.; Chattopadhyay, A.; Mardani, M.; Kurth, T.; Hall, D.; Li, Z.; Azizzadenesheli, K.; et al. FourCastNet: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv 2022, arXiv:2202.11214. [Google Scholar]
Bi, K.; Xie, L.; Zhang, H.; Chen, X.; Gu, X.; Tian, Q. Accurate medium-range global weather forecasting with 3D neural network. Nature 2023, 619, 533–538. [Google Scholar] [CrossRef]
Kutz, J.N.; Brunton, S.L.; Brunton, B.W.; Proctor, J.L. Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, USA, 2016. [Google Scholar]
Schmid, P.J. Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 2010, 656, 5–28. [Google Scholar] [CrossRef]
Schmid, P.J.; Li, L.; Juniper, M.P.; Pust, O. Applications of the dynamic mode decomposition. Theor. Comput. Fluid Dyn. 2011, 25, 249–259. [Google Scholar] [CrossRef]
Tu, J.H.; Rowley, C.W.; Luchtenburg, D.M.; Brunton, S.L.; Kutz, J.N. On Dynamic Mode Decomposition: Theory and Applications. J. Comput. Dyn. 2014, 1, 391–421. [Google Scholar] [CrossRef]
Schmid, P.J. Application of the dynamic mode decomposition to experimental data. Exp. Fluids 2011, 50, 1123–1130. [Google Scholar] [CrossRef]
Duke, D.; Honnery, D.; Soria, J. Experimental investigation of nonlinear instabilities in annular liquid sheets. J. Fluid Mech. 2012, 691, 594–604. [Google Scholar] [CrossRef]
Seena, A.; Sung, H.J. Dynamic mode decomposition of turbulent cavity flows for self-sustained oscillations. Int. J. Heat Fluid Flow 2011, 32, 1098–1110. [Google Scholar] [CrossRef]
Arbabi, H.; Mezic, I. Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the Koopman operator. SIAM J. Appl. Dyn. Syst. 2017, 16, 2096–2126. [Google Scholar] [CrossRef]
Curtis, C.W.; Alford-Lago, D.J.; Bollt, E.; Tuma, A. Machine Learning Enhanced Hankel Dynamic-Mode Decomposition. arXiv 2023, arXiv:2303.06289. [Google Scholar] [CrossRef] [PubMed]
Fujii, K.; Takeishi, N.; Kibushi, B.; Kouzaki, M.; Kawahara, Y. Data-driven spectral analysis for coordinative structures in periodic human locomotion. Sci. Rep. 2019, 9, 16755. [Google Scholar] [CrossRef] [PubMed]
Jiang, H.; Chen, J.; Dong, G.; Liu, T.; Chen, G. Study on Hankel matrix-based SVD and its application in rolling element bearing fault diagnosis. Mech. Syst. Signal Process. 2015, 52, 338–359. [Google Scholar] [CrossRef]
Vasconcelos Filho, E.; dos Santos, P.L. A dynamic mode decomposition approach with Hankel blocks to forecast multi-channel temporal series. IEEE Control. Syst. Lett. 2019, 3, 739–744. [Google Scholar] [CrossRef]
Yang, D.; Gao, H.; Cai, G.; Chen, Z.; Wang, L.; Ma, J.; Li, D. Synchronized ambient data-based extraction of interarea modes using Hankel block-enhanced DMD. Int. J. Electr. Power Energy Syst. 2021, 128, 106687. [Google Scholar] [CrossRef]
Frame, P.; Towne, A. Space-time POD and the Hankel matrix. arXiv 2022, arXiv:2206.08995. [Google Scholar] [CrossRef] [PubMed]
Hess, M.W.; Quaini, A.; Rozza, G. A data-driven surrogate modeling approach for time-dependent incompressible Navier-Stokes equations with dynamic mode decomposition and manifold interpolation. Adv. Comput. Math. 2023, 49, 22. [Google Scholar] [CrossRef]
Girfoglio, M.; Ballarin, F.; Infantino, G.; Nicoló, F.; Montalto, A.; Rozza, G.; Scrofani, R.; Comisso, M.; Musumeci, F. Non-intrusive PODI-ROM for patient-specific aortic blood flow in presence of a LVAD device. Med. Eng. Phys. 2022, 107, 103849. [Google Scholar] [CrossRef] [PubMed]
Hajisharifi, A.; Romano’, F.; Girfoglio, M.; Beccari, A.; Bonanni, D.; Rozza, G. A non-intrusive data-driven reduced order model for parametrized CFD-DEM numerical simulations. J. Comput. Phys. 2023, 491, 112355. [Google Scholar] [CrossRef]
Demo, N.; Tezzele, M.; Mola, A.; Rozza, G. An efficient shape parametrisation by free-form deformation enhanced by active subspace for hull hydrodynamic ship design problems in open source environment. In Proceedings of the ISOPE International Ocean and Polar Engineering Conference. ISOPE, Sapporo, Japan, 10–15 June 2018; pp. 565–572. [Google Scholar]
Demo, N.; Tezzele, M.; Gustin, G.; Lavini, G.; Rozza, G. Shape optimization by means of proper orthogonal decomposition and dynamic mode decomposition. In Proceedings of the Technology and Science for the Ships of the Future: NAV 2018: 19th International Conference on Ship & Maritime Research, Trieste, Italy, 20–22 June 2018; pp. 212–219. [Google Scholar]
Ripepi, M.; Verveld, M.J.; Karcher, N.; Franz, T.; Abu-Zurayk, M.; Görtz, S.; Kier, T. Reduced-order models for aerodynamic applications, loads and MDO. CEAS Aeronaut. J. 2018, 9, 171–193. [Google Scholar] [CrossRef]
Gonzalez, F.J.; Balajewicz, M. Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems. arXiv 2018, arXiv:1808.01346. [Google Scholar]
Maulik, R.; Lusch, B.; Balaprakash, P. Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders. Phys. Fluids 2021, 33, 037106. [Google Scholar] [CrossRef]
Mohan, A.; Daniel, D.; Chertkov, M.; Livescu, D. Compressed convolutional LSTM: An efficient deep learning framework to model high fidelity 3D turbulence. arXiv 2019, arXiv:1903.00033. [Google Scholar]
Shi, X.; Chen, Z.; Wang, H.; Yeung, D.Y.; Wong, W.K.; Woo, W.C. Convolutional LSTM network: A machine learning approach for precipitation nowcasting. Adv. Neural Inf. Process. Syst. 2015, 28. [Google Scholar]
RBniCSx—Reduced Order Modelling in FEniCSx. Available online: https://github.com/rbnics/rbnicsx (accessed on 28 June 2024).
libROM—Library for Reduced Order Models. Available online: https://www.librom.net/ (accessed on 28 June 2024).
Moser, D.R.; Kim, J.; Mansour, N.N. Direct numerical simulation of turbulent channel flow up to Re_τ = 590. Phys. Fluids 1999, 11, 943–945. [Google Scholar] [CrossRef]

Figure 1. Overall view of the energy cascade, from injection to dissipation of energy, and associated types of modeling.

Figure 2. LES schematic showing the input flow variable,

u

, that cannot be represented on a given coarse mesh, and the filtered flow variable,

\bar{u}

, that can be accurately represented on the coarse mesh.

Figure 2. LES schematic showing the input flow variable,

u

, that cannot be represented on a given coarse mesh, and the filtered flow variable,

\bar{u}

, that can be accurately represented on the coarse mesh.

Figure 3. Schematic of the concept proposed in [99].

Figure 4. Images of a patient-specific AoD showing the true lumen and the false lumen.

Figure 5. Simulation in a patient-specific AoD. Top left: pressure. Top right: velocity (in cm/s) in the descending aorta and at the entrance of the false lumen. The two bottom panels outline the complexity of the flow induced by the entry tear for the velocity (left) and the wall shear stress (right).

Figure 6. Anatomies of several AoDs, pinpointing the diversity of the possible morphologies. Geometries reconstructed at Emory University with Vascular Modeling ToolKit [110].

Figure 7. EFR simulation of the hemodynamics in a patient-specific AoD: TAWSS in different regions of the false lumen.

Figure 8. Sobol’ indexes in a patient-specific geometry for the sensitivity of the TAWSS and the OSI to the radius

α

(left), the inflow rate Q (center), and the geometry (right). Blue regions identify parts of the domain weakly affected by variations in the input in comparison with the other uncertainties.

Figure 8. Sobol’ indexes in a patient-specific geometry for the sensitivity of the TAWSS and the OSI to the radius

α

Figure 9. Rising thermal bubble: perturbation of potential temperature

θ^{'}

t = 1020

s computed with the EFR and

a_{S}

and

a_{D}

with the coarser mesh (first two panels) and the finer mesh (last two panels).

Figure 9. Rising thermal bubble: perturbation of potential temperature

θ^{'}

t = 1020

s computed with the EFR and

a_{S}

and

a_{D}

with the coarser mesh (first two panels) and the finer mesh (last two panels).

Figure 10. Density potential temperature fluctuation

θ^{'}

(left) and indicator function (right) for the EFR method with

a_{S}

(top) and

a_{D}

(bottom). The mesh size is

h = 50

Figure 10. Density potential temperature fluctuation

θ^{'}

(left) and indicator function (right) for the EFR method with

a_{S}

(top) and

a_{D}

(bottom). The mesh size is

h = 50

Figure 11. Lift coefficient

C_{L}

computed by the FOM and the projection/data-driven ROM from [197] for different thresholds of cumulative energy.

Figure 11. Lift coefficient

C_{L}

computed by the FOM and the projection/data-driven ROM from [197] for different thresholds of cumulative energy.

Figure 12. Pareto plots for the velocity and pressure: time-averaged relative

L^{2}

error versus relative wall time when the number of basis functions for the velocity is varied for the 2D (left) and 3D (right) cylinder tests.

Figure 12. Pareto plots for the velocity and pressure: time-averaged relative

L^{2}

error versus relative wall time when the number of basis functions for the velocity is varied for the 2D (left) and 3D (right) cylinder tests.

Figure 13. T-junction test case: instantaneous temperature field at

R e

= 10,000 for an investigation on thermal striping, which is critical in nuclear engineering [200].

Figure 13. T-junction test case: instantaneous temperature field at

R e

= 10,000 for an investigation on thermal striping, which is critical in nuclear engineering [200].

Figure 14. Near-wall temperature history at

y = 0

z = 0.45

for several x locations in the outlet branch.

Figure 14. Near-wall temperature history at

y = 0

z = 0.45

for several x locations in the outlet branch.

Figure 15. T-junction at

R e = 10, 000

: comparison of the near-wall temperature history at

z = 0.45

between the FOM, the G-ROM, and the LES-ROMs for

x = 2, 3, 4, 5

(outlet branch).

Figure 15. T-junction at

R e = 10, 000

: comparison of the near-wall temperature history at

z = 0.45

between the FOM, the G-ROM, and the LES-ROMs for

x = 2, 3, 4, 5

(outlet branch).

Figure 16. T-junction at

R e = 10, 000

: comparison of the near-wall temperature history at

z = 0.45

between the FOM, the G-ROM, and the LES-ROMs for

x = 6, 7, 8, 9

(outlet branch).

Figure 16. T-junction at

R e = 10, 000

: comparison of the near-wall temperature history at

z = 0.45

between the FOM, the G-ROM, and the LES-ROMs for

x = 6, 7, 8, 9

(outlet branch).

Figure 17. Simplified geometry of a TCPC. The vertical vessel is the vena cava (VC: superior at the top—SVC, inferior at the bottom—IVC). The pulmonary artery (PA) is the horizontal vessel. The inflow sections are at the SVC and at the IVC. This generates colliding fronts. The picture reports the results corresponding to two different surgical options. The difference is in the flow distribution from the IVC (the so-called hepatic flow distribution):

F D_{L P A}

is the fraction of hepatic flow directed to the left PA. An even flow distribution (i.e.,

F D_{L P A} \approx 50 %

) is desirable. Notice that the different

F D_{L P A}

are created by different offsets between the SVC and IVC.

F D_{L P A}

is the fraction of hepatic flow directed to the left PA. An even flow distribution (i.e.,

F D_{L P A} \approx 50 %

) is desirable. Notice that the different

F D_{L P A}

are created by different offsets between the SVC and IVC.

Figure 18. Rising thermal bubble:

θ^{'}

given by the ROMs and the FOM at time values within (left) and outside (right) the training dataset.

Figure 18. Rising thermal bubble:

θ^{'}

given by the ROMs and the FOM at time values within (left) and outside (right) the training dataset.

Figure 19. Density current:

θ^{'}

given by the ROMs and the FOM at time values within (top) and outside (bottom) the training dataset.

Figure 19. Density current:

θ^{'}

given by the ROMs and the FOM at time values within (top) and outside (bottom) the training dataset.

Table 1. Minimal value of the relaxation parameter that prevents backflows for different values of the radius and the time step for two meshes with two different refinements, coarse (C) and fine (F). See [121] for more details.

$α = h_{m i n}$			$α = 0.9$
Mesh	C	F	Mesh	C	F
$Δ t = \{\begin{matrix} 0.01 \\ 0.005 \end{matrix}$	$\begin{matrix} 0.07 \\ 0.04 \end{matrix}$	$\begin{matrix} 0.08 \\ 0.06 \end{matrix}$	$Δ t = \{\begin{matrix} 0.01 \\ 0.005 \end{matrix}$	$\begin{matrix} 0.06 \\ 0.04 \end{matrix}$	$\begin{matrix} 0.06 \\ 0.05 \end{matrix}$

Table 2. Density current: computational time taken by the evolve step and filter step per time step and total simulation time for the EFR algorithm with indicator functions

a_{S}

and

a_{D}

for mesh

h = 50

Table 2. Density current: computational time taken by the evolve step and filter step per time step and total simulation time for the EFR algorithm with indicator functions

a_{S}

and

a_{D}

for mesh

h = 50

Model	Evolve (s)	Filter (s)	Total (s)
EFR, $a_{S}$	0.06	0.023	772
EFR, $a_{D}$	0.06	0.029	790

Table 3. Rising thermal bubble: computational time needed to construct the reduced basis offline and to perform a simulation online for DMD, HDMD, and PODI.

Method	Basis Construction	Online Run
DMD	0.085 s	0.02 s
HDMD	2.865 s	2.95 s
PODI	0.1 s	0.018 s

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Quaini, A.; San, O.; Veneziani, A.; Iliescu, T. Bridging Large Eddy Simulation and Reduced-Order Modeling of Convection-Dominated Flows through Spatial Filtering: Review and Perspectives. Fluids 2024, 9, 178. https://doi.org/10.3390/fluids9080178

AMA Style

Quaini A, San O, Veneziani A, Iliescu T. Bridging Large Eddy Simulation and Reduced-Order Modeling of Convection-Dominated Flows through Spatial Filtering: Review and Perspectives. Fluids. 2024; 9(8):178. https://doi.org/10.3390/fluids9080178

Chicago/Turabian Style

Quaini, Annalisa, Omer San, Alessandro Veneziani, and Traian Iliescu. 2024. "Bridging Large Eddy Simulation and Reduced-Order Modeling of Convection-Dominated Flows through Spatial Filtering: Review and Perspectives" Fluids 9, no. 8: 178. https://doi.org/10.3390/fluids9080178

APA Style

Quaini, A., San, O., Veneziani, A., & Iliescu, T. (2024). Bridging Large Eddy Simulation and Reduced-Order Modeling of Convection-Dominated Flows through Spatial Filtering: Review and Perspectives. Fluids, 9(8), 178. https://doi.org/10.3390/fluids9080178

Article Menu

Bridging Large Eddy Simulation and Reduced-Order Modeling of Convection-Dominated Flows through Spatial Filtering: Review and Perspectives

Abstract

1. Introduction

2. LES as a Full-Order Model

2.1. Nonlinear Spatial Filtering for LES

2.1.1. The EFR Method for the Incompressible Navier–Stokes Equations

2.1.2. The EFR Method for the Weakly Compressible Euler Equations

2.1.3. Indicator Function

Physical Phenomenology-Based Indicator Functions

Mathematics-Based Indicator Functions

2.2. Machine Learning for LES

3. Applications of LES for FOM

3.1. Incompressible Flows

3.1.1. Computational Hemodynamics in Type B Aortic Dissections

Geometry

Boundary Conditions and Backflows

EFR in Action

Filter Radius and Relaxation Parameter

Results

Sensitivity Analysis for the Filtering Radius

3.1.2. Open Problems

3.2. Compressible Flows

Open Problems

4. LES for Reduced-Order Models

4.1. LES-ROMs

4.1.1. ROM Filters and Approximate Deconvolution

ROM Differential Filter

ROM Higher-Order Algebraic Filter

ROM Projection

ROM Filter Radius

Approximate Deconvolution

4.1.2. EFR-ROM

4.1.3. Leray ROM, Approximate Deconvolution Leray ROM, and Time-Relaxation ROM

Leray ROM

Approximate Deconvolution Leray ROM

Time-Relaxation ROM

4.1.4. Other LES-ROMs

4.2. LES-ROM Consistency

4.3. LES-ROM Numerical Analysis

5. Applications of LES for ROM

5.1. Incompressible Flows

5.1.1. Flow Past a Cylinder

5.1.2. T-Junction

5.1.3. Hemodynamics Applications

5.1.4. Wind Energy Applications

5.2. Compressible Flows

6. Concluding Remarks and Outlook

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI